Robust vehicle routing problem with hard time windows under demand and travel time uncertainty

: Due to an increase in customer-oriented service strategies designed to meet more complex and exacting customer requirements, meeting a scheduled time window has become an important part of designing vehicle routes for logistics activities. However, practically, the uncertainty in travel times and customer demand often means vehicles miss these time windows, increasing service costs and decreasing customer satisfaction. In an eﬀort to ﬁnd a solution that meets the needs of real-world logistics, we examine the vehicle routing problem with hard time windows under demand and travel time uncertainty. To address the problem, we build a robust optimization model based on novel route-dependent uncertainty sets. However, due to the complex nature of the problem, the robust model is only able to tackle small-sized instances using standard solvers. Therefore, to tackle large instances, we design a two-stage algorithm based on a modiﬁed adaptive variable neighborhood search heuristic. The ﬁrst stage of the algorithm minimizes the total number of vehicle routes, while the second stage minimizes the total travel distance. Extensive computational experiments are conducted with modiﬁed versions of Solomon’s benchmark instances. The numerical results show that the proposed two-stage algorithm is able to ﬁnd optimal solutions for small-sized instances and good-quality robust solutions for large-sized instances with little increase to the total travel distance and/or the number of vehicles used. A detailed analysis of the results also reveals several managerial insights for decision-makers in the logistics industry.


Introduction
The vehicle routing problem (VRP) is a combinatorial optimization problem, first introduced by Dantzig and Ramser (1959), that aims to find the optimal set of routes for a fleet of vehicles delivering goods or services to a given set of customers. Due to its broad application in a variety of practical contexts (e.g., logistics, internet routing, and freight transportation), the VRP and its numerous variants have been studied extensively (Baldacci et al., 2012;Kumar and Panneerselvam, 2012;Toth and Vigo, 2014).
The vehicle routing problem with time windows (VRPTW), as an important variant of the VRP, assumes that each customer must be served within a given time window. For example, some customers might only accept deliveries between 11:00 am and 1:00 pm. Due to increasing and more exacting customer requirements, incorporating customer time windows into today's logistic and freight transportation activities has become more and more popular over the last two decades, and the VRPTW and its variants are likely to continue to be a hot research topic among researchers for the foreseeable future (Bräysy and Gendreau, 2005a,b;Baldacci et al., 2012;Miranda and Conceição, 2016). Several popular real-world applications of the VRPTW include waste collection, postal deliveries, school bus routing, and security patrol services (Bräysy and Gendreau, 2005a). The VRPTW assumes that all the input data, such as customer demands and travel times, are both deterministic and are known in advance. However, the solutions derived by deterministic models are often infeasible when applied to real-world situations because the level of customer demand and the travel times are uncertain (Taş et al., 2013). For example, in school bus routing problems, a single bus collects students from several pre-assigned stops and must arrive at the school within a specified time window. However, not all students use the bus on a daily basis and the likelihood that a student will take the bus on any given day varies significantly in practice (Caceres et al., 2017). Thus, the number of students at each stop is highly uncertain. Moreover, the levels of traffic congestion in urban areas means the travel time between any two stops can vary greatly. Hence, ignoring the uncertainty of travel times or student demand may mean students are late for school or that the bus does not have room for all the students. Real-world waste collection problems are also plagued by these two types of uncertainty. The amount of waste at each collection point can vary daily and vehicle speeds are influenced by the prevailing traffic conditions. Indeed, uncertainty in demand and travel times are common issues in many real-world VRPTW applications. And when vehicles fail to serve customers within an agreed upon time window, a reputation for poor service and customer dissatisfaction is often the result.
In an attempt to find an efficient routing strategy for such real-world applications, we study the VRPTW with demand and travel time uncertainty. The robust optimization theory outlined by Ben-Tal et al. (2009) is used to deal with the uncertainties. Robust optimization treats the uncertain parameters of an optimization problem as random variables. Early studies on robust optimization focused on deriving a robust counterpart for the original optimization problem using predefined uncertainty sets, such as ellipsoids or polyhedrons, to obtain a feasible least-cost solution for any realization of the uncertain parameters over the uncertainty sets (Ben-Tal and Nemirovski, 2002;Sim, 2003, 2004). However, more recent works on robust optimization have developed several solution frameworks that can exploit the additional distribution information available in the uncertain parameters to produce less conservative solutions (Bertsimas et al., 2011). One advantage of robust optimization is that the consistency between the tractability of the robust counterpart and the original problem can be maintained as long as the uncertainty sets are designed properly, e.g., by using polyhedral uncertainty sets (Bertsimas and Sim, 2004). Additionally, the uncertainty sets can be derived with only partial information about the distributions of the uncertain parameters, such as supports and moments (Han et al., 2013). Thus, robust optimization is advantageous because of its computational tractability and its ability to tackle real-world applications where only partial distribution information or a small amount of historical data about the uncertain parameters is available.
Our aim through this research is to explore the joint impact of demand and travel time uncertainty on route planning, and produce a cost-effective, robust a priori routing strategy that is insensitive to uncertainty. The resulting contributions of this paper are as follows: (1) a robust optimization model based on novel route-dependent uncertainty sets for the VRPTW with demand and travel time uncertainty; (2) a two-stage algorithm based on a modified adaptive variable neighborhood search (AVNS) that is able to deal with large-sized instances; and (3) an extensive computational study using adapted benchmark instances that shows our two-stage algorithm can generate optimal solutions for small-sized instances and goodquality robust solutions for large-sized instances with little increase to the total travel distance and/or the number of vehicles used.
The remainder of this paper is structured as follows. The relevant literature is reviewed in Section 2. Section 3 describes the VRPTW and a robust version of the VRPTW with demand and travel time uncertainty. The corresponding mathematical models for these two problems are also provided. In Section 4, we describe the details of the proposed two-stage algorithm. Section 5 presents an extensive computational study using the modified versions of Solomon's benchmark instances (Solomon, 1987) and a comprehensive analysis of the results. Section 6 concludes the paper and offers several managerial insights.

Related work
The VRPTW is one of several well-known extensions to the VRP. The literature related to this study includes the VRPTW and its two important variants: stochastic versions of the VRPTW and robust versions of the VRPTW. For a comprehensive background on the VRP, we refer readers to Toth and Vigo (2014).
In the VRPTW, each customer must be served by one vehicle within a scheduled time window. Customer time windows can be "hard", where earliness or tardiness are not allowed, or "soft", where earliness and tardiness are allowed but with a penalty (Miranda and Conceição, 2016). Over the last two decades, scholars have intensively proposed methods of solving the VRPTW with both exact and heuristic algorithms. The most successful exact algorithms rely on a column generation technique to compute the lower bound and a branch-and-price or a branch-and-cut-and-price algorithm to find an optimal integer solution (Jepsen et al., 2008;Desaulniers et al., 2008;Baldacci et al., 2011). Some of the more well-known heuristics include Tabu-search (TS) (Cordeau et al., 2001), large neighborhood search (LNS) (Pisinger and Ropke, 2007), variable neighborhood search (VNS) (Mladenović and Hansen, 1997), iterated local search (Ibaraki et al., 2005) and evolutionary algorithms (Mester and Bräysy, 2005). Detailed surveys of the various methods for solving the VRPTW are provided in Bräysy and Gendreau (2005a,b). As extensions to the VRPTW, a number of more complex routing problems have been explored. One interesting extension is the VRPTW with multiple depots and heterogeneous vehicles. In this scenario, the vehicles have different capacities and multiple vehicle depots are available (Dondo and Cerdá, 2007).
The pickup and delivery problem with time windows is another interesting extension. Here, a number of customer requests, that each involves shipping some goods from a pickup point to a delivery point, need to be satisfied by a fleet of vehicles. Two successful heuristics for solving this problem are the reactive TS heuristic (Nanry and Barnes, 2000) and the adaptive LNS (ALNS) heuristic (Ropke and Pisinger, 2006).
More recently, researchers have begun to associate issues with data uncertainty and the VRPTW, giving rise to stochastic versions of the VRPTW (Ritzinger et al., 2016;Gendreau et al., 2016) and robust versions of the VRPTW (Agra et al., 2012(Agra et al., , 2013Braaten et al., 2017). Reviews of these two important extensions follow.
Stochastic versions of the VRPTW extend the VRPTW by assuming that the key problem parameters, such as customer demands, travel times, and service times, are stochastic. Over the last 20 years, many researchers have studied different stochastic versions of the VRPTW with the majority focusing on stochastic demands and stochastic travel times. Of the studies that focus on the VRPTW with stochastic demands, several different mathematical models and algorithms have been proposed and explored: Ong et al. (1997) proposed a chanceconstrained model; Chang (2005) proposed a two-stage stochastic programming with recourse model and developed a heuristic based on the integer L-shaped method; Lei et al. (2011) also formulated a two-stage stochastic recourse model and devised an ALNS heuristic; and Zhang et al. (2016) proposed three probabilistic models to address on-time deliveries from different perspectives of the carrier and customers. The VRPTW with stochastic travel times was introduced more recently and has attracted greater researcher attention than the VRPTW with stochastic demands ever since. Ando and Taniguchi (2006) used a genetic algorithm to solve the VRPTW with stochastic travel times. Russell and Urban (2008) addressed the VRPTW with stochastic travel times based on the assumption that the travel time between any two customers obeys an Erlang distribution. Li et al. (2010) considered the VRPTW with both stochastic travel times and service times and introduced a chanceconstrained model and a stochastic recourse model. Taş et al. (2013) addressed the VRP with soft time windows and stochastic travel times by assuming the travel times obey a gamma distribution. Taş et al. (2014) also studied the VRP with soft time windows and considered time-dependent and stochastic travel times. A TS heuristic or an ALNS heuristic is used to solve the problem depending on whether or not service times are a factor. More recently, Ehmke et al. (2015) solved the VRPTW with stochastic travel times by a chanceconstraint approach. Recent surveys of stochastic versions of the VRPTW can be found in Ritzinger et al. (2016) and Gendreau et al. (2016).
Our literature review of stochastic versions of the VRPTW reveals that researchers often solve these problems using a two-stage stochastic programming with recourse model or a chance-constrained programming model. However, both these methods have two distinct disadvantages. First, both methods assume that the precise distributions of the uncertain parameters are known or can be estimated in advance, which is rare in real-world VRPTW applications. Second, since both methods suffer from the curse of dimensionality, their computational tractability is seriously affected by the dimensions of the problem at hand.
As a result, both stochastic programming and chance-constrained programming are suitable for small-and medium-sized VRPTW applications where there is enough historical data to fairly accurately estimate the distributions of the uncertain parameters.
Consequently, a few researchers have turned to robust optimization as a new framework for overcoming these disadvantages. Robust optimization treats uncertain parameters as random variables that take their values from predefined uncertainty sets. As such, this framework has great merit in terms of its computational tractability and its ability to solve real-world VRPTW applications with only partial distribution information or little historical data about the uncertain parameters (Gounaris et al., 2014). For a detailed introduction to robust optimization, we refer readers to Ben-Tal et al. (2009).
Unlike stochastic versions of the VRPTW, robust versions of the VRPTW, and even robust versions of the VRP that do not consider time windows, have received much less research attention. According to our review, Sungur et al. (2008) were the first to discuss a robust version of the VRP with demand uncertainty. Gounaris et al. (2013) also considered a robust version of the VRP under demand uncertainty. They developed robust counterparts of several known deterministic VRP formulations and proposed a branch-and-cut solution procedure that is able to optimally solve instances of up to 50 nodes (i.e., customers). Gounaris et al. (2014) presented an adaptive memory programming meta-heuristic that can provide high-quality solutions for the same problem discussed in Gounaris et al. (2013). Researchers have also used robust optimization to deal with some variants of the VRP that consider other types of uncertainty. Ordóñez (2010) showed that the Miller-Tucker-Zemlin model could be modified to handle various robust versions of the VRP that include uncertain travel times, travel costs, and customer demands. Lee et al. (2012) tackled a robust version of the VRP with deadlines under demand and travel time uncertainty using an exact solution algorithm based on a column generation method. However, their algorithm is limited to solving instances of up to 25 nodes and takes around an hour to provide a routing strategy. techniques that can significantly reduce the extreme points of an uncertainty set. However, their exact algorithm is only able to solve instances of up to 40 nodes. Jaillet et al. (2016) proposed a novel mathematical framework based on distributionally robust optimization to address a specific class of routing problems with uncertainty. To illustrate the efficacy of their framework, they solved the VRP with soft time windows given uncertain travel times.

Solano
Heuristic methods are a more recent introduction to solving robust versions of the VRPTW. Toklu et al. (2014) proposed an ant colony optimization heuristic to solve a robust version of the VRPTW under travel time uncertainty. Braaten et al. (2017) developed an ALNS heuristic to solve the same problem discussed in Agra et al. (2013).
Our review of the relevant literature only reveals a limited number of papers that deal with robust versions of the VRPTW. Most only consider one kind of uncertainty, either travel time or customer demand; few consider multiple kinds of uncertainty or explore the joint impact of multiple uncertain parameters on route planning. Moreover, the existing exact algorithms can only tackle robust versions of the VRPTW with small-to-medium-sized instances. Efficient algorithms that are able to solve real-world scenarios with large-sized instances are still needed. Thus, in this paper, we consider a robust version of the VRPTW with both demand and travel time uncertainty and develop a two-stage algorithm that can solve the problem with large-sized instances.
3 Problem description and model formulation

The deterministic VRPTW
The VRPTW is defined on a complete digraph G = (N, A), where the set of nodes is represented by N = {0, 1, · · · , n, n+1}, and the set of arcs is represented by A = {(i, j)|i, j ∈ N, i = j}. N c = {1, · · · , n} denotes the set of customers and K = {1, · · · , m} denotes the set of homogeneous vehicles with a capacity of Q. The vehicle depot is represented by two nodes 0 and n+1. Every vehicle starts from the depot node 0, visits a subset of customer nodes, and ends its route at depot node n+1. Each customer node i ∈ N c has a non-negative demand q i , a service time s i , and a time window [a i , b i ]. a i and b i denote the earliest and latest allowable arrival time at customer node i, respectively. If a vehicle arrives at customer node i before a i , it must wait until a i , and the time window is missed if it arrives at node i after b i . A time window is also associated with the depot nodes 0 and n + 1, i.e., [a 0 , b 0 ] = [a n+1 , b n+1 ]. The travel time between node i and node j is denoted as t ij and, similarly, the travel distance between these nodes is denoted as d ij , where (i, j) ∈ A. Each customer can only be served by exactly one vehicle and the vehicle fleet is big enough to serve all customers. The primary objective of the problem is to minimize the number of vehicle routes; the secondary objective is to minimize the total travel distance.
The mathematical formulation for the VRPTW contains two types of variables. The first are flow variables x ijk for all (i, j) ∈ A and k ∈ K. Let x ijk be 1 if arc (i, j) is used by vehicle k, and 0 otherwise. The second type of variables are time variables y i for all i ∈ N . y i specifies the arrival time of a vehicle at node i. The (deterministic) VRPTW is formulated as follows: (1) k∈K j∈N Objective (1) lexicographically minimizes the number of vehicles (routes) and the total travel distance. Constraints (2) guarantee that each customer can only be visited by one vehicle.
Constraints (3)-(5) are the flow conservation constraints, which ensure that each vehicle's route starts from depot 0 and ends at depot n + 1. Constraints (6) are the vehicle capacity constraints. Constraints (7) calculate the arrival time of a vehicle at node i (i ∈ N ). Constraints (8) ensure that each time window is respected. Constraints (9) deal with the nature of the variables.

A robust VRPTW under demand and travel time uncertainty
As discussed in Section 1, real-world VRPTW applications are often subject to a variety of uncertainties. Therefore, a deterministic VRPTW model, which ignores the uncertainty in data, may not be an appropriate choice; a robust optimization model, which does consider uncertainty, may be more suitable. We use the robust optimization theory discussed in Ben-Tal et al. (2009) and consider a robust version of the VRPTW with both demand and travel time uncertainty. However, to effectively represent the uncertain parameters, robust optimization requires a careful and practical definition of the uncertainty sets. In real-world VRPTW applications, each vehicle route services a different subset of customers. Thus, each vehicle suffers different levels of uncertainty. For example, if a vehicle has a longer pre-planned route and/or more customers on its route, it is likely to experience higher uncertainty. Further, not every customer on a route will have the same level of demand uncertainty, and not every route segment will have the same level of travel time uncertainty.
In reality, a route may only contain a few customers with high demand uncertainty and a few crowed route segments with high travel time uncertainty. Lee et al. (2012) made similar observations and, hence, restricted the number of customers or segments with high uncertainty on each route to control the overall robustness of the routing strategy.
Based on the above observations and discussions, we define two types of uncertainty sets for each vehicle used k: the customer demand uncertainty set U k q and the travel time uncertainty set U k t . Each is defined as a budget uncertainty polytope, as discussed in Bertsimas and Sim (2003,2004). The following polytopes formally define the overall demand uncertainty set U q and travel time uncertainty set U t for the considered robust version of the VRPTW.
Equation (10) reflects that the overall demand uncertainty set U q is the Cartesian product of the demand uncertainty set U k q for each vehicle. In equation (11), N k denotes the set of customers on the route of vehicle k; q i represents the nominal value of uncertain demand q i ; andq i denotes the maximum deviation from the nominal value for each node i ∈ N k . α i is the auxiliary variable, and Γ k q is the uncertainty budget that controls the level of demand uncertainty. Given the level of demand uncertainty experienced by each vehicle relates to the number of customers on its route, Γ k q is defined as equaling θ q |N k | . θ q |N k | is the least integer that is greater than or equal to θ q |N k |. θ q is the demand uncertainty budget coefficient, and has a value of between 0 and 1. If θ q = 0, Γ k q = 0 and q i = q i , there is no demand uncertainty. If θ q = 1, Γ k q = |N k |, each customer demand q i can take any value in the interval [q i −q i , q i +q i ]. Broadly speaking, Γ k q imposes an upper bound on the number of customers with high demand uncertainty and is determined by θ q and |N k |. Equation (12) shows that the overall travel uncertainty set U t is the Cartesian product of the travel time uncertainty set U k t for each vehicle. The parameters in equation (13) have similar meanings to those in equation (11). A k denotes the set of arcs on the route of vehicle k. t ij denotes the nominal value of uncertain travel time t ij , andt ij represents the maximum deviation from the nominal value for each arc (i, j) ∈ A k . β ij is the auxiliary variable. θ t is the travel time uncertainty budget coefficient, and has a value of between 0 and 1. The travel time uncertainty budget Γ k t controls the degree of travel time uncertainty and is defined as equaling θ t |A k | . |A k | denotes the number of arcs on the route of vehicle k.
Based on the above defined route-dependent uncertainty sets, we extend the deterministic model outlined in equations (1)-(9) into a robust model by adapting the resource inequalities formulation from Agra et al. (2013). Agra et al. (2013) proposed the resource inequalities formulation for a robust version of the VRPTW with travel time uncertainty based on adjustable robust optimization discussed in Ben-Tal et al. (2009). In adjustable robust optimization, some decision variables are allowed to adapt themselves as uncertain parameters vary in uncertainty sets. As such, we change each decision variable y i in the deterministic Agra et al. (2013), only the demand vector q ∈ ext(U q ) and travel time vector t ∈ ext(U t ) need to be considered in the robust formulation. ext(U q ) and ext(U t ) are sets that contain all the extreme points of sets U q and U t , respectively. Thus, the considered robust version of the VRPTW given the route-dependent uncertainty sets can be formulated as follows: s.t.
(2), (3), (4), (5) Constraints (7) and (8) are rewritten for every q ∈ ext(U q ) and t ∈ ext(U t ) as shown in constraints (15) and (16). One key observation from the above formulation is that the number of decision variables y i (q, t) and the number of constraints (15) and (16) increase with the size of the sets N , ext(U q ), and ext(U t ). Notably, finding all the extreme points in the sets U q and U t is a non-trivial task. However, the column-and-row generation algorithm discussed in Agra et al. (2013) can be used to solve the above formulation with small-sized instances.

A two-stage algorithm for the robust VRPTW
The above adapted robust formulation always finds it difficult to handle large-sized VRPTW instances using standard solvers. Therefore, to generate high-quality solutions with largesized instances, we devise a simple two-stage algorithm based on a modified AVNS heuristic.
Subsection 4.1 provides an overview of the two-stage algorithm, followed by a detailed description of the modified AVNS heuristic in Subsection 4.2.

Overview of the two-stage algorithm
Given there are two lexicographical objectives in the considered robust version of the VRPTW, the proposed two-stage algorithm separates the optimization process into two stages. The first stage minimizes the total number of routes and, therefore, the number of vehicles used.
The second stage minimizes the total travel distance in the overall routing strategy. Both stages rely on a modified AVNS heuristic for optimization. An overview of the two-stage algorithm is shown in Figure 1.

A modified AVNS heuristic
The AVNS heuristic was first proposed by Stenger et al. (2013). It incorporates an adaptive mechanism that biases the random shaking step in the basic VNS algorithm (Mladenović and Hansen, 1997). The AVNS heuristic has two major advantages. First, several selected routes, and the sequence of nodes in those routes, can be modified during the adaptive shaking step in one iteration, which produces strong diversification possibilities to escape the local minimum. Second, the algorithm can adapt based on recent search performance, so running times with large-sized instances are reasonable. As discussed in Subsection 4.1, we use a modified AVNS heuristic in both stages of the two-stage algorithm. The modified AVNS heuristic considers the characteristics of the studied robust version of the VRPTW.
The structure of the modified AVNS heuristic is shown in Algorithm 1 in pseudo-code.
Algorithm 1 The structure of the modified AVNS heuristic 1: Define adaptive shaking neighborhood structures N shake I with I = 1, · · · , I max and the local search neighborhood structures N search J with J = 1, · · · , J max 2: Generate an initial solution R initial and improve R initial by local search 3: Set current solution R current ← R initial and I ← 1 4: repeat Update the value of parameters in the adaptive selection mechanism 21: until a given time limit T or a given number Iter of iterations without improvement The main components of the modified AVNS heuristic consist of an adaptive shaking step and a local search step that are repeated until a given time limit T is met or a given number of main iterations Iter without improvement has been reached. Initially, two sets of neighborhood structures {N shake I |I = 1, · · · , I max } and {N search J |J = 1, · · · , J max } are defined. The initial solution R initial is generated, then improved through a local search. In the adaptive shaking step, a shaking solution R shake is generated in the I-th neighborhood N shake I given the current solution R current . In the local search step, a local optimum solution R search is determined by a local search method with neighborhoods N search Specifically, a randomized variable neighborhood descent (RVND) heuristic is used in the local search step. If R search is accepted, it replaces the current solution R current and I is reset to one. Otherwise, R search is discarded and I is increased by one. Since infeasible solutions are allowed in the search process, the next subsection explains the feasibility check and how the value of the objective function is calculated for a given solution, followed by a detailed description of each of the major steps of the modified AVNS heuristic in subsequent subsections.

Feasibility check and objective calculation
Infeasible solutions are allowed in the search process because, in tightly constrained problems, the local search process often quickly converges to local optimal solutions (e.g., the considered robust version of the VRPTW), but penalties are added to the objective function for violated constraints. However, allowing for infeasible solutions creates the challenge of checking route feasibility and calculating the objective function value of a given solution. As mentioned in Subsection 3.2, uncertain travel times and customer demands take their values from the route-dependent uncertainty sets. Thus, the largest possible vehicle load and the latest vehicle arrival time at each customer node on a route can be used to check a route's feasibility.
The intuition behind this notion is that if the largest possible vehicle load on a given route is greater than the vehicle capacity or the latest possible vehicle arrival time for any customer node falls behind the schedule time window, the route is infeasible.
To calculate the largest possible vehicle load and the latest possible vehicle arrival time, consider a given solution (routing strategy) R that consists of m routes and R = {r 1 , · · · , r m }.
Let r k denote a route served by vehicle k ∈ K, and r k ∈ R. The total number of nodes on route r k is |r k | and the ith node on route r k is denoted as p k i . Thus, route r k can be presented as an ordered sequence of nodes, and r k = {p k 1 , p k 2 , · · · , p k |r k |−1 , p k |r k | }. Note that p k 1 is node 0 and p k |r k | is node n+1, which represent the starting and ending depots, respectively. {p k 2 , · · · , p k |r k |−1 } is included in the customer set N c . As discussed in Subsection 3.2, Γ k q and Γ k t control the number of nodes and segments on route r k that have high levels of demand and travel time uncertainty. Thus, the latest possible vehicle arrival time at the ith node on route r k is denoted as Λ k (i, Γ k t ), and the vehicle's largest possible load is denoted as Θ k (Γ k q ).
Λ k (i, Γ k t ) is calculated using a function originally proposed by Agra et al. (2013), formulated as follows: Equation (18) shows that the latest possible vehicle arrival time at each customer node on route r k is a recursive function of the sequence of the nodes on the route and the parameter Γ k t defined in the uncertainty set U k t . In addition to Λ k (i, Γ k t ), we can also calculate the largest possible vehicle load Θ k (Γ k q ) on route r k by arranging the sequence of nodes in route r k in decreasing order of customer demand and generating a new sequence Note that o k i represents the node with the ith largest customer demand on route r k . Based on this sequence, the largest possible vehicle load Θ k (Γ k q ) on route r k can be calculated by equation (19).
The objective function value of a given solution R can be calculated as follows: where Cost(R) represents the total travel distance in solution R, δ t denotes the penalty factor for a time window violation, and δ d denotes the penalty factor for a vehicle capacity violation in equation (20).

Initial solution
As outlined in Subsection 4.1, the modified AVNS heuristic is used in both stages of the two-stage algorithm. The initial solution for the modified AVNS heuristic is constructed using a sequential best insertion heuristic in the first stage of the algorithm. In the second stage of the algorithm, the final output from the first stage is used as the initial solution for the modified AVNS heuristic.

Adaptive shaking
In the shaking step of the modified AVNS heuristic, several selection methods are employed to initially determine the routes and node sequences of the current solution R current to be involved in the shaking. Then, the shaking solution R shake is generated based on the predefined shaking neighborhood operators. An adaptive mechanism is also incorporated into this step to improve the convergence rate of the heuristic toward good-quality solutions.

Shaking neighborhoods
The choice of shaking neighborhood structures is critical since the shaking step helps the AVNS explore a larger search space. In the modified AVNS, the shaking neighborhoods are generated based on two operators: a cyclic exchange and a sequence reinsert. Note that the sequence reinsert operator can be seen as a special case of the cyclic exchange operator. The cyclic exchange operator, originally introduced by Thompson and Psaraftis (1993), moves nodes among routes in a cyclic way. There are two important parameters for the cyclic exchange operator: the number of routes involved N R and the maximum number of nodes to be exchanged N P . The cyclic exchange operator works as follows. For each route r k , the operator moves the node sequence S k i,N P k starting with node i at a length of N P k to route r k+1 at the former position of sequence S k+1 j,N P k+1 (Ibaraki et al., 2005). An example of a cyclic exchange with three routes is shown in Figure 2(a). Note that, if the total number of current routes is less than the number of routes to be moved, N R is reduced accordingly. Similarly, N P is adjusted if it is greater than the number of customer nodes on a given route. The sequence reinsert operator relocates a node sequence from one route to another. An example of the sequence reinsert operator with two routes is shown in Figure 2(b). Table 1 shows the neighborhood operators used in the shaking process. The modified AVNS considers neighborhood operators that exchange up to four routes simultaneously and move up to 10 nodes among them. Sequence lengths of five or fewer nodes are randomly chosen within the interval [0, min(N P, |r k |)]. Sequence lengths of more than five nodes are assumed to be fixed.

Selection methods
To determine the routes and node sequences to be involved in the shaking, we perform a route selection process and a node sequence selection process. We develop some route and node sequence selection methods and apply several successful methods originating from Stenger et al. (2013) to make the adaptive shaking step more effective. The developed selection methods are based on a potential type of uncertainty associated with a vehicle route, e.g., the largest possible vehicle load or the latest vehicle arrival time at any node.
Route selection. Given the current shaking neighborhood N shake I , the route selection process needs to determine N R routes that will be involved in the shaking. The first of the N R routes is selected using one of the five route selection methods below.
1. Random. Each route has the same probability of being selected.
2. Route length. The probability of selecting a route is proportional to the route length.
The idea is to redistribute nodes on long routes into shorter routes in the hope of reducing the total traveling distance.
3. Largest possible vehicle load on a route. The probability of selecting a route is proportional to the largest possible vehicle load on a route. This is intended to reduce the number of routes at risk of exceeding vehicle capacity because the larger the vehicle load, the more likely it is that the route will violate the vehicle capacity constraint. After selecting the first route, the remaining routes are chosen based on the following rule: the next route is randomly chosen from the three spatially closest routes to the previously selected route. The distance between any two routes is assumed to be the distance between the gravity points of these two routes.
Node sequence selection. Once the routes involved in the shaking have been determined, the node sequence to be exchanged for each selected route is determined using one of the following four methods.
1. Random. Each node sequence is randomly chosen.
2. Distance to the next route. The probability of selecting a node sequence is inversely proportional to the distance between the gravity points of the node sequence and the route they will be inserted into.
3. Customer slack time. The probability of selecting a node sequence is inversely proportional to the corresponding total customer slack time of the sequence. This is intended to remove time-sensitive nodes from the route, which may increase the route's robustness.

4.
Missed time windows. The probability of selecting a node sequence is proportional to the total number of missed time windows in that sequence. The intention is to remove nodes with violated time windows and reinsert them into other routes, which may result in an overall feasible solution. If the route is feasible, the node sequence is chosen at random.

Adaptive mechanism
At each shaking step, we use the roulette wheel selection procedure described in Pisinger and Ropke (2007) to determine the route and node selection methods. In the beginning, each route or node sequence selection method is assigned the same probability. Then, the probability of each method is dynamically updated after γ main iterations based on a scoring system. The scoring system has the following rules. If a new overall best solution is achieved after applying a selection method, then a score of nine is added to the method. If the current solution is improved, a score of three is added. If the solution is worse than the current one, but accepted according to the acceptance criterion, a score of one is added to the method.
Suppose the weights of s selection methods are w i (i = 1, · · · , s) and the current score of method i is π i . If the number of times the method was selected since the last weights update is φ i , the new weight w i can be computed by w i (1 − ρ) + ρ(π i /φ i ) every γ main AVNS iterations, where ρ ∈ [0, 1]. Note that π i and φ i are reset to zero once the corresponding parameters are updated.

Local search
The solution generated within the shaking step R shake is then improved by a RVND heuristic. The RVND heuristic employs a set of neighborhood operators in random order and implements the best improvement strategy to prevent poor-quality local optimal solutions.
Six neighborhood operators are used in the RVND: intra-route swap, intra-route reinsert, intra-route 2-opt, inter-route swap, inter-route reinsert, and inter-route 2-opt. A graphic illustration of the six neighborhood operators appears in Figure 3.  Figure 3, the intra-route swap operator exchanges a node i with another node j located in the same route, while the inter-route swap operator exchanges a node i with another node j located in two different routes. The intra-route reinsert operator removes one node i from one route and inserts it into the middle of node j − 1 and j of the same route, while the inter-route reinsert operator removes one node i from one route and inserts it into the middle of node j − 1 and j of another route. Note that the inter-route reinsert operator is able to remove the only node in a route. The resulting empty route is then deleted because it is no longer necessary. The intra-route 2-opt operator disconnects (i − 1, i) and (j, j + 1) in one route, then reconnects (i − 1, j) and (i, j + 1). The inter-route 2-opt operator removes the connection between (i − 1, i) and (j − 1, j) where i, j are in two different routes; and (i − 1, j) and (j − 1, i) are reconnected, which exchanges the second sections of the two affected routes.

Acceptance decision
The solution generated from the local search R search is compared to the current solution R current . If R search is accepted, it becomes the current solution R current and I is reset to one.
The acceptance criterion used in the heuristic was proposed by Hemmelmayr et al. (2009), and is inspired by a simulated annealing mechanism. Specifically, if f (R search ) < f (R current ), . η is the temperature parameter which is decreased from its initial value η 0 by factor η − after every AVNS iteration. We also reset η to η 0 after non-improving iterations to increase diversification of the solutions.

Computational experiments
In this section, we describe the design of the test instances for the considered robust version of the VRPTW, followed by a discussion on the parameter settings for the proposed twostage algorithm in Subsection 5.1. Subsection 5.2 presents the results of each experiment along with a detailed analysis.

Experiment description and parameter settings
The test instances were derived from the well-known Solomon's instances (Solomon, 1987).
These instances were designed to test the VRPTW with deterministic customer demands and travel times. The instances are grouped into six datasets, called R1, R2, C1, C2, RC1, and RC2. We refer to the names of each data set and the instances they contain interchangeably. The R1 and R2 instances include random customer locations, the C1 and C2 instances contain clustered customer locations, and the RC1 and RC2 instances contain a mix of random and clustered customer locations. However, the time windows in the R1, C1 and RC1 instances tend to be narrow whereas, in the R2, C2, and RC2 instances, they tend to be wide. We made no changes to the customer locations and time windows from the original Solomon's instances. However, the customer demands and travel times in the considered robust version of the VRPTW are uncertain and they can take any possible value from the route-dependent uncertainty sets. Thus, descriptions of the parameter values for the uncertainty sets U k q and U k t follow. Each nominal value q i (i ∈ N ) and t ij ((i, j) ∈ A) was assumed to equal to the corresponding customer demand and travel time in each Solomon's instance, respectively. In addition, we assumed that the maximal demand deviationq i was 0.2q i , and the maximal travel time deviationt ij was 0.2t ij . As shown in equations (11), the uncertainty budget coefficient θ q determines the uncertainty budget Γ k q , and θ t determines Γ k t as per equation (13). We assumed that θ q = θ t = 0.3 for all R1, C1 and RC1 instances since the time windows in these instances are narrow. We set θ q = θ t = 0.2 for all R2, C2, and RC2 instances because these instances contain wide time windows.
The proposed two-stage algorithm was coded using Matlab 2015a, and the experiments were run on a laptop with a 2.4 GHz dual processor and 8G RAM. We conducted preliminary tests to find the ideal parameter settings for all adapted instances. In the two-stage algorithm, the modified AVNS heuristic was used to do the optimization tasks in both stages. Therefore, we set the stopping criteria T and Iter for the modified AVNS heuristic to 1200 seconds and 200 iterations in the first stage of the algorithm, and to 2400 seconds and 400 iterations in the second stage of the algorithm. All other parameter values for the modified AVNS heuristic were the same in both stages of the algorithm, and were set as follows: the penalty factors were δ t = 1000 and δ d = 500; γ was set to 20 and ρ was set to 0.3 in the adaptive shaking step; and η 0 was set to 20, η − was set to 0.1, and was set to 50 for the solution acceptance decision. All instances were tested ten times, and the best solution for each instance was selected as the result for inclusion in Subsection 5.2.
Evaluating the robustness of the generated solutions was an important issue in these experiments. Therefore, we used a Monte-Carlo simulation to assess the robustness of the final solutions. In these simulation tests, we assumed that each uncertain travel time parametert ij ((i, j) ∈ A) followed a normal distribution with a mean equal to t ij and a standard deviation of 0.2t ij . Each uncertain demand parameterq i (i ∈ N ) was also assumed to have a normal distribution with a mean equal to t ij and a standard deviation of 0.2t ij . We generated 1000 scenarios which have random travel time matrices and demand vectors. The robustness of a solution was determined by calculating how many scenarios: succeeded in serving all customers; missed one customer time window at most; and missed two customer time windows at most among the 1000 generated. The number of missed customer time windows for a given solution can be counted by examining customer time windows and vehicle capacities using the current scenario's travel time matrix and demand vector. It is worth noting that, although we evaluated the robustness of the solutions based on the assumption of uncertain parameters with normal distributions, it is not necessary to know the exact distributions to produce a robust solution with the proposed robust formulation and the two-stage algorithm.
The developed methods can tackle the VRPTW with uncertain parameters of any symmetric distribution as long as the parameter values are drawn from the proposed route-dependent uncertainty sets.

Computational Results
The following subsections present the results of a series of experiments we conducted to evaluate how a range of factors affect the solutions of the considered robust version of the VRPTW. By analyzing the effect of these factors, it is possible to investigate whether data variations in a factor or a type of uncertainty are likely to generate riskier solutions. The factors examined include demand and travel time uncertainty, customer location patterns, time window length and vehicle capacity, and some key parameters of the route-dependent uncertainty sets. A performance analysis of the proposed two-stage algorithm solving several adapted small-sized instances has also been included.

The effect of demand and travel time uncertainty
To see the detailed effects of demand and travel time uncertainty on vehicle routes, we compared the solutions of the considered robust version of the VRPTW with the solutions of the deterministic VRPTW using the R1 instances. The results of this comparison are shown in Table 2. For each test instance in Table 2, "Det" row shows the details of the deterministic solution and "Rob" row shows the details of the robust solution. The columns "N.V.", "OBJ", and "Increase" respectively represent the number of vehicles used, the total travel distance, and the percentage increase in the total travel distance of the robust solution vs. the deterministic solution in each instance. Each solution's robustness, as measured by the Monte-Carlo simulation test, is reported under the columns "Total", "Demand", and "Travel time". The columns under "Total" show the results of the simulation tests that consider both demand and travel time uncertainty. The columns under "Demand" show the results of the simulation tests that only consider demand uncertainty, while the columns under "Travel time" show the results of the simulation tests that only consider travel time uncertainty. Columns "V 0 ", "V 1 ", and "V 2 " respectively mean the probability of serving every customer, the probability of not serving one customer at most, and the probability of not serving two customers at most within a solution. Note that the deterministic solutions for all large-sized test instances were obtained by the proposed two-stage algorithm. However, we made some small modifications to the feasibility check and objective calculation procedure discussed in Subsection 4.2.1 because uncertainty is not considered in the deterministic VRPTW.  In analyzing the detailed results for the R1 instances, one clear observation is that the probability of servicing every customer with the deterministic solutions is very low. This shows that deterministic solutions are very fragile when travel times and customer demands are uncertain. The R1 instances contain randomly distributed customers with narrow time windows, short scheduling horizons, and small-capacity vehicles, which is likely the reason for this fragility. However, we observe that, in most instances, the robustness of the routes improved significantly with the robust solutions. For example, in instance R112, the probability of servicing all customers within the prescribed time window increased to 80.6% compared to 0.2% with the deterministic solution. In instance R101, this same probability only increased to 30.3% with the robust solution, but the probability of not serving one customer at most increased from 1.4% to 64.9%. We consider this to be a satisfactory result given the most of the customers with narrow time windows in the R101 instance. One interesting finding is that travel time uncertainty seems to be the main cause of route infeasibility for most of the R1 instances, whereas customer demand uncertainty does not appear to have a strong impact. For example, in the deterministic solutions for the R101, R102, R105, and R106 instances, the probability of serving every customer was 100% when only demand uncertainty was considered. This may be because some of the R1 instances contain many narrow time windows with short scheduling horizons. Hence, fewer customers are serviced on each route and demand uncertainty can be offset by vehicle capacity, yet travel time uncertainty can still easily result in violations to these narrow customer time windows.

The effect of customer location patterns
As      Table 5. Results for the C1 instances So, although the robust solution greatly reduces the risk of a missed customer time window, decision-makers may still need to assess the trade-off between total travel distance and route robustness in situations with clustered location patterns.

The effect of time window length and vehicle capacity
Contrary to the instances in the R1, C1 and RC1 datasets, the R2, RC2, and C2 instances tend to contain wide customer time windows with a long scheduling horizon and largecapacity vehicles. These characteristics allow many customers to be served along the same route. To show the effect of wide customer time windows and large-capacity vehicles on the planned routes, we present the details of the deterministic and robust solutions for the R2, RC2 and C2 instances in Tables 6, 7, and 8, respectively.   Table 7. Results for the RC2 instances   instances are more robust than the deterministic solutions for the R1, RC1 and C1 instances.
For example, the average probability of servicing every customer with the deterministic solutions for the C2 instances was 54.2%, which indicates that these deterministic solutions are not very fragile. In the robust solutions for the R2, RC2, and C2 instances, the total travel distance only increased a little, and the number of vehicles used stayed the same in most of the instances. However, the route robustness with these solutions improved significantly. For example, the average probability of serving every customers was 95.3%, 94.7%, and 99.5% for the R2, RC2 and C2 instances, respectively. These results indicate that the solutions are very robust given both travel time and demand uncertainty. We can therefore conclude that a high level of route robustness can be reached for instances with wide time windows and large-capacity vehicles using the robust solutions at almost no additional cost (i.e., little increase in the total travel distance and the number of vehicles used). However, failing to serve all customers is more prevalent in deterministic solutions given instances with narrow time windows, and ensuring a high level of route robustness is much more expensive, as shown in Tables 3, 4, and 5.
From Tables 6 and 7, one interesting finding is that the robust solutions use more vehicles but travel less total distance than the deterministic solutions, and yet still improve the route robustness significantly in some of the R2 and RC2 instances. For example, in instance R202, the robust solution used one more vehicle than the deterministic solution, but the total travel distance was 6.69% less and the probability of serving all customers increased to 88.4%. Instances R207, R211, and RC202 show similar results. Thus, robust solutions may be very attractive in some instances because they can decrease the total travel distance and increase the route robustness with the addition of only a small number of vehicles. It is also worth noting that adding more vehicles does not necessarily decrease the total travel distance, but this can be effective strategy for improving the robustness of the planned routes.

The effect of uncertainty set parameters
A solution's robustness is affected by the parameters of the route-dependent uncertainty sets outlined in equations (10)-(13). Hence, our next set of experiments were designed to analyze the effects of the uncertainty budget coefficients θ q and θ t and the maximal deviation of uncertain demandq i (i ∈ N c ) and travel timet ij ((i, j) ∈ A). Instances RC101 and RC201 were selected for analysis, and the assumptions made during these experiments follow. The ratio ν q between the maximal demand deviationq i and its nominal value q i stayed the same for all customer nodes i ∈ N c ; the ratio ν t between the maximal travel time deviationt ij and its nominal value t ij stayed the same for all arcs (i, j) ∈ A. These ratios are referred to as the uncertainty ranges. Note that we analyzed the uncertainty ranges and the uncertainty budget coefficients separately. When analyzing the effect of the uncertainty budgets, the uncertainty ranges were assumed to be fixed, and were set to ν q = ν t = 0.2. When analyzing the effect of the uncertainty ranges, the uncertainty budget coefficients were assumed to be fixed, and were set to θ q = θ t = 0.3 for the RC101 instance and to θ q = θ t = 0.2 for the RC201 instance. The robust solutions for the RC101 and RC201 instances given different uncertainty ranges are shown in Table 9. The robust solutions given different uncertainty budget coefficients are shown in Table 10. Let us first focus on the effect of different uncertainty ranges with the robust solution for the RC101 instance, shown in Table 9. As the value of the uncertainty ranges increased from 0 to 0.5, the number of vehicles used increased from 14 to 20 and the total travel distance increased from 1696.95 to 2139.63 with the robust solution. Instance RC101 contains many narrow time windows, small-capacity vehicles, and a short scheduling horizon. Thus, more vehicles need to be scheduled to meet the prescribed time windows as the value of the uncertainty ranges increases and this results in a greater total travel distance. However, the effect of varying uncertainty ranges with the RC201 instance is quite different. As the value of the uncertainty ranges increased from 0 to 0.4, the number of vehicles used stayed the same but the total travel distance increased from 1413.52 to 1603.16 with the robust solution.
Instance RC201 contains wide customer time windows, large-capacity vehicles, and a long scheduling horizon. Hence, the sequence in which each customer is visited along a route can simply be reordered to achieve a high level of route robustness without increasing the number of vehicles used when the value of the uncertainty ranges is relatively small. However, as the value of the uncertainty ranges increased from 0.4 to 0.5, the number of vehicles used increased by one and the total travel distance decreased from 1603.16 to 1422.77. The reason is that one more vehicle is needed to deal with the greater levels of uncertainty and the vehicle routes are rearranged which significantly decreases the total travel distance.   Table 10, the number of vehicles used stayed the same and the total travel distance increased from 1413.52 to 1477.13 as the value of the uncertainty budget coefficients increased from 0 to 0.4, after which the robust solution stayed the same in all respects. Although, notably, the resulting solution with a 90% probability that every customer will be served had a high level of route robustness.
Based on the above analysis, we can conclude that adjusting the uncertainty ranges and/or the uncertainty budget coefficients can affect the solutions generated for the considered robust version of the VRPTW. The results also highlight two insights for decisionmakers: (1) decision-makers should select routing strategies with a high level of robustness for practical applications characterized by wide customer time windows and large-capacity vehicles; and (2) decision-makers may have to balance the total travel distance and the number of vehicles used (i.e., the cost) with the robustness of the planned routes when selecting an ideal routing strategy for practical applications characterized by narrow customer time windows and small-capacity vehicles.

Algorithm performance
To test the performance of our proposed two-stage algorithm, we adapted 24 of Solomon's instances by selecting the first ten customers in each instance. Given the vehicle capacity of each Solomon's instance is too large to be meaningful, we reset the vehicle capacity to 75 for the R1 and R2 instances, to 100 for the C1 and C2 instances, and to 150 for the RC1 and RC2 instances. We also reset the uncertainty budget coefficients to θ q = θ t = 0.6. The stopping criteria T and Iter for the modified AVNS heuristic were set to 50 seconds and 20 iterations in both stages of the two-stage algorithm due to the small size of the instances.
We retained the values described in Subsection 5.1 for all other parameters of the modified AVNS heuristic. We tested the proposed two-stage algorithm ten times with each instance and compared the resulting solutions to the optimal solutions generated by solving the robust formulation outlined in Subsection 3.2. The comparison results are shown in Table 11. Note that the column-and-row generation algorithm discussed in Agra et al. (2013) was adapted to solve the robust formulation. The algorithm was coded in GAMS software with the CPLEX 12.3 solver. To further assess the quality of the proposed algorithm, we solved a robust version of the VRP with deadlines under demand and travel time uncertainty discussed by Lee et al. (2012). Lee et al. (2012) developed an exact algorithm to solve the problem with the adapted small-sized Solomon's R1 instances. However, the robust version of the VRP with deadlines considers only one objective and different uncertainty sets. Thus, we adapted the modified AVNS heuristic described in Subsection 4.2 to solve the problem. Each instance was also tested ten times. The comparison between the exact algorithm in Lee et al. (2012) and the modified AVNS heuristic is shown in Table 12. In Table 11, the columns under "AC&RG algorithm" show the details of the solutions obtained by adapting the column-and-row generation algorithm discussed in Agra et al. (2013). The columns "N.V.", "Optimal", and "Time" represent the number of vehicles used, the optimal objective value, and the computation time of the algorithm for each instance, respectively. Similarly, the columns under "Two-stage algorithm" show the details of the solutions obtained using the proposed two-stage algorithm. The columns "Best", "Avg", and "Time Avg" represent the best objective value, the average objective value, and the average computational time of ten runs, respectively. Table 11 shows that the proposed two-stage algorithm was able to find the optimal solution for the considered robust VRPTW with each small-sized instance. The two-stage algorithm was very stable in terms of computation time, taking less time to generate solutions than the adapted column-and-row generation algorithm for more than half the test instances. However, the adapted column-and-row generation algorithm was faster when solving several instances including R101, R111, R201, R210, C101, C108, C201, and RC201. The reason is related to the instance structures.
For example, the R101 and C101 instances contain many narrow time windows and small-capacity vehicles. This type of instance structure is particularly suited to the CPLEX solver, and hence it is able to find an optimal solution in every iteration of the algorithm very quickly.
For instances with many wide time windows and large-capacity vehicles such as instances R201, C201 and RC201, the adapted column-and-row generation algorithm only requires a few iterations to generate the optimal robust solution. In Table 12, the columns under " Lee et al. (2012)" show the details of the solutions obtained by the exact algorithm proposed in Lee et al. (2012). The columns "Optimal" and "Time" represent the optimal objective value and the computation time of the exact algorithm for each instance, respectively. The columns under "Modified AVNS" summarize the details of the solutions obtained by the modified AVNS heuristic. The columns "Best", "Avg", and "Time Avg" represent the best objective value, the average objective value, and the average computational time of ten runs, respectively. Table 12 clearly shows that the best solution produced by the heuristic is the same as the optimal solution generated by the exact algorithm for each instance. Moreover, with the exception of the first instance, the modified AVNS heuristic was much faster than the exact algorithm.

Conclusions
In this paper, we considered a robust version of the VRPTW with demand and travel time uncertainty. To address the problem, we assumed the uncertain parameters took values from the route-dependent uncertainty sets and adapted a robust mathematical formulation. Due to the complexity of the problem, the adapted robust model was only able to tackle smallsized instances using standard solvers. To solve large-sized instances, we proposed a two-stage algorithm based on a modified AVNS heuristic. Extensive computational experiments were conducted using modified versions of Solomon's benchmark instances. The results show that the two-stage algorithm was able to produce optimal solutions for small-sized instances and the planned routes generated for large-sized instances were significantly more robust at the cost of a small increase in the total travel distance and the number of vehicles used.
We performed a comprehensive analysis of the results and several managerial insights are revealed as follows: (1) Generating vehicle routes without considering demand and travel time uncertainty yields very fragile routing strategies that often result in missed time windows.
(2) Incorporating more vehicles into a routing schedule and reordering the sequence of customer visits within each vehicle route can sometimes increase the robustness of a given routing strategy.
(3) Highly robust routing strategies can be generated for settings characterized by wide customer time windows and large-capacity vehicles at little additional cost.
(4) Failing to serve a customer within the prescribed time window in situations characterized by narrow time windows and small-capacity vehicles is more prevalent when the routing strategy is generated with a deterministic model. And reaching a high level of robustness is much more expensive -more vehicles and much longer travel distances are required.