# A Scenario Decomposition Algorithm for Strategic Time Window Assignment
Vehicle Routing Problems

Anirudh Subramanyam  Department of Chemical Engineering, Carnegie Mellon
University, United States Akang Wang  Department of Chemical Engineering,
Carnegie Mellon University, United States Chrysanthos E. Gounaris  Department
of Chemical Engineering, Carnegie Mellon University, United States

###### Abstract

We study the strategic decision-making problem of assigning time windows to
customers in the context of vehicle routing applications that are affected by
operational uncertainty. This problem, known as the Time Window Assignment
Vehicle Routing Problem, can be viewed as a two-stage stochastic optimization
problem, where time window assignments constitute first-stage decisions,
vehicle routes adhering to the assigned time windows constitute second-stage
decisions, and the objective is to minimize the expected routing costs. To
that end, we develop in this paper a new scenario decomposition algorithm to
solve the sampled deterministic equivalent of this stochastic model. From a
modeling viewpoint, our approach can accommodate both continuous and discrete
sets of feasible time window assignments as well as general scenario-based
models of uncertainty for several routing-specific parameters, including
customer demands and travel times, among others. From an algorithmic
viewpoint, our approach can be easily parallelized, can utilize any available
vehicle routing solver as a black box, and can be readily modified as a
heuristic for large-scale instances. We perform a comprehensive computational
study to demonstrate that our algorithm strongly outperforms all existing
solution methods, as well as to quantify the trade-off between computational
tractability and expected cost savings when considering a larger number of
future scenarios during strategic time window assignment.

Keywords: vehicle routing under uncertainty, time window assignment, service
consistency, stochastic programming, decomposition, branch-and-bound.

##  1 Introduction

The commitment to deliver (or pickup) goods within scheduled time windows is a
common practice in several real world distribution networks. In many
industries, these time windows are mutually agreed upon by the distributor and
customer through long-term delivery contracts. For example, in a distribution
network of retailers, it is common that deliveries to a retail store are
always made on the same day of the week (at about the same time) for an entire
year [33, 39]. Likewise, in maritime distribution of liquefied natural gas, a
central planning activity is to design and negotiate contractual agreements of
annual delivery plans that specify delivery dates and corresponding delivery
quantities to customers [40]. From the customer’s point of view, this is
crucial for efficient inventory management and scheduling of personnel to
process the delivery. From the distributor’s point of view, it reduces the
variability across repetitive deliveries and exposes efficiencies that add up
to significant cost savings. Short- and medium-term contracts of similar
nature can be also found in small-package shipping where, for instance,
courier companies provide a delivery time window to customers receiving
sensitive packages [20]. Other examples of applications where such operations
are typical include, among others, attended home delivery in e-commerce
businesses [1] and internet installation services [37].

Once a time window has been agreed upon and communicated to the customer, the
distributor must attempt to meet it on an operational (e.g., daily) basis as
well as possible. This is done by solving a Vehicle Routing Problem with Time
Windows (VRPTW) to determine a delivery schedule that adheres to the agreed
time windows. The assigned time windows strongly influence the structure of
feasible delivery schedules and, hence, the daily incurred distribution costs.
Therefore, a natural choice is to assign time windows based on the arrival
times at customer locations in the optimal (i.e., minimum cost) vehicle
routing schedule. However, this seemingly optimal decision may become highly
suboptimal in the presence of operational uncertainty.

In reality, operational level information (such as customer demands or travel
times) is often not known with certainty at the strategic level when time
windows are to be decided. For example, the demand volume of a customer
typically fluctuates per delivery. Similarly, travel times vary on a day-to-
day basis (e.g., because of unpredictable traffic conditions). The true values
of these operational parameters are not known far in advance, and often may
become known only on the day of delivery before the vehicles are dispatched.
This makes the strategic assignment of time windows a non-trivial task.
Indeed, if one utilizes only nominal values of the uncertainty when assigning
time windows, then it will often lead to situations in which the distribution
costs are unacceptably high, since the nominal delivery schedule may no longer
be feasible, let alone optimal, in such cases. Fortunately, with the
increasing availability of data, distributors can readily obtain forecasts of
uncertain operational parameters (e.g., as perturbations from their nominal
values). It is possible to take advantage of this information and assign time
windows in a way that will lead to low distribution costs in the long run. The
goal of this paper is to study the problem of strategic time window assignment
in the presence of operational uncertainty.

Our paper builds upon the work of [33], which introduced the Time Window
Assignment Vehicle Routing Problem (TWAVRP). The TWAVRP consists of assigning
time windows of pre-specified width within some exogenous time windows to a
set of known customers. The exogenous time windows typically correspond to
operating hours of the customer but may also arise from hours-of-work or other
government regulations. The work of [33] studies the TWAVRP under situations
in which the demand volume of the customers is unknown and subject to
uncertainty. However, a finite set of “scenarios,” each describing a possible
realization of demand for every customer, is assumed to be given with known
probability of occurrence. This information is used to formulate a two-stage
stochastic program, in which the first-stage decisions are to assign time
windows, while the second-stage decisions are to design vehicle routing
schedules satisfying the assigned time windows, one for each of the demand
scenarios. The objective is to minimize the total routing costs, averaged over
the postulated scenarios. A similar modeling approach is followed in [32],
with the only difference that the first-stage time windows are selected from a
finite set of a priori constructed windows; this problem is referred to as the
discrete TWAVRP to distinguish it from the original continuous TWAVRP. In this
paper, we consider both cases, and we shall in fact allow also for the
generalized case in which feasible time window assignments lie in a continuous
set for some portion of the customer base and in a discrete set for the
remaining portion.

Algorithms to solve the aforementioned stochastic programming models have been
proposed in [33, 10] for the continuous version, and in [32] for the discrete
version of the problem. The algorithms of [33, 32] are based on branch-price-
and-cut and can solve instances with 25 customers and 3 demand scenarios to
optimality, while the algorithm of [10] is based on branch-and-cut and can
address instances containing 40 customers and 3 scenarios. Several heuristics
have also been proposed in [32] for the discrete setting that can address
instances containing up to 60 customers. Recently, [31] studied a variant of
the TWAVRP with time-dependent travel times and proposed a branch-price-and-
cut algorithm that can solve instances with 25 customers and 3 demand
scenarios.

A problem that is closely related to the strategic TWAVRP is the Consistent
Vehicle Routing Problem (ConVRP) [18], which is motivated in the context of
operational level planning. The ConVRP aims to design minimum cost vehicle
routes over a finite, multi-day horizon to serve a set of customers with known
demands. The goal is to design routes that are _consistent_ over time; this
translates to satisfying any of the following requirements each time service
is provided to a customer: (i) arrival-time consistency, wherein the customer
should be visited at roughly the same time during the day, (ii) person-
oriented consistency, in which the customer should be visited by the same
driver, and whenever applicable, (iii) load consistency, for which a customer
should receive roughly the same quantity of goods. We refer the reader to [22]
for an overview of this problem and its applications.

Conceptually, the assigned time windows in the TWAVRP (which are also referred
to as the endogenous time windows) serve to satisfy the arrival-time
consistency requirement of the ConVRP, which requires that every customer be
visited at roughly the same time whenever service is requested. Formally, the
ConVRP requires that the difference between the earliest and the latest
arrival times at each customer location must differ by no more than some pre-
specified constant bound, which is referred to as the maximum allowable
arrival-time differential. This bound is analogous to the pre-specified width
of the endogenous time window in the continuous TWAVRP, and this equivalence
between the two problems has been previously acknowledged in [31, 33].

The equivalence between the TWAVRP and the arrival-time ConVRP has two
important consequences. First, we observe that, in the most general case, the
ConVRP allows for the possibility that not all customers require service in
all time periods and that operational parameters (such as customer demands or
travel times) differ from one time period to the other. Translated in the
context of the TWAVRP, this allows us to address applications in which a
fraction of the customer base does not require frequent (e.g., daily) service
(by considering scenarios where certain subsets of customers have no demand),
as well as to treat uncertainty in a wider variety of parameters such as
travel and service times (by considering scenarios in which their values
represent perturbations from some nominal value). We therefore study the
TWAVRP under a more general definition than what has been previously
considered in the literature, in which it is possible to incorporate
uncertainty in several operational parameters at once. However, we do remark
that, as is the case with traditional stochastic programming models, the
simultaneous treatment of uncertainty in several parameters may come at the
cost of an explosion in the number of scenarios that have to be considered.

The second consequence of the equivalence between the TWAVRP and the arrival-
time ConVRP is that any algorithm developed for the latter can be used to
obtain solutions for the former. In this paper, we adapt a decomposition
algorithm proposed in [35] for the Consistent Traveling Salesman Problem
(ConTSP), the single-vehicle variant of the ConVRP that focuses purely on the
aspect of arrival-time consistency, to obtain a new algorithm for the TWAVRP.
Our method can be viewed as a scenario decomposition algorithm in the language
of stochastic programming, and is not based on branch(-price)-and-cut that has
been the de facto approach for solving TWAVRP models. Our algorithm has the
attractive features of modularity and scalability. It is modular in
accommodating (i) continuous and discrete time windows, (ii) any VRPTW solver
(exact or heuristic), (iii) routing-specific constraints (e.g., heterogeneous
fleets), and (iv) generic scenario descriptions. Moreover, it can be readily
parallelized which allows postulating a large number of scenarios of the
uncertainty. Together with (ii) above, this means that our algorithm is also
scalable. The distinct contributions of our work may be summarized as follows.

  * •

We generalize the definition of the TWAVRP. In particular, we study problems
with scenario-based models of uncertainty in which any operational parameter
may be uncertain and in which the endogenous time windows may be chosen from
either continuous or discrete sets.

  * •

We propose an exact, scenario decomposition algorithm for the TWAVRP. The
algorithm outperforms all state of the art methods, solving 54 out of 81
previously open instances. Furthermore, it can be easily parallelized, utilize
any available VRPTW solver in a “black box” fashion, and be readily modified
as a heuristic to solve large-scale instances.

  * •

We propose a new class of path-based disjunctions for the TWAVRP. We show via
numerical experiments that these disjunctions, which manifest as branching
rules in our algorithm, significantly improve the performance of the latter.
In view of the relationship between the TWAVRP and ConVRP, these disjunctions
are new for the ConVRP as well.

  * •

We conduct experiments with a parallel implementation of our algorithm to
solve instances consisting of up to fifteen scenarios, representing a five-
fold increase compared to existing literature. We use these solutions to
elucidate, via out-of-sample simulations, the cost savings that are to be
expected when considering more scenarios during time window assignment.

The rest of this paper is organized as follows. Section 2 reviews the relevant
literature; in Section 3, we provide a general mathematical definition of the
TWAVRP; in Section 4, we describe our solution algorithm for this problem;
Section 5 elaborates on important implementation details of our algorithm;
Section 6 presents computational results on existing as well as new datasets;
and, finally, we conclude in Section 7.

##  2 Related Literature

This section reviews papers in the vehicle routing literature, other than
those mentioned in the introduction, which deal with aspects of (i)
endogenously imposed time windows, (ii) consistent service considerations, and
(iii) stochastic or uncertain parameters. We choose not to review the
extensive literature on the VRPTW; instead, we refer interested readers to
[11].

The authors of [20] study a time window assignment problem that is encountered
by courier companies who must quote delivery time windows to customers
receiving sensitive packages. In this problem, which they refer to as the VRP
with Self-Imposed Time Windows, travel times are uncertain but all customers
and their demands are known a priori. Using this information, the service
provider must simultaneously determine (i) a single routing plan to serve all
customers, and (ii) time window assignments that will be quoted to the
customers before the vehicles depart from the depot. The objective is to
minimize the sum of (deterministic) routing costs and (expected) overtime and
tardiness penalty costs. Since travel times are uncertain, the key challenge
is to determine the optimal placement of time windows (along each vehicle
route) in step (ii) so as to avoid penalties due to delays. The uncertainty in
travel time along each arc is modeled via a discrete set of “disruption”
scenarios (each representing a deviation from some nominal value). Under the
assumption that at most one arc will be disrupted on any vehicle route, the
authors propose a Tabu Search heuristic for route generation and a linear
programming approach for time window placement that inserts “time buffers”
along each vehicle route. Recently, with a goal to solve the same problem, the
authors of [38] extended the work of [20]. On the one hand, they relax the
assumption that only one arc will be disrupted and use probabilistic chance
constraints to guarantee reliable service. On the other hand, they propose an
alternative model of uncertainty in which the stochastic deviations in travel
times are modeled as continuous gamma-distributed random variables. Finally,
by considering also the width of the time window (along with its placement) as
a decision variable, they propose an Adaptive Large Neighborhood Search
solution procedure.

A related line of work is the so-called Time Slot Management Problem that is
motivated in the context of attended home delivery in e-commerce businesses
[1]. Here, customers place online orders for products (e.g., groceries) and,
during this process, they select one time window (amongst a number of
available ones) in which they want their product to be delivered. From the
service provider’s point of view, the challenge is to design a finite set of
time windows (instead of just one) to offer to potential customers in
different zip code areas. The problem is complicated by the fact that, during
the design phase, the set of customers as well as their demand is not known
with certainty. The objective is to design time windows that would not only
yield low distribution costs in the long run, but also satisfy marketing or
regulatory requirements. In general, existing approaches (e.g., see [1, 19,
8]) deal with uncertainty by simply using expected values of customer demand
whose temporal distribution is either assumed to be uniform over the offered
time windows or determined via simulations. The expected routing cost
associated with a candidate set of time windows is then estimated via coarse
continuous approximation methods (e.g., see [9, 14]) or detailed vehicle
routing models. These estimates are embedded within some heuristic procedure
(e.g., local search) to determine the final set of time windows. We refer the
reader to [2] for an overview of research problems in the area of attended
home delivery, including time slot management. Finally, we mention the work of
[37] who also study a time window assignment problem that is motivated in the
context of home-attended services (e.g., cable installation). Here, as is the
case in attended home delivery, customers dynamically place orders for some
service. However, instead of the customer choosing a time window from a number
of available ones, the service provider must quote a service time window to
the customer at the time of request. Similar to the VRP with Self-Imposed Time
Windows, travel and service times are stochastic and the objective is to
minimize expected delays. However, unlike the latter problem, not all
customers who will be serviced are known at the time when a particular request
is received, and thus the customer base is also stochastic. The authors use
approximate dynamic programming techniques to obtain time window assignments
in real time.

We now review papers that study the ConVRP or its variants. The problem was
introduced in [18] motivated by operations encountered in small package
shipping services offered by courier companies. Since then, there have been
many studies that have focused on trying to design solution algorithms for the
ConVRP (e.g., see [18, 36, 23, 25]). To the best of our knowledge, most of
these approaches are heuristic in nature, and are based on the concept of
generating a “template” routing plan that services only the most frequent
customers; the daily routes are derived from the template by appropriately
modifying it in a way that satisfies all service consistency requirements, in
particular, driver and arrival-time consistency. With a focus on trying to
address purely the requirement of arrival-time consistency, exact solution
approaches for the ConTSP (the single-vehicle variant of the ConVRP) have been
proposed in [34, 35]. Given our result from later in this paper that
establishes equivalence between the TWAVRP and arrival-time ConVRP, one can,
in principle, adapt any of the aforementioned methods to obtain algorithms for
the TWAVRP by ignoring aspects of driver consistency, whenever applicable. We,
however, shall choose the decomposition algorithm of [35] for this purpose,
due to a number of reasons. First, because of its decomposition principle, it
solves the original problem by breaking it down into period-specific routing
problems with time windows. This has the advantage that users can use their
own routing solver as a module to obtain time window assignments and that it
allows easy parallelization to solve instances with many time periods (or
scenarios in the case of the TWAVRP). Second, it can be run in either exact or
heuristic mode by changing the underlying routing solver. Third, as we will
demonstrate, it can be readily extended to solve also the discrete TWAVRP, for
which other ConVRP approaches would not be directly applicable.

We conclude our literature review by mentioning the relationship of the TWAVRP
to Stochastic Vehicle Routing Problems (SVRP). The latter class of problems
also treats parameter uncertainty in the context of vehicle routing. However,
unlike the TWAVRP which is inherently a strategic decision-making problem, the
SVRP is an operational problem. Specifically, in the TWAVRP, the exact values
of all parameters are assumed to be known before the vehicle routes are to be
determined on a particular day. In contrast, in the SVRP, the vehicle routes
must be determined before the parameter values become known, which are only
gradually revealed during the execution of the routing plan. This requires
fundamentally different modeling considerations and corresponding solution
approaches. The most common modeling paradigms are (i) recourse models, (ii)
reoptimization models, and (iii) chance-constrained models. In (i), a planned
solution is designed in the first stage and recourse actions based on a
predetermined policy are taken in the second stage when the uncertainties are
revealed. For example, the capacity of a vehicle may get exceeded en route, if
demands are stochastic at the time of vehicle dispatch; in such cases, a
recourse policy, such as a detour to the depot to empty the vehicle, must be
explicitly incorporated in the model [13, 16]. In (ii), the planned solution
is dynamically modified as the uncertain parameters (e.g., demands or travel
times) become gradually revealed during the execution of the vehicle routes
[29]. Finally, in (iii), probabilistic or chance constraints are used to
explicitly control the level of risk that is acceptable to the decision-maker
[24]. We refer interested readers to [17], who provide an excellent overview
of applications, models and solution algorithms for the SVRP and its variants.

##  3 Notation and Problem Definition

Let $G=(V,A)$ denote a directed graph with nodes $V=\\{0,1,\ldots,n\\}$ and
arcs $A$. Node $0\in V$ represents the unique depot, and each node $i\in
V_{C}\coloneqq V\setminus\\{0\\}$ represents a customer. The operating hours
of the depot are represented by the time window $[e_{0},\ell_{0}]$, where an
unlimited number of vehicles, each of capacity $Q$, are available for service.
Each vehicle incurs a transportation cost $c_{ij}\in\mathbb{R}_{\geq 0}$ and a
travel time $t_{ij}\in\mathbb{R}_{\geq 0}$ if it traverses the arc $(i,j)\in
A$. Furthermore, each customer $i\in V_{C}$ features a demand
$q_{i}\in\mathbb{R}_{\geq 0}$, service time $u_{i}\in\mathbb{R}_{\geq 0}$ and
_exogenous_ time window $[e_{i},\ell_{i}]$ (e.g., representing operating
hours). The key decisions in the TWAVRP are to decide the _endogenous_ time
windows $\tau_{i}\in TW_{i}$ to be assigned to each customer $i\in V_{C}$. The
definition of the feasible time window set $TW_{i}$ may be either of the
following (refer to Figure 1):

  * •

In the _continuous_ setting, the assigned time window must have a pre-
specified width $w_{i}\in\mathbb{R}_{\geq 0}$; that is,
$TW_{i}=\left\\{[y_{i},y_{i}+w_{i}]:e_{i}\leq
y_{i}\leq\ell_{i}-w_{i}\right\\}$, where we assume, without loss of
generality, that $e_{i}\leq\ell_{i}-w_{i}$.

  * •

In the _discrete_ setting, the assigned time windows must belong to a pre-
specified finite set; that is,
$TW_{i}=\left\\{[\underaccent{\bar}{y}_{i1},\bar{y}_{i1}],\ldots,[\underaccent{\bar}{y}_{iN_{i}},\bar{y}_{iN_{i}}]\right\\}$,
where we can assume, without loss of generality, that all $N_{i}$ candidate
time windows are pairwise either non-overlapping or partially
overlapping.111Two completely overlapping time windows $[a,b]$ and $[c,d]$
with $a\leq c\leq d\leq b$ can be replaced with the larger of the two time
windows $[a,b]$. Therefore, the set $TW_{i}$ can be ordered so that
$e_{i}=\underaccent{\bar}{y}_{i1}\leq\ldots\leq\underaccent{\bar}{y}_{iN_{i}}$
and $\bar{y}_{i1}\leq\ldots\leq\bar{y}_{iN_{i}}=\ell_{i}$.222We remark that
$e_{i}=\underaccent{\bar}{y}_{i1}$ and $\ell_{i}=\bar{y}_{iN_{i}}$ can be
achieved by preprocessing. If $e_{i}<\underaccent{\bar}{y}_{i1}$, then $e_{i}$
can be shifted forward to match $\underaccent{\bar}{y}_{i1}$, and if
$e_{i}>\underaccent{\bar}{y}_{ib}$, then $\underaccent{\bar}{y}_{ib}$ can be
shifted forward to match $e_{i}$. A similar argument applies for $\ell_{i}$.

We denote by $V_{\text{cont}}\subseteq V_{C}$ and $V_{\text{disc}}\subseteq
V_{C}$ the subset of customers whose feasible time window sets are continuous
and discrete respectively. We note that $V_{\text{cont}}\cap
V_{\text{disc}}=\emptyset$ and $V_{\text{cont}}\cup V_{\text{disc}}=V_{C}$.

time$\vphantom{\ell_{i}}e_{i}$$\vphantom{\ell_{i}}y_{i}$$\vphantom{\ell_{i}}y_{i}+w_{i}$$\ell_{i}$$w_{i}$time$\vphantom{\ell_{i}}e_{i}=\underaccent{\bar}{y}_{i1}$$\vphantom{\ell_{i}}\bar{y}_{i1}$$\vphantom{\ell_{i}}\underaccent{\bar}{y}_{i2}$$\vphantom{\ell_{i}}\bar{y}_{i2}$$\vphantom{\ell_{i}}\underaccent{\bar}{y}_{i3}$$\bar{y}_{i3}=\ell_{i}$
Figure 1: Illustration of continuous (left) and discrete (right) time window
sets $TW_{i}$. In the continuous case, the assigned time window can be any
sub-interval of $[e_{i},\ell_{i}]$ of length $w_{i}$, while in the discrete
case, the assigned window must be one of the intervals
$[\underaccent{\bar}{y}_{i1},\bar{y}_{i1}]$,
$[\underaccent{\bar}{y}_{i2},\bar{y}_{i2}]$ or
$[\underaccent{\bar}{y}_{i3},\bar{y}_{i3}]$.

In practice, operational parameters such as those related to the
transportation network (costs $c$, travel times $t$) or the customers (demands
$q$, service times $u$) are often not known with certainty at the strategic
level when time windows must be allocated. Let $\theta$ denote the set of all
operational parameters, and let $\mathbb{P}$ denote the joint probability
distribution of $\theta$. The goal of the TWAVRP is to assign the time windows
$\tau_{i}\in TW_{i}$ in a way that minimizes the expected cost of routing:

$\begin{array}[]{r@{\;\;}l}\displaystyle\mathop{\text{minimize}}_{\tau}&\displaystyle\mathbb{E}_{\theta\sim\mathbb{P}}\left[\mathtt{VRPTW}(\tau;\theta)\right]\\\
\text{subject to}&\tau_{i}\in TW_{i}\quad\forall i\in V_{C}.\end{array}$ (1)

In the above _stochastic programming_ formulation,
$\mathtt{VRPTW}(\tau;\theta)$ denotes the minimum cost of the vehicle routing
problem with time windows $\tau_{i}$, $i\in V_{C}$ and operational parameters
$\theta$. A formal mathematical definition of $\mathtt{VRPTW}(\tau;\theta)$
follows.

###  3.1 Mathematical Definition of $\mathtt{VRPTW}(\tau;\theta)$

In this section, the time windows $\tau$ and operational parameters $\theta$
are assumed to be fixed to certain values in their domain and support
respectively. As far as the routing operation is concerned, we shall assume
that each customer with non-zero demand must be visited exactly once by a
single vehicle; that is, split deliveries are not allowed. In this regard, a
_route set_ $\mathbf{R}=(R_{1},\ldots,R_{m})$, where $m\geq 1$, represents a
partition of the customer set $V_{C}$. Here,
$R_{k}=(R_{k,1},\ldots,R_{k,n_{k}})$ represents the $k^{\text{th}}$ vehicle
route, $R_{k,l}$ the $l^{\text{th}}$ customer and $n_{k}$ the number of
customers visited by vehicle $k$. The cost of a route $R_{k}$ is evaluated as
$c(R_{k})=\sum_{l=0}^{n_{k}}c_{R_{k,l},R_{k,l+1}}$, where we define
$R_{k,0}=R_{k,n_{k}+1}=0$, and the cost of $\mathbf{R}$ is defined as
$c(\mathbf{R})=\sum_{k=1}^{m}c(R_{k})$. The route set $\mathbf{R}$ is
feasible, if (i) all capacity constraints are satisfied, i.e., $\sum_{i\in
R_{k}}q_{i}\leq Q$ for all $k\in\\{1,\ldots,m\\}$, and (ii) all time window
constraints are satisfied, i.e., there exists a vector of arrival times,
$\bm{a}\in\mathcal{X}(\mathbf{R},\tau;\theta)$, where
$\mathcal{X}(\mathbf{R},\tau;\theta)$ is the feasible solution set of the
following linear system of inequalities:

$\mathcal{X}(\mathbf{R},\tau;\theta)=\begin{aligned}
\left\\{\bm{a}\in\mathbb{R}_{\geq 0}^{n}\left|\begin{aligned} &a_{R_{k,1}}\geq
e_{0}+t_{0,R_{k,1}}&&\forall k\in K\coloneqq\\{1,\ldots,m\\}\\\
&a_{R_{k,l+1}}-a_{R_{k,l}}\geq t_{R_{k,l},R_{k,l+1}}+u_{R_{k,l}}&&\forall
l\in\left\\{1,\ldots,n_{k}-1\right\\},\forall k\in K\\\
&a_{R_{k,n_{k}}}\leq\ell_{0}-t_{R_{k,n_{k}},0}-u_{R_{k,n_{k}}}&&\forall k\in
K\\\ &a_{i}\in\tau_{i}&&\forall i\in
V_{C}\end{aligned}\right.\right\\}\end{aligned}$ (2)

In this definition, $a_{R_{k,l}}$ denotes the arrival time at location
$R_{k,l}$, i.e., the arrival time at the $l^{\text{th}}$ location on the
$k^{\text{th}}$ vehicle route. The first three inequalities essentially
require that the arrival time at any location must be at least as large as the
sum of the arrival and service times in the previous location and the time to
travel from the previous to the current location. The last inequality requires
the arrival time at customer location $i\in V_{C}$ to be within the time
window $\tau_{i}$. Observe that, by this definition, if a vehicle arrives at
customer location $i\in V_{C}$ at a time earlier than
$\underaccent{\bar}{\tau}_{i}\coloneqq\min_{t\in\tau_{i}}t$, then it is
allowed to wait until $\underaccent{\bar}{\tau}_{i}$. However, arriving later
than $\bar{\tau}_{i}\coloneqq\max_{t\in\tau_{i}}t$ is not permitted. We denote
by $\mathcal{R}(\tau;\theta)$ the set of all feasible route sets for the given
realization of operational parameters $\theta$ and time window assignment
$\tau$. The value of $\mathtt{VRPTW}(\tau;\theta)$ can now be defined as the
optimal value of the following optimization problem:

$\begin{array}[]{r@{\;\;}l}\displaystyle\mathop{\text{minimize}}_{\mathbf{R}}&c\left(\mathbf{R}\right)\\\
\text{subject to}&\mathbf{R}\in\mathcal{R}(\tau;\theta).\end{array}$
($\mathtt{VRPTW}(\tau;\theta)$)

###  3.2 Deterministic Equivalent of Stochastic Programming Formulation

In practice, the exact joint probability distribution $\mathbb{P}$ is either
not explicitly available or is hard to obtain. Indeed, even if it is known
exactly, computing the objective function involves multi-dimensional
integration of the function $\mathtt{VRPTW}(\tau;\cdot)$, which is practically
impossible considering that the solution of the deterministic problem
$\mathtt{VRPTW}(\tau;\theta)$ is itself challenging and cannot be obtained in
closed form. Instead, we assume that we are given a finite set of $S$
scenarios $\theta_{1},\ldots,\theta_{S}$ along with associated probabilities
of occurrence $p_{1},\ldots,p_{S}$, where $p_{s}>0$,
$s\in\mathcal{S}\coloneqq\\{1,\ldots,S\\}$ and
$\sum_{s\in\mathcal{S}}p_{s}=1$. In this situation, we seek to optimize the
following _deterministic equivalent_ of problem (1), obtained by replacing the
expectation with a sample average:

$\begin{array}[]{r@{\;\;}l}\displaystyle\mathop{\text{minimize}}_{\tau}&\displaystyle\sum_{s\in\mathcal{S}}p_{s}\mathtt{VRPTW}(\tau;\theta_{s})\\\
\text{subject to}&\tau_{i}\in TW_{i}\quad\forall i\in V_{C}.\end{array}$ (3)

Following the definition of $\mathtt{VRPTW}(\tau;\theta)$, the above sample
average formulation (3) can be equivalently represented as follows:

$\begin{array}[]{r@{\;\;}l}\displaystyle\mathop{\text{minimize}}_{\tau,\mathbf{R}}&\displaystyle\sum_{s\in\mathcal{S}}p_{s}c\left(\mathbf{R}_{s}\right)\\\
\text{subject to}&\tau_{i}\in TW_{i}\quad\forall i\in V_{C}\\\
&\mathbf{R}_{s}\in\mathcal{R}(\tau;\theta_{s})\quad\forall
s\in\mathcal{S}.\end{array}$ (4)

The optimization problem (4) shall be our primary focus for the rest of the
paper. We shall denote by $(\tau,\\{{\mathbf{R}_{s}}\\}_{s\in\mathcal{S}})$ a
feasible solution to this problem.

##  4 Solution Approach

Our solution approach for the TWAVRP is motivated by the observation that, in
the continuous setting (where $V_{\text{disc}}=\emptyset$), problem (4) can be
reduced to an instance of the arrival-time ConVRP. Consequently, any algorithm
to solve the latter class of problems can be used to solve problem (4).
Section 4.1 presents an exact branch-and-bound algorithm for this purpose; we
note, however, that the presented algorithm is more general-purpose, since it
can also address the setting where $V_{\text{disc}}\neq\emptyset$. Section 4.2
presents new valid disjunctions that can be used as alternative branching
rules in the algorithm; Section 4.3 presents upper bounding procedures (i.e.,
generating good time window assignments) in the context of our algorithm; and
finally, Section 4.4 shows how the (exact) algorithm can be modified as a
heuristic to solve large-scale instances.

###  4.1 Overview of Exact Algorithm

We adapt the decomposition algorithm of [35] developed for the Consistent
Traveling Salesman Problem (the single-vehicle variant of the ConVRP) to solve
the TWAVRP to optimality. Notably, we extend the algorithm to incorporate also
the case of discrete time windows. This section presents the main ingredients
of the algorithm translated into the TWAVRP context.

The algorithm uses a branch-and-bound tree search to identify the optimal time
window assignments by solving within each node a set of VRPTW instances. The
tree is initialized with the original problem instance enforcing only the
exogenous time windows $[e,\ell]$. If valid time window assignments $\tau\in
TW$ cannot be constructed using the optimal solution of the current node, the
algorithm creates new nodes by using disjunctions (5a) and (5b) as branching
rules. The resulting branching rules are valid because, for every feasible
solution $(\tau,\\{\mathbf{R}_{s}\\}_{s\in\mathcal{S}})$ in problem (4), there
exist arrival-time vectors
$\bm{a}_{s}\in\mathcal{X}(\mathbf{R}_{s},[e,\ell];\theta_{s})$ for each
$s\in\mathcal{S}$ that satisfy disjunctions (5).

$\displaystyle\left[a_{si}\leq\beta+w_{i}/2\;\forall
s\in\mathcal{S}\right]\;\vee\;$ $\displaystyle\left[a_{si}\geq\beta-
w_{i}/2\;\forall s\in\mathcal{S}\right]$
$\displaystyle\qquad\forall\beta\in\mathbb{R},\;\forall i\in V_{\text{cont}}$
(5a) $\displaystyle\left[a_{si}\leq\bar{y}_{ib}\;\forall
s\in\mathcal{S}\right]\;\vee\;$
$\displaystyle\left[a_{si}\geq\underaccent{\bar}{y}_{i,b+1}\;\forall
s\in\mathcal{S}\right]$ $\displaystyle\qquad\forall
b\in\\{1,\ldots,N_{i}-1\\},\;\forall i\in V_{\text{disc}}.$ (5b)

Conversely, if there exist route sets
$\mathbf{R}_{s}\in\mathcal{R}([e,\ell];\theta_{s})$ and arrival-time vectors
$\bm{a}_{s}\in\mathcal{X}(\mathbf{R}_{s},[e,\ell];\theta_{s})$ for each
$s\in\mathcal{S}$ satisfying disjunctions (5), then there exists a time window
assignment $\tau\in TW$ such that
$(\tau,\\{\mathbf{R}_{s}\\}_{s\in\mathcal{S}})$ is feasible in problem (4).
Specifically, a feasible time window assignment is

$\tau_{i}^{\star}=\begin{dcases}\left[y_{i},y_{i}+w_{i}\right],\text{where
}y_{i}=\min\left\\{\ell_{i}-w_{i},\min_{s\in\mathcal{S}}a_{si}\right\\}&\text{if
}i\in V_{\text{cont}}\\\
\left[\underaccent{\bar}{y}_{ib_{i}},\bar{y}_{ib_{i}}\right],\text{where
}b_{i}\in\mathop{\arg\min}_{b\in\\{1,\ldots,N_{i}\\}}\left\\{\bar{y}_{ib}:\bar{y}_{ib}\geq\max_{s\in\mathcal{S}}a_{si}\right\\}&\text{if
}i\in V_{\text{disc}}\end{dcases}\quad\forall i\in V_{C}.$ (6)

For given route sets $\mathbf{R}_{s}\in\mathcal{R}([e,\ell];\theta_{s})$,
verifying the existence of arrival-time vectors
$\bm{a}_{s}\in\mathcal{X}(\mathbf{R}_{s},[e,\ell];\theta_{s})$,
$s\in\mathcal{S}$, that satisfy disjunctions (5) is equivalent to verifying
that the optimal objective value $\delta$ of the following mixed-integer
linear optimization problem is non-positive.

$\begin{array}[]{l@{\quad}ll}\displaystyle\mathop{\text{minimize}}_{\delta,a,z,\bar{\mu},\underaccent{\bar}{\mu}}&\delta\\\
\text{subject to}&\displaystyle\delta\in\mathbb{R}_{\geq
0},\;\;\bm{a}_{s}\in\mathcal{X}(\mathbf{R}_{s},\tau;\theta_{s}),\;s\in\mathcal{S}\\\
&\left.\begin{array}[]{l}\mspace{-11.0mu}\displaystyle\bar{\mu}_{i},\underaccent{\bar}{\mu}_{i}\in\mathbb{R}_{\geq
0},\;\;z_{ib}\in\\{0,1\\},\;b\in\\{1,\ldots,N_{i}\\}\\\
\mspace{-11.0mu}\displaystyle\sum_{b=1}^{N_{i}}z_{ib}=1\\\
\mspace{-11.0mu}\displaystyle z_{ib}=1\;\;\Rightarrow\;\;\bar{\mu}_{i}\geq
a_{si}-\bar{y}_{ib}\quad\forall(s,b)\in\mathcal{S}\times\\{1,\ldots,N_{i}\\}\\\
\mspace{-11.0mu}\displaystyle
z_{ib}=1\;\;\Rightarrow\;\;\underaccent{\bar}{\mu}_{i}\geq\underaccent{\bar}{y}_{ib}-a_{si}\quad\forall(s,b)\in\mathcal{S}\times\\{1,\ldots,N_{i}\\}\end{array}\right\\}&\forall
i\in V_{\text{disc}}\\\ &\displaystyle\delta\geq
a_{s_{1}i}-a_{s_{2}i}-w_{i}\quad\forall(s_{1},s_{2})\in\mathcal{S}\times\mathcal{S}&\forall
i\in V_{\text{cont}}\\\ &\displaystyle\delta\geq\sum_{i\in
V_{\text{disc}}}\left(\bar{\mu}_{i}+\underaccent{\bar}{\mu}_{i}\right).\end{array}$
(7)

In this problem, $\delta$ records the minimum possible violation of
disjunctions (5) across all feasible arrival-time vectors
$\bm{a}_{s}\in\mathcal{X}(\mathbf{R}_{s},[e,\ell];\theta_{s})$,
$s\in\mathcal{S}$. The second-to-last constraint ensures that $\delta$ is at
least as large as the maximum violation across all members of (5a), while the
last constraint ensures that $\delta$ is at least as large as the sum of
violations across all members of (5b). In the latter case, the binary variable
$z_{ib}$ indicates whether the $b^{\text{th}}$ member of (5b) is minimally
violated for given $i\in V_{\text{disc}}$. In other words, if $z_{ib}=1$, then
$\left[\underaccent{\bar}{y}_{ib},\bar{y}_{ib}\right]$ is the best time window
for customer $i$, and
$\bar{\mu}_{i}=\left[\max_{s\in\mathcal{S}}a_{si}-\bar{y}_{ib}\right]_{+}$ and
$\underaccent{\bar}{\mu}_{i}=\left[\underaccent{\bar}{y}_{ib}-\min_{s\in\mathcal{S}}a_{si}\right]_{+}$
respectively record the arrival-time violations with respect to the start and
end of this time window (refer to Figure 2). Here,
$[\cdot]_{+}=\max\\{\cdot,0\\}$.

time$\vphantom{\ell_{i}}e_{i}$$\vphantom{\ell_{i}}a_{1i}$$\vphantom{\ell_{i}}a_{2i}$$\vphantom{\ell_{i}}a_{3i}$$\ell_{i}$$\leq\delta+w_{i}$time$\vphantom{\ell_{i}}e_{i}$$\vphantom{\ell_{i}}a_{1i}$$\vphantom{\ell_{i}}\underaccent{\bar}{y}_{ib}$$\vphantom{\ell_{i}}\bar{y}_{ib}$$\vphantom{\ell_{i}}a_{2i}$$\vphantom{\ell_{i}}a_{3i}$$\ell_{i}$$\underaccent{\bar}{\mu}_{i}$$\bar{\mu}_{i}$
Figure 2: Decision variables in problem (7) for the continuous (left) and
discrete (right) cases.

#### Algorithm.

The above observations suggest that we can solve the TWAVRP without having to
explicitly encode its time window assignments. In fact, any TWAVRP instance
decomposes into its individual scenarios if we track, within each node of a
branch-and-bound search tree, a vector of applicable time windows (one for
each customer), which we shall denote by $\tau$. The root node enforces only
the exogenous time windows (see Step 1). Processing a node amounts to solving
a set of VRPTW instances (one for each scenario) with time window constraints
enforced by that node (see Step 3). It is important to remark that these VRPTW
instances are uncoupled, and can thus be solved independently of each other.
This is because none of the expressions within each disjunct in (5) (upon
which our branching rules are based) link arrival-times from different
scenarios in the same inequality. The resulting optimal route sets
$\\{\mathbf{R}_{s}\\}_{s\in\mathcal{S}}$ are then used as inputs to the
separation problem (7) (see Step 5). If the optimal objective value satisfies
$\delta\leq 0$, then a new, improved time window assignment $\tau^{\star}$ is
recorded, as per (6). Otherwise, the algorithm creates two new nodes with
tightened time windows for customer $i^{\star}\in V_{C}$ (see Step 6). The
overall algorithm is as follows.

  1. 1.

Initialize. Set root node
$\tau^{0}\leftarrow\left(\left[e_{1},\ell_{1}\right],\ldots,\left[e_{n},\ell_{n}\right]\right)$,
node queue $\mathcal{N}\leftarrow\left\\{\tau^{0}\right\\}$, upper bound
$UB\leftarrow+\infty$ and optimal time window assignment
$\tau^{\star}\leftarrow\emptyset$.

  2. 2.

Check convergence. If $\mathcal{N}=\emptyset$, then stop: $\tau^{\star}$ is
the optimal time window assignment with (expected) cost $UB$. Otherwise,
select a node $\tau$ from $\mathcal{N}$, and set
$\mathcal{N}\leftarrow\mathcal{N}\setminus\\{\tau\\}$.

  3. 3.

Process node. For each $s\in\mathcal{S}$, solve
$\mathtt{VRPTW}(\tau;\theta_{s})$; and let $\mathbf{R}_{s}$ denote an optimal
solution.

  4. 4.

Fathom by bound. If $\sum_{s\in\mathcal{S}}p_{s}c(\mathbf{R}_{s})\geq UB$,
then go to Step 2.

  5. 5.

Check feasibility. Let $(\delta,a,z,\bar{\mu},\underaccent{\bar}{\mu})$ be an
optimal solution of problem (7). If $\delta\leq 0$, then a new, improved time
window assignment is found: set $\tau^{\star}$ as per (6), set
$UB\leftarrow\sum_{s\in\mathcal{S}}p_{s}c(\mathbf{R}_{s})$ and go to Step 2.

  6. 6.

Branch. Instantiate two children nodes, $\tau^{\text{L}}$ and
$\tau^{\text{R}}$, from the parent node: $\tau^{\text{L}}\leftarrow\tau$,
$\tau^{\text{R}}\leftarrow\tau$. If $\delta>\sum_{i\in
V_{\text{disc}}}(\bar{\mu}_{i}+\underaccent{\bar}{\mu}_{i})$, then do Step 6a;
otherwise, do Step 6b:

    1. (a)

Branch as per (5a). Let $i^{\star}\in V_{\text{cont}}$ be any customer for
which
$\delta=\max_{s\in\mathcal{S}}a_{si^{\star}}-\min_{s\in\mathcal{S}}a_{si^{\star}}-w_{i^{\star}}$,
and let
$\beta^{\star}=\left(\max_{s\in\mathcal{S}}a_{si^{\star}}+\min_{s\in\mathcal{S}}a_{si^{\star}}\right)/2$.
Tighten the time window for $i^{\star}$ as follows: (i)
$\tau^{\text{L}}_{i^{\star}}\leftarrow\left[\min_{t\in\tau_{i^{\star}}^{\text{L}}}t,\beta^{\star}+\frac{1}{2}w_{i^{\star}}\right]$,
(ii)
$\tau^{\text{R}}_{i^{\star}}\leftarrow\left[\beta^{\star}-\frac{1}{2}w_{i^{\star}},\max_{t\in\tau_{i^{\star}}^{\text{R}}}t\right]$.

    2. (b)

Branch as per (5b). Let $i^{\star}\in V_{\text{disc}}$ be any member of
$\mathop{\arg\max}_{i\in
V_{\text{disc}}}\\{\bar{\mu}_{i}+\underaccent{\bar}{\mu}_{i}\\}$. If
$\bar{\mu}_{i^{\star}}\geq\underaccent{\bar}{\mu}_{i^{\star}}$, let
$b^{\star}=\sum_{b=1}^{N_{i}}b\mathbbm{1}_{\left[z_{ib}=1\right]}$; else, let
$b^{\star}=\sum_{b=1}^{N_{i}}(b-1)\mathbbm{1}_{\left[z_{ib}=1\right]}$.
Tighten the time window for $i^{\star}$ as follows: (i)
$\tau^{\text{L}}_{i^{\star}}\leftarrow\left[\min_{t\in\tau_{i^{\star}}^{\text{L}}}t,\bar{y}_{i^{\star}b^{\star}}\right]$,
(ii)
$\tau^{\text{R}}_{i^{\star}}\leftarrow\left[\underaccent{\bar}{y}_{i^{\star}\,b^{\star}+1},\max_{t\in\tau_{i^{\star}}^{\text{R}}}t\right]$.

Set
$\mathcal{N}\leftarrow\mathcal{N}\cup\\{\tau^{\text{L}},\tau^{\text{R}}\\}$,
and go to Step 2.

We remark that any node selection strategy can be used in Step 2 to guarantee
convergence.

An illustration on a small example is now presented to aid understanding and
give intuition about the algorithm. Consider the TWAVRP instance shown in
Figure 3. This example features $n=4$ customers, with
$V_{\text{cont}}=\\{1,2,3\\}$ and $V_{\text{disc}}=\\{4\\}$. Only customer
demands are uncertain and they are represented using $S=2$ scenarios.

301244226547$\left[0,15\right]$$\left[0,30\right]$$\left[0,10\right]$$\left[0,20\right]$$\left[0,7\right],\left[14,22\right]$

$(n,S,Q)=(4,2,3)$ $p_{s}$ $q_{s1}$ $q_{s2}$ $q_{s3}$ $q_{s4}$ $s=1\quad$ 1/2 1
3 1 1 $s=2\quad$ 1/2 3 1 1 1 $w_{1}=w_{2}=w_{3}=3$

Figure 3: Instance parameters for the illustrative example. The arc weights
denote both travel times and costs, and they are such that the triangle
inequality is satisfied. Note that the graph is assumed to be complete, with
nodes 0–3 lying along a straight line and costs/times to travel along this
line being cumulative (not all arcs are shown for convenience). All service
times are zero. The depicted intervals denote exogenous time windows for nodes
1–3 and feasible ones for node 4.

The search tree of our algorithm to solve the illustrative example of Figure 3
is shown in Figure 4. Each “rectangle” denotes a node of our search tree.
Within each rectangle, for each scenario $s\in\\{1,2\\}$, the optimal route
set $\mathbf{R}_{s}$ is shown. The $x$-coordinate of each customer $i\in
V_{C}$ denotes its arrival-time $a_{si}$, which is computed by solving the
separation problem (7) to optimality (note that there may be multiple optimal
solutions in each case). Finally, on top of each rectangle, the time windows
$\tau$ enforced by our algorithm and the objective value of the corresponding
optimal route sets (equal to $\sum_{s\in\mathcal{S}}p_{s}c(\mathbf{R}_{s})$)
are shown. In the root node, the separation problem (7) certifies (in Step 5)
that no valid time window assignment can be constructed from its optimal
solution. In particular, if we focus on customer $3\in V_{\text{cont}}$, then
$a_{13}=14$ and $a_{23}=5$. These arrival times clearly do not fall within a
time window of width $w_{3}=3$. Therefore, as per Step 6a,
$\beta^{\star}=(14+5)/2=9.5$. The time window of customer $3$ in the left
child is tightened to $[e_{3},\beta^{\star}+w_{3}/2]=[0,11]$, while in the
right child to $[\beta^{\star}-w_{3}/2,\ell_{3}]=[8,15]$. Similarly, if we
focus on customer $4\in V_{\text{disc}}$ in this right child node, then
$a_{14}=9$ and $a_{24}=22$. These arrival times also do not simultaneously
satisfy either of the two candidate windows, $[0,7]$ or $[14,22]$. Therefore,
a branch is made, as per Step 6b, tightening the time window of customer $4$
in the left child to $[0,7]$ and in the right one to $[14,22]$.

0510152025$s=2$$s=1$12341234Branch as per (5a):
$(i^{\star},\beta^{\star})=(3,9.5)$$\tau=\big{(}[0,10],\,[0,20],\,[0,15],\,[0,22]\big{)}$,
obj =
25.50510152025$s=2$$s=1$12341234$\tau_{3}^{\star}$$\tau_{1}^{\star}$$\leftarrow\tau_{4}^{\star}$$\tau_{2}^{\star}$(optimal)$\tau=\big{(}[0,10],\,[0,20],\,[0,{11}],\,[0,22]\big{)}$,
obj = 26.00510152025$s=2$$s=1$12341234Branch as per (5b):
$(i^{\star},b^{\star})=(4,1)$$\tau=\big{(}[0,10],\,[0,20],\,[{8},15],\,[0,22]\big{)}$,
obj =
26.00510152025$s=2$$s=1$12341234$\leftarrow\tau_{4}^{\star}$$\tau_{3}^{\star}$$\tau_{1}^{\star}$$\tau_{2}^{\star}$(feasible)$\tau=\big{(}[0,10],\,[0,20],\,[8,15],\,[0,{7}]\big{)}$,
obj =
26.50510152025$s=2$$s=1$12341234$\tau_{3}^{\star}$$\tau_{1}^{\star}$$\leftarrow\tau_{4}^{\star}$$\tau_{2}^{\star}$(feasible)$\tau=\big{(}[0,10],\,[0,20],\,[8,15],\,[{14},22]\big{)}$,
obj = 26.5$a_{3}\leq 11$$a_{3}\geq 8$$a_{4}\leq 7$$a_{4}\geq 14$ Figure 4: The
search tree of our algorithm for the illustrative example of Figure 3.

We remark that, in a given node of our search tree (except the root node), one
does not need to solve a VRPTW subproblem for every scenario as required by
Step 3 of the algorithm, and the optimal VRPTW route sets for some of the
scenarios can be directly transferred from the parent subproblems. In fact,
after branching has occurred in Step 6 of the algorithm, at most $S$ (out of a
total of $2S$) VRPTW subproblems have to be solved across the two children
nodes. This is because, by construction, any arrival time vectors
corresponding to the optimal solution, $\mathbf{R}_{s}$, of a given scenario
$s$ cannot simultaneously violate both disjuncts of the applied branching
disjunction (either (5a) or (5b)), although it may satisfy both disjuncts
simultaneously. Therefore, as far as a given scenario $s$ is concerned, a
VRPTW subproblem needs to be solved only at most once across the two children
nodes. To illustrate this, see Figure 4. After branching has occurred in the
root node, $\mathbf{R}_{1}$ in the right child node is exactly the same as
that in the parent, since it already satisfies the applied branching
constraint $[a_{3}\geq 8]$. For the same reason, $\mathbf{R}_{2}$ in the left
child node is exactly the same as that in the parent. Once the branching rule
has been established, whether a scenario-specific set of routes remains
feasible (and hence, optimal) for the VRPTW instance of a child node can be
inferred trivially by inspection, and hence, the corresponding instance need
not be solved, as the routes can be copied over. For the VRPTW instances that
indeed warrant a new route set to be computed, Section 5 describes methods for
solving the corresponding subproblems.

###  4.2 Path-based Disjunctions

The algorithm described in the previous section converges in finite time (the
argument is similar to that in [35]). However, the branching Step 6 may not
necessarily “cut off” the VRPTW route set
$\\{\mathbf{R}_{s}\\}_{s\in\mathcal{S}}$ found in the parent node. Indeed, it
is only guaranteed to cut off the arrival time vector
$\\{\bm{a}_{s}\\}_{s\in\mathcal{S}}$ corresponding to
$\\{\mathbf{R}_{s}\\}_{s\in\mathcal{S}}$. More specifically, an optimal route
set of the (right) child node $\tau^{\text{R}}$ may also be optimal for its
parent node $\tau$. This is because the only difference between
$\tau^{\text{R}}$ and its parent $\tau$ is that the former features a tighter
earliest start time for customer $i^{\star}\in V_{C}$. It is possible for an
optimal route set of $\tau$ to be unchanged, if the tighter earliest start
time constraint can be satisfied by simply allowing the vehicle to wait longer
at $i^{\star}$. To see this, consider again the example of Figure 3. If the
exogenous time window of customer $2$ is increased by one unit to $[0,21]$,
then an optimal route set $\mathbf{R}_{2}$ in the root node of the algorithm
(see Figure 4) will also be optimal in its right child node since the vehicle
visiting customer $3$ will simply wait longer until it satisfies the branching
constraint, $a_{3}\geq 8$. Along with the fact that $\mathbf{R}_{1}$ is also
the same (refer to the discussion in the previous section), this means that
the computed VRPTW route set $\\{\mathbf{R}_{s}\\}_{s\in\mathcal{S}}$ in the
root node has not changed in its right child node. The impact of this is non-
improving lower bounds: the objective value of the node $\tau^{\text{R}}$ will
be exactly the same as that of its parent, leading to slow convergence and
poor numerical performance. This observation motivates us to investigate
branching rules which are guaranteed to cut off the parent route set.

Our motivation for the new class of disjunctions comes from the _path
precedence inequalities_ proposed in [10] for the continuous TWAVRP and the
_inconsistent path elimination constraints_ proposed in [34] for the ConTSP.
Consider a feasible solution $(\tau,\\{\mathbf{R}_{s}\\}_{s\in\mathcal{S}})$
to problem (4). Suppose that there is a vehicle route in the solution
$\mathbf{R}_{s_{1}}$ of scenario $s_{1}\in\mathcal{S}$ in which customer $i\in
V_{C}$ is visited before customer $j\in V_{C}\setminus\\{i\\}$, and that there
is a vehicle route in the solution $\mathbf{R}_{s_{2}}$ of scenario
$s_{2}\in\mathcal{S}\setminus\\{s_{1}\\}$ in which $j$ is visited before $i$.
Since both $i$ and $j$ are visited within their respective time windows
$\tau_{i}$ and $\tau_{j}$ in both scenarios, it must be the case that the sum
of the travel times from $i$ to $j$ in scenario $s_{1}$ and from $j$ to $i$ in
scenario $s_{2}$ is at most the sum of the widths of their time windows,
$w_{i}+w_{j}$, as shown in Figure 5. Consequently, if this condition is not
satisfied by a route set $\\{\mathbf{R}_{s}\\}_{s\in\mathcal{S}}$ for all
possible pairs $(i,j)$, then there cannot exist a feasible time window
assignment $\tau\in TW$ such that
$(\tau,\\{\mathbf{R}_{s}\\}_{s\in\mathcal{S}})$ is feasible in problem (4).

time$a_{s_{2}j}$$a_{s_{1}i}$$a_{s_{2}i}$$a_{s_{1}j}$width = $w_{i}$width =
$w_{j}$$(a_{s_{2}i}-a_{s_{2}j})$$(a_{s_{1}j}-a_{s_{1}i})$$j$$i$$i$$j$scenario
$s_{1}$scenario $s_{2}$ Figure 5: Motivation for the path-based disjunctions.
The above shown sub-paths must satisfy
$(a_{s_{2}i}-a_{s_{2}j})+(a_{s_{1}j}-a_{s_{1}i})\leq w_{i}+w_{j}$, if they are
to be part of a feasible solution.

To construct valid disjunctions based on the above observation, we first
introduce some notation. Let $\pi=(v_{1},\ldots,v_{p})$ denote a directed
$v_{1}-v_{p}$ path in graph $G$ that is formed by the arcs in the set
$\\{(v_{i},v_{i+1}):i=1,\ldots,p-1\\}$, where $(v_{i},v_{i+1})\in A$ for all
$i=1,\ldots,p-1$. We shall only consider paths which are open and simple,
i.e., $p>1$ and $v_{i}\neq v_{j}$ for $i\neq j$. For a given realization of
the travel and service times, the travel time along $\pi$ is defined to be
$t(\pi)=\sum_{i=1}^{p-1}\left(t_{v_{i}v_{i+1}}+u_{v_{i}}\right)$, where we
define $u_{0}=0$. Note that as per this definition, the travel time along a
path does not include any waiting time that might potentially be incurred at
its nodes. Finally, a route set $\mathbf{R}=(R_{1},\ldots,R_{m})$, where
$R_{k}=(R_{k,1},\ldots,R_{k,n_{k}})$ for $k=1,\ldots,m$, is said to contain
$\pi$, if $\pi$ appears as a sub-path in any of its routes, i.e., if for some
$k$, we have
$R_{k}=(R_{k,1},\ldots,R_{k,l}=v_{1},R_{k,l+1}=v_{2},\ldots,R_{k,l+p-1}=v_{p},\ldots,R_{k,n_{k}})$.
The key result of this section is that every feasible solution
$(\tau,\\{{\mathbf{R}_{s}}\\}_{s\in\mathcal{S}})$ to problem (4) satisfies the
following disjunctions:

$\displaystyle\left[\begin{array}[]{c}\mathbf{R}_{s}\text{ does not contain
any }i-j\\\ \text{path with travel time}\geq d_{1}\end{array}\;\;\forall
s\in\mathcal{S}\right]\vee\left[\begin{array}[]{c}\mathbf{R}_{s}\text{ does
not contain any }j-i\\\ \text{path with travel time}\geq
d_{2}\end{array}\;\;\forall s\in\mathcal{S}\right]$ (8)
$\displaystyle\forall(d_{1},d_{2})\in\mathbb{R}^{2}:d_{1}+d_{2}>w_{i}+w_{j},\;\;\forall(i,j)\in
V_{C}\times V_{C}:i\neq j,$

where $w_{k}$ for any $k\in V_{\text{disc}}$ is defined to be
$w_{k}=\max_{b\in\\{1,\ldots,N_{k}\\}}\\{\bar{y}_{kb}-\underaccent{\bar}{y}_{kb}\\}$.
However, the converse is not true. To see this, consider the following
counter-example.

  * •

$(n,S,Q)=(4,2,3)$. $V_{\text{cont}}=V_{C}$ and $w_{i}=1$ for all $i\in V_{C}$.

  * •

$[e,\ell]=\left([0,6],[0,6],[3,4],[4,5]\right)$. Also,
$[e_{0},\ell_{0}]=[0,10]$.

  * •

$G=(V,A)$ is complete. $c_{ij}=t_{ij}=1$ for all $(i,j)\in A$ and $u_{i}=0$
for all $i\in V_{C}$.

  * •

Demand is uncertain. $(q_{s1},q_{s2},q_{s3},q_{s4})=(1,1,1,3)$ for $s=1$ and
$(1,1,3,1)$ for $s=2$.

Consider the route sets $\\{\mathbf{R}_{s}\\}_{s=1,2}$ shown in Figure 6. The
$x$-coordinates correspond to arrival times. This solution satisfies all path-
based disjunctions (8). However, it is not a feasible TWAVRP solution since
there are no valid time window assignments $\tau\in TW$ for customers $1$ and
$2$. In contrast, observe that this solution does indeed violate the
(necessary and sufficient) time window-based disjunctions (5a) corresponding
to $(i,\beta)=(1,4)$ as well as $(i,\beta)=(2,4)$.

02468$s=2$$s=1$12341234 Figure 6: Counter-example to show that the
disjunctions (8) are not sufficient for the TWAVRP.

The above observations suggest that we can use the path-based disjunctions (8)
as the basis of a branching rule in addition to the time window-based
disjunctions (5). The corresponding branching constraints, i.e., constraints
within each individual disjunct of (8), are equivalent to _path elimination
constraints_ (e.g., see [3]). Consequently, the subproblems to be solved in
Step 3 of the algorithm are VRPTW instances with additional path elimination
constraints. To describe these subproblems formally, if $\mathcal{F}$ is a
finite collection of triples $(i,j,d)\in V_{C}\times V_{C}\times\mathbb{R}$,
such that $i\neq j$, then we let each member of $\mathcal{F}$ represent a
family of _forbidden paths_. Specifically, the member $(i,j,d)\in\mathcal{F}$
represents the set of all $i-j$ paths with travel time greater than or equal
to $d$. Given a collection $\mathcal{F}$ of representative forbidden paths,
time windows $\tau$ and operational parameters $\theta$, we define
$\mathtt{VRPTWFP}(\tau,\mathcal{F};\theta)$ to be the optimal value of the
following optimization problem:

$\begin{array}[]{r@{\;\;}l}\displaystyle\mathop{\text{minimize}}_{\mathbf{R}}&c\left(\mathbf{R}\right)\\\
\text{subject to}&\mathbf{R}\in\mathcal{R}(\tau;\theta)\\\ &\mathbf{R}\text{
does not contain any $i$-$j$ path }\pi\text{ with }t(\pi)\geq
d\quad\forall(i,j,d)\in\mathcal{F}.\end{array}$
($\mathtt{VRPTWFP}(\tau,\mathcal{F};\theta)$)

Note that
$\mathtt{VRPTWFP}(\tau,\emptyset;\theta)=\mathtt{VRPTW}(\tau;\theta)$ so our
definition is consistent with the one in Section 3.1.

Before we can incorporate the path-based disjunctions (8) as branching rules
in our algorithm, we also need a separation algorithm, which will take as
input route sets $\\{\mathbf{R}_{s}\\}_{s\in\mathcal{S}}$, and return either a
violated member of (8) or a certificate that all of its members are satisfied.
We compute the following quantities, where we assume $(i,j)\in V_{C}\times
V_{C}$, such that $i\neq j$, is a given pair of customers.

  * •

$\mathcal{S}_{ij}$: scenarios containing an $i-j$ path; that is,
$\mathcal{S}_{ij}=\left\\{s\in\mathcal{S}:\mathbf{R}_{s}\text{ contains an
}i-j\text{ path}\right\\}$.

  * •

$d_{ij}^{s}$: travel time of $i-j$ path in $\mathbf{R}_{s}$, where
$s\in\mathcal{S}_{ij}$.

  * •

$\nu_{ij}$: no. of violating scenario pairs; that is,
$\nu_{ij}=\left\lvert\left\\{(s_{1},s_{2})\in\mathcal{S}_{ij}\times\mathcal{S}_{ji}:d_{ij}^{s_{1}}+d_{ji}^{s_{2}}>w_{i}+w_{j}\right\\}\right\rvert$.

  * •

$\Delta_{ij}$: minimum value of sum of path travel times (across violating
scenario pairs); that is,
$\Delta_{ij}=\inf_{(s_{1},s_{2})\in\mathcal{S}_{ij}\times\mathcal{S}_{ji}}\left\\{d_{ij}^{s_{1}}+d_{ji}^{s_{2}}:d_{ij}^{s_{1}}+d_{ji}^{s_{2}}>w_{i}+w_{j}\right\\}$.

We are now in a position to incorporate the path-based-disjunctions (8) in our
algorithm. To do so, we store the set of forbidden paths $\mathcal{F}$ as part
of a node’s characteristic data (along with $\tau$) at the time of node
creation (initialization and branching steps), and we use this set as input to
the VRPTW subproblem at the time of node processing. We update $\mathcal{F}$
in the branching step based on a ranked list $\mathcal{L}$ of pairs $(i,j)$
for which the corresponding path-based disjunctions (8) can be used as
branching rules. In all, the following modifications are made to the main
algorithm.

  1. 1′

Initialize. Set root node
$\tau^{0}\leftarrow\left(\left[e_{1},\ell_{1}\right],\ldots,\left[e_{n},\ell_{n}\right]\right)$,
$\mathcal{F}^{0}\leftarrow\emptyset$, node queue
$\mathcal{N}\leftarrow\left\\{(\tau^{0},\mathcal{F}^{0})\right\\}$, upper
bound $UB\leftarrow+\infty$ and optimal time window assignment
$\tau^{\star}\leftarrow\emptyset$.

  2. 3′

Process node. For each $s\in\mathcal{S}$, solve
$\mathtt{VRPTWFP}(\tau,\mathcal{F};\theta_{s})$; and let $\mathbf{R}_{s}$
denote its optimal solution.

  3. 6′

Branch. Instantiate two children nodes from $(\tau,\mathcal{F})$:
$(\tau^{\text{L}},\mathcal{F}^{\text{L}})\leftarrow(\tau,\mathcal{F})$,
$(\tau^{\text{R}},\mathcal{F}^{\text{R}})\leftarrow(\tau,\mathcal{F})$. Set
$\mathcal{L}\leftarrow\emptyset$. For each pair $(i,j)\in V_{C}\times V_{C}$,
such that $i\neq j$ and $\nu_{ij}\geq 1$, set
$\mathcal{L}\leftarrow\mathcal{L}\cup\\{(i,j)\\}$.
If $\mathcal{L}\neq\emptyset$, do Step 6′c; else if $\delta>\sum_{i\in
V_{\text{disc}}}(\mu^{+}_{i}+\mu^{-}_{i})$, do Step 6a; else, do Step 6b:

    1. (c)

Sort $\mathcal{L}$ in decreasing order of $\nu$ (breaking ties in increasing
order of $\Delta$). Let $(i,j)$ be the first element of $\mathcal{L}$ and let
$(s_{1},s_{2})=\mathop{\arg\min}_{(s,s^{\prime})\in\mathcal{S}_{ij}\times\mathcal{S}_{ji}}\Big{\\{}d_{ij}^{s}+d_{ji}^{s^{\prime}}:d_{ij}^{s}+d_{ji}^{s^{\prime}}>w_{i}+w_{j}\Big{\\}}$.
Set $d_{1}^{\star}$ and $d_{2}^{\star}$ as follows: if $d_{ij}^{s_{1}}\leq
d_{ji}^{s_{2}}$, set $d_{1}^{\star}\leftarrow d_{ij}^{s_{1}}$,
$d_{2}^{\star}\leftarrow w_{i}+w_{j}-d_{1}^{\star}+\varepsilon$; otherwise,
set $d_{2}^{\star}\leftarrow d_{ji}^{s_{2}}$, $d_{1}^{\star}\leftarrow
w_{i}+w_{j}-d_{2}^{\star}+\varepsilon$; here, $\varepsilon$ is a small
positive number. Set $i^{\star}\leftarrow i$, $j^{\star}\leftarrow j$. Add
path elimination constraints: (i)
$\mathcal{F}^{\text{L}}\leftarrow\mathcal{F}^{\text{L}}\cup\left\\{(i^{\star},j^{\star},d_{1}^{\star})\right\\}$,
(ii)
$\mathcal{F}^{\text{R}}\leftarrow\mathcal{F}^{\text{R}}\cup\left\\{(j^{\star},i^{\star},d_{2}^{\star})\right\\}$.

Set
$\mathcal{N}\leftarrow\mathcal{N}\cup\left\\{\left(\tau^{\text{L}},\mathcal{F}^{\text{L}}\right),\left(\tau^{\text{R}},\mathcal{F}^{\text{R}}\right)\right\\}$,
and go to Step 2.

The modified algorithm gives preference to the branching Step 6′c over Steps
6a and 6b, because unlike the latter, the former would necessarily cut off the
VRPTW route set of the parent node in at least one scenario (in both children
nodes). However, as discussed earlier, the path-based disjunctions (8) (upon
which branching rule 6′c is based) are only necessary but not sufficient. This
is in contrast to the time window-based disjunctions (5) (upon which the
branching rules 6a and 6b are based), which are both necessary and sufficient.

Figure 7 shows the effect that the new branching rules have on the branch-and-
bound search tree for the case of our illustrative example from Figure 3. Note
how, in the root node, instead of branching via the time window-based
disjunctions (5), the modified branching Step 6′ certifies that it is
impossible to have both $3-4$ and $4-3$ paths in different scenario route
sets, and branches using the disjunction (8) instead. Observe that the search
tree is much smaller than in the case of Figure 4. Our numerical experiments
confirm that this is generally true, i.e., that incorporating the path-based
disjunctions (8) results in fewer nodes being explored (see Section 6.3).

0510152025$s=2$$s=1$12341234Branch as per (8):
$(i^{\star},j^{\star},d_{1}^{\star},d_{2}^{\star})=(4,3,7,5)$$\mathcal{F}=\emptyset$,
obj = 25.5 No $4-3$ paths with travel time $\geq 7$ No $3-4$ paths with travel
time $\geq 5$
0510152025$s=2$$s=1$12341234$\tau_{3}^{\star}$$\tau_{1}^{\star}$$\leftarrow\tau_{4}^{\star}$$\tau_{2}^{\star}$(optimal)$\mathcal{F}=\left\\{(4,3,7)\right\\}$,
obj =
26.00510152025$s=2$$s=1$12341234$\leftarrow\tau_{4}^{\star}$$\tau_{3}^{\star}$$\tau_{1}^{\star}$$\tau_{2}^{\star}$(also
optimal)$\mathcal{F}=\left\\{(3,4,5)\right\\}$, obj = 26.0 Figure 7: The
search tree of our algorithm (utilizing path-based disjunctions) for the
illustrative example of Figure 3.

###  4.3 Generating Upper Bounds

A small value of the global upper bound $UB$ can significantly speed up the
solution process by fathoming more nodes of the search tree using the
fathoming Step 4. The algorithm presented in Section 4.1 updates $UB$ in Step
5 only. However, it is possible to update $UB$ more frequently by generating
candidate (feasible) TWAVRP solutions using the route sets
$\\{\mathbf{R}_{s}\\}_{s\in\mathcal{S}}$ obtained in the node processing Step
3. The basic idea is to select some scenario and use its solution as a
“template” [18], assigning the time windows based on the arrival times in this
template solution. Specifically, for every scenario $s^{\star}\in\mathcal{S}$
in which a new route set $\mathbf{R}_{s^{\star}}$ is computed in Step 3 of the
algorithm, we can attempt to generate a feasible TWAVRP solution
$(\tau^{h},\\{\mathbf{R}_{s}^{h}\\}_{s\in\mathcal{S}})$ using the following
procedure.

  1. 1.

Initialize arrival times. Let $\mathbf{R}\leftarrow\mathbf{R}_{s^{\star}}$ and
$\theta\leftarrow\theta_{s^{\star}}$. Let
$\bm{a}\in\mathcal{X}(\mathbf{R},[e,\ell];\theta)$ be the vector of arrival
times with minimum cumulative waiting time. That is, we recursively define
$a_{R_{k,0}}=e_{0}$ and
$a_{R_{k,l+1}}=\max\\{e_{R_{k,l+1}},a_{R_{k,l}}+t_{R_{k,l},R_{k,l+1}}\\}$ for
all $l\in\\{1,\ldots,n_{k}-1\\}$ and $k\in\\{1,\ldots,m\\}$. For $i\in V_{C}$
such that $q_{i}=0$ (that is, $i\neq R_{k,l}$ for any $k,l$), we define
$a_{i}=e_{i}$.

  2. 2.

Assign time windows. Let $\tau^{h}$ be defined as follows.

$\tau_{i}^{h}=\begin{dcases}\left[x_{i},x_{i}+w_{i}\right],\text{where
}x_{i}=\begin{dcases}e_{i}&\text{if }a_{i}-w_{i}/2\leq e_{i}\\\
\ell_{i}-w_{i}&\text{if }a_{i}+w_{i}/2\geq\ell_{i}\\\
a_{i}-w_{i}/2&\text{otherwise}\end{dcases}&\text{if }i\in V_{\text{cont}}\\\
\left[x_{i,b_{i}},y_{i,b_{i}}\right],\text{where
}b_{i}\in\mathop{\arg\min}_{b\in\\{1,\ldots,N_{i}\\}}\left\\{\min_{\omega\in[x_{i,b},y_{i,b}]}\left\lvert
a_{i}-\omega\right\rvert\right\\}&\text{if }i\in
V_{\text{disc}}\end{dcases}\quad\forall i\in V_{C}.$

  3. 3.

Compute upper bound. For each $s\in\mathcal{S}$, let $\mathbf{R}^{h}_{s}$ be a
(possibly suboptimal) solution of $\mathtt{VRPTW}(\tau^{h};\theta_{s})$. Let
$ub^{h}\leftarrow\sum_{s\in\mathcal{S}}p_{s}c(\mathbf{R}^{h}_{s})$. If
$ub^{h}\leq UB$, set $UB\leftarrow ub^{h}$ and
$\tau^{\star}\leftarrow\tau^{h}$.

We remark that it is not necessary to solve the VRPTW instances to optimality
in Step 3 above, since the generated upper bounds $ub^{h}$ are guaranteed to
still be valid. Consequently, we can utilize any (possibly heuristic) VRPTW
solver to quickly compute candidate upper bounds. For example, the
computational results reported in Section 6 were obtained by implementing
Solomon’s sequential insertion construction heuristic “I1” [30] in combination
with a local search procedure [7], which used the Relocate, 2-opt, 2-opt⋆ and
Or-opt moves within a deterministic Variable Neighborhood Descent algorithm
(e.g., see [6]).

###  4.4 Modification as a Heuristic Algorithm

The exact algorithm of Section 4.1 can be readily modified as a heuristic
algorithm. Indeed, if one uses a heuristic VRPTW solver (e.g., one based on a
metaheuristic) in place of an exact one in the node processing Step 3, then
the time window assignment $\tau^{\star}$ determined by the algorithm is still
guaranteed to be feasible, although not necessarily optimal. Overall TWAVRP
optimality can no longer be guaranteed, since the fathoming Step 4 may
incorrectly prune a node whose descendant contains the optimal time window
assignment. Nevertheless, this modification as a heuristic is particularly
suited from a practical viewpoint, as it allows one to utilize any available
VRPTW heuristic solver “out of the box.”

##  5 Exact Solution of VRPTW Subproblems

Various exact solution schemes have been proposed for the VRPTW over the last
several decades. These include algorithms based on branch-and-cut, Lagrangean
relaxation and column generation, among others. We refer the reader to [11]
for a recent survey. The most successful of these are branch-price-and-cut
algorithms, which correspond to branch-and-bound algorithms in which the
bounds are obtained by solving linear relaxations of a _set partitioning_
model by column generation, and are further strengthened by generating cutting
planes.

In order to solve VRPTW instances in Steps 3 and 3′ of the main algorithm, we
implemented the branch-price-and-cut algorithm described in [27], as well as
procedures to warm start the latter exact method in the context of our
algorithm. Our implementation incorporates several elements of the algorithm
described in [27], including $ng$-routes, bidirectional labeling, variable
fixing, route enumeration, and limited-memory subset row cuts. In what
follows, we highlight only the most important of these ingredients; for
details, we refer the reader to [27].

###  5.1 Branch-Price-and-Cut Implementation

For ease of notation, we shall drop the subscript $s$ referencing a scenario
and assume that all operational parameters $\theta$ are fixed to certain given
values. We shall also assume that the time windows have been fixed at $\tau$
and that the set of forbidden paths is given to be $\mathcal{F}$. Note that we
only describe the solution approach for
$\mathtt{VRPTWFP}(\tau,\mathcal{F};\theta)$; the approach for
$\mathtt{VRPTW}(\tau;\theta)$ is obtained by simply setting
$\mathcal{F}=\emptyset$. We remark that, before we call the exact solution
method, we modify the VRPTW instance to obtain an equivalent one by tightening
the time windows $\tau$ and reducing the arc set $A$ using the preprocessing
routines described in [3, 21]. In addition to this, it is possible to further
reduce the arc set $A$ using members of $\mathcal{F}$. Indeed, if for some
$(i,j,d)\in\mathcal{F}$, the shortest travel time from $i$ to $j$ in graph $G$
exceeds $d$, then we can remove the arc $(i,j)$ from $A$.

In the following, $\mathcal{P}_{ij}$ denotes the set of all $i-j$ paths in
graph $G$ (after preprocessing), where $(i,j)\in V_{C}\times V_{C}$, such that
$i\neq j$. We use $\Omega$ to denote the set of all (elementary) vehicle
routes that are feasible with respect to capacity and time window constraints.
For a given route $r\in\Omega$, $\lambda_{ir}$ denotes the number of times
customer $i\in V_{C}$ is visited in route $r$, $\eta_{ijr}$ denotes the number
of times arc $(i,j)\in A$ is traversed by route $r$, while $c_{r}$ denotes its
cost, i.e., $c_{r}\equiv c(r)$. The set partitioning model is described in the
following. In this model, $x_{r}$ is a binary path-flow variable that encodes
whether route $r\in\Omega$ is part of the optimal route set.

$\begin{array}[]{r@{\;\;}l}\displaystyle\mathop{\text{minimize}}_{x}&\displaystyle\sum_{r\in\Omega}c_{r}x_{r}\\\
\text{subject to}&x_{r}\in\\{0,1\\},\quad\forall r\in\Omega,\\\
&\displaystyle\sum_{r\in\Omega}\lambda_{ir}x_{r}=1,\quad\forall i\in V_{C},\\\
&\displaystyle\sum_{r\in\Omega}\sum_{(i^{\prime},j^{\prime})\in\pi}\eta_{i^{\prime}j^{\prime}r}x_{r}\leq\lvert\pi\rvert-1,\quad\forall\pi\in\mathcal{P}_{ij}:t(\pi)\geq
d,\;\;\forall(i,j,d)\in\mathcal{F}.\end{array}$ (9)

In the above, the last set of inequalities are _infeasible path elimination
constraints_ (e.g., see [3, 21]) that forbid the occurrence of $i-j$ paths
with travel time greater than or equal to $d$. In our implementation, we
replace the subscript of the innermost summation with
$(i^{\prime},j^{\prime})\in\text{tr.cl.}(\pi)$, where $\text{tr.cl.}(\pi)$
denotes the transitive closure of $\pi$.333If $\pi=(v_{1},\ldots,v_{p})$
denotes an elementary path, then its transitive closure is the set of arcs
$(v_{k},v_{l})$ such that $v_{l}$ can be reached from $v_{k}$ using only arcs
in $\pi$, i.e, $\text{tr.cl.}(\pi)=\\{(v_{k},v_{l})\in
A:(k,l)\in\\{1,\ldots,p\\}\times\\{1,\ldots,p\\},k<l\\}$. This so-called
_tournament form_ of the inequality is stronger than the version presented
above, see [3] for a proof. Furthermore, it is also well known that one can
relax the feasible space of the above set partitioning model by including non-
elementary vehicle routes in $\Omega$ without sacrificing optimality. In our
implementation, we replace $\Omega$ with the set of so-called $ng$-routes
$\Omega^{\text{ng}}\supseteq\Omega$, which are not necessarily elementary [5].

We now describe the branch-price-and-cut algorithm to solve the set
partitioning model over $\Omega^{\text{ng}}$. The root node, which solves the
linear relaxation of this set partitioning model, is initialized with a subset
of $\Omega^{\text{ng}}$ (single-customer vehicle routes) but no infeasible
path elimination constraints. A pricing subproblem is used to generate other
members of $\Omega^{\text{ng}}$ (also referred to as columns), as necessary.
After column generation has converged, if the gap between the current node
lower bound and global upper bound is sufficiently small ($\leq 1\%$ in our
implementation), we employ route enumeration to generate all feasible vehicle
routes with reduced costs smaller than this gap [4]. [4] have shown that this
subset must contain the routes of all optimal solutions. Therefore, we can
solve the resulting “reduced” set partitioning model by using a standard
integer programming solver. We remark that this is done only if the number of
generated routes is less than a threshold, which we set to $3\times 10^{6}$,
as suggested in [27]. On the other hand, if the current node gap is large or
if route enumeration generated too many routes, then we attempt to separate
the infeasible path elimination constraints (as well as other valid
inequalities) in order to tighten the linear relaxation. The above process is
iterated until we cannot generate any more columns or inequalities. At this
stage, if the current node solution is fractional, we create additional
children nodes by branching on the number of used vehicles or by branching on
edges/arcs.

#### Pricing subproblem.

The pricing subproblem is a shortest path problem with resource constraints,
with customer demands and arc travel times considered as resources constrained
by the vehicle capacity and time windows, respectively. We utilize the dynamic
programming algorithm described in [27] to solve this problem. In order to
speed up the solution of the pricing subproblem, we apply various techniques
including bidirectional labeling and variable fixing (based on reduced costs).
In addition, we implemented the bucket pruning heuristic (e.g., see [15]) to
find candidate columns and use the dynamic programming algorithm only if the
former fails to generate columns.

#### Route enumeration.

We use the dynamic programming algorithm of [27] with a modification to
account for the presence of the infeasible path elimination constraints. In
particular, we consider two different “partial routes” that visit the same set
of customers (but possibly in different sequences) to be undominated
irrespectively of their resource consumptions, and we do not perform any
associated dominance checks. This prevents incorrectly “pruning” a vehicle
route that satisfies the infeasible path elimination constraints in favor of
one that does not. Among all enumerated routes, we only consider those which
are elementary and satisfy all infeasible path elimination constraints (which
can be done in $\mathcal{O}(\left\lvert\mathcal{F}\right\rvert)$ time per
route) to include in the final “reduced” set partitioning model.

#### Cutting planes.

We use the tournament form of the infeasible path elimination constraints as
“necessary cuts” (a violating member is guaranteed to be separated), and the
extended capacity cuts [28] and limited-node-memory subset row cuts [27] as
“strengthening cuts” (a violating member may not necessarily be separated).
The separation algorithms for these cuts are described next.

  * •

The infeasible path elimination constraints are separated by utilizing the
polynomial-time path-growing scheme of [3]. Specifically, suppose
$G_{\bar{x}}=(V,A_{\bar{x}})$ is the so-called support graph of the current
linear programming solution $\bar{x}$, where
$A_{\bar{x}}=\left\\{(i^{\prime},j^{\prime})\in
A:\sum_{r\in\Omega^{\text{ng}}}\eta_{i^{\prime}j^{\prime}r}\bar{x}_{r}>0\right\\}$.
Then, for every $(i,j,d)\in\mathcal{F}$, the scheme of [3] is used to obtain
the set of all $i-j$ paths $\pi$ in $G_{\bar{x}}$, which satisfy
$\sum_{r\in\Omega^{\text{ng}}}\sum_{(i^{\prime},j^{\prime})\in\text{tr.cl.}(\pi)}\eta_{i^{\prime}j^{\prime}r}\bar{x}_{r}>\lvert\pi\rvert-1$.
Amongst all such paths, we choose the ones for which the travel time $t(\pi)$
is greater than $d$ and add the corresponding tournament form of the
infeasible path elimination constraints to the current linear relaxation.

  * •

The extended capacity cuts are separated by first using the CVRPSEP package
[26] to separate the so-called rounded capacity inequalities and then lifting
these, as described in [28].

  * •

The limited-node-memory subset row cuts are separated by first identifying all
subset row cuts violated by node sets with cardinality up to $5$ and then
identifying their so-called “node memory sets,” as described in [27].

We remark that, since the infeasible path elimination inequalities are defined
over arcs, they are “robust” and affect the pricing subproblem only through a
corresponding term in the arc cost, i.e., their dual value. In other words,
the addition of these inequalities does not affect the complexity of the
pricing subproblem. In contrast, the extended capacity and limited-node-memory
subset row cuts are “non-robust” as their addition increases the complexity of
the pricing subproblem [27].

###  5.2 Warm Starting

Initializing the branch-price-and-cut algorithm described in the previous
section with a feasible set of columns can speed up the convergence of its
column generation process, leading to small computation times. In addition to
this, providing a valid initial upper bound on the optimal objective value can
also significantly speed up the search, both in the context of route
enumeration, where it can result in fewer routes being enumerated, as well as
branch-and-bound, where more parts of the search tree can be fathomed early in
the process. In this section, we describe warm starting procedures that can be
employed in the context of our algorithm of Section 4.

Consider a parent node $(\tau,\mathcal{F})$ (already processed) with optimal
solution $\\{\mathbf{R}_{s}\\}_{s\in\mathcal{S}}$ as obtained in Step 3. Also,
consider one of its child nodes $(\tau^{h},\mathcal{F}^{h})$, where
$h\in\\{\text{L},\text{R}\\}$, created in the branching Step 6′. Finally,
consider a scenario $s\in\mathcal{S}$ such that the corresponding optimal
route set in the parent node, $\mathbf{R}_{s}$, is not feasible for the child
node $(\tau^{h},\mathcal{F}^{h})$. This means that $\mathbf{R}_{s}$ cannot be
directly transferred to the latter, and warm starting is sought to aid the
search towards the optimal one.

#### Generating an initial set of columns.

Let $\Omega^{\text{ng}}_{s}$ be the set of columns that were generated during
the branch-price-and-cut process in scenario $s$ of the parent node
$(\tau,\mathcal{F})$. Since the child node $(\tau^{h},\mathcal{F}^{h})$
differs from its parent in exactly one constraint, it is likely that several
members of $\Omega^{\text{ng}}_{s}$ are also feasible in the set partitioning
model of the child’s scenario$-s$ VRPTW subproblem. Therefore, we can simply
loop through the members of $\Omega^{\text{ng}}_{s}$ and filter out all
infeasible columns to generate the initial linear relaxation in the branch-
price-and-cut algorithm.

#### Generating an initial upper bound.

The following procedure computes a valid upper bound, $ub_{s}^{h}$, on the
cost of the child’s scenario$-s$ VRPTW subproblem. The procedure attempts to
“repair” the scenario$-s$ optimal route set of the parent node and generate
one that is feasible for the child.

  1. 1.

Set $ub_{s}^{h}\leftarrow
UB-\sum_{s^{\prime}\in\mathcal{S}\setminus\\{s\\}}c(\mathbf{R}_{s^{\prime}})$,
where $UB$ is the currently applicable upper bound from the algorithm of
Section 4; and set $\mathbf{R}^{\prime}\leftarrow\mathbf{R}_{s}$.

  2. 2.

Let $H\subseteq V_{C}$ be defined as follows: if $(\tau^{h},\mathcal{F}^{h})$
was created using branching Steps 6a or 6b, then $H=\\{i^{\star}\\}$;
otherwise $H=\\{i^{\star},j^{\star}\\}$.

  3. 3.

For each $i\in H$:

    1. (a)

Remove customer $i$ from its current position in $\mathbf{R}^{\prime}$ and
insert it into a new vehicle route.

    2. (b)

Apply local search on $\mathbf{R}^{\prime}$, ensuring that each accepted move
satisfies all time window and path elimination constraints in
$(\tau^{h},\mathcal{F}^{h})$.

    3. (c)

If $c(\mathbf{R}^{\prime})<ub_{s}^{h}$, set $ub_{s}^{h}\leftarrow
c(\mathbf{R}^{\prime})$.

In our implementation of local search, we considered the Relocate, 2-opt,
2-opt⋆ and Or-opt moves within a deterministic Variable Neighborhood Descent
algorithm (e.g., see [6]). We remark that the solution of (up to $S$) VRPTW
instances in Step 3 of the algorithm can be easily parallelized.
Alternatively, one can do this serially on a single CPU thread (e.g., by
starting with the lowest indexed scenario). In the former setting, no
information can be exchanged among the VRPTW instances. In the latter setting,
however, one can capitalize on instances that have already been solved so as
to obtain improved upper bounds in Step 1 of the above procedure as follows:

$ub_{s}^{h}\leftarrow
UB-\sum_{s^{\prime}\in\mathcal{S}:s^{\prime}<s}c(\mathbf{R}_{s^{\prime}}^{h})-\sum_{s^{\prime}\in\mathcal{S}:s^{\prime}>s}c(\mathbf{R}_{s^{\prime}}),$

where $\mathbf{R}_{s^{\prime}}^{h}$ is the just computed optimal solution in
scenario $s^{\prime}$ of the node under consideration.

##  6 Computational Results

This section presents computational results obtained by our algorithm on
benchmark instances from the literature. Specifically, in Section 6.1, we
present the characteristics of the test instances; in Section 6.2, we present
a summary of the numerical performance of our algorithm and compare it to
existing solution methods; in Section 6.3, we present detailed tables of
results outlining the performance of each component of our algorithm; and,
finally in Section 6.4, we present results of a parallel implementation on
instances containing a large number of scenarios.

The algorithm was coded in C++ and the runs were conducted on an Intel Xeon
E5-2687W 3.1 GHz processor with 4 GB of available RAM. The nodes in Step 2 of
the algorithm were selected using a simple depth-first rule that backtracked
whenever the gap between the objective value of the current node and the
current upper bound exceeded 50% of the gap between the global lower and
current upper bound. All subordinate linear and mixed-integer linear programs
were solved using default settings of the IBM ILOG CPLEX Optimizer 12.7.
Finally, except for the results presented in Section 6.4, all runs were
restricted to a single CPU thread. This facilitates a fair comparison with
existing algorithms from the literature in Section 6.2 and between different
configurations of our algorithm in Section 6.3. The results presented in
Section 6.4 were obtained with OpenMP by using up to $\min\\{S,10\\}$ threads
in parallel, where $S$ is the number of scenarios. In all cases, an overall
“wall clock” time limit of one hour per instance was imposed.

###  6.1 Benchmark Instances

Existing benchmark instances for the TWAVRP focus solely on demand
uncertainty. For the continuous setting, the authors of [33] introduced 40
randomly generated instances. The number of customers ($n$) in these instances
varies from 10 to 25. Subsequently, [10] proposed 50 additional instances with
$n$ varying from 30 to 50. Each of the 90 available instances consists of 3
demand scenarios (low, medium, high), each with equal probability of
occurrence. The average demand (for each customer) across the three scenarios
is about $1/6$ of the vehicle capacity $Q$. The exogenous time windows are
designed to be much wider than the endogenous time windows; in particular, the
average (across customers) exogenous time window width ($\ell_{i}-e_{i}$) is
10.8, compared to an endogenous time window width ($w_{i}$) of just 2.0. For
the discrete setting, the authors of [32] introduced 80 randomly generated
instances with $n$ varying from 10 to 60. Except for the structure of the
feasible time window set, the instances share similar characteristics as in
the continuous setting, each consisting of 3 demand scenarios. In each
instance, the number of candidate time windows ($N_{i}$) is equal to 3 for
about 10% of the customers, 5 for about 60% of the customers, and 7 for the
remaining 30%. Similarly to the continuous setting, the customer locations are
generated using a uniform distribution over a square with the depot located in
the center. Moreover, the time windows and vehicle capacities are chosen such
that no more than eight customers can be visited on a single vehicle route in
any scenario. All of the aforementioned test instances are inspired from the
Dutch retail sector, and can be found online at
http://people.few.eur.nl/spliet.

To test our algorithm on instances containing a large number of scenarios, we
generated 80 additional benchmark instances. Specifically, for each existing
(continuous and discrete) TWAVRP instance with $n\leq 25$, we used a similar
procedure as described in [32] to generate 15 additional demand scenarios. For
a particular instance, we first generate a nominal demand $\bar{q}_{i}$ for
each customer $i\in V_{C}$ using a normal distribution with mean $5.0$ and
variance $1.5$. To generate additional scenarios, we draw additive
disturbances $\epsilon_{si}$ from a uniform distribution on $[-1.5,1.5]$ for
each $i\in V_{C}$ and $s\in\mathcal{S}$, and multiplicative factors $f_{s}$
from a uniform distribution on $[0.625,1.375]$ for each $s\in\mathcal{S}$. The
demand of customer $i$ in scenario $s$ is then computed as
$q_{si}=\left\lceil\max\left\\{f_{s}\left(\bar{q}_{i}+\epsilon_{si}\right),10^{-6}\right\\}\right\rceil$.
The multiplicative factors $f_{s}$ determine the level of correlation among
the customer demands. For instance, high (low) values of $f_{s}$ may represent
the behavior that demands increase (decrease) uniformly for all retailers in a
supply chain, whereas values close to one represent the nominal situation in
which the demands are uncorrelated. The additional benchmark instances can be
downloaded from http://gounaris.cheme.cmu.edu/datasets/twavrp.

###  6.2 Comparison with Existing Methods

We first compare the performance of our algorithm with the results published
in [10] for the case of the continuous TWAVRP. We do not compare with the
algorithm of [33], since the authors of [10] have demonstrated that their
algorithm is superior to the former. Table 1 summarizes the comparison of the
numerical performance across all 90 instances that are available for the
continuous TWAVRP. The column # denotes the number of test instances that
contain $n$ customers. For each algorithm, Optimal denotes the number of test
instances that it could solve to optimality in one hour while Time (sec)
denotes the average time in seconds to solve these instances to optimality.
For those instances which could not be solved to optimality in one hour, the
column Gap (%) reports the average optimality gap, defined as $(UB-
LB)/UB\times 100\%$, where $LB$ and $UB$ are respectively the global lower and
upper bounds determined by the algorithm after one hour. The two methods are
also compared in the performance profiles [12] of Figure 8(a). Our proposed
algorithm is able to solve all but one (89 out of 90) benchmark instances to
optimality, utilizing an average computation time of 169 seconds; of these, 32
instances were unsolved by the best previous method, while the one unsolved
instance was determined to be within 0.8% of optimality. These results
demonstrate that our algorithm strongly outperforms the existing method,
solving more problems and achieving (or matching) the fastest computation time
in all instances.

Table 1: Computational comparison of the proposed algorithm against the existing state-of-the-art algorithm (DS18) [10] on all 90 benchmark instances of the continuous TWAVRP. |  | DS18 | Proposed algorithm
---|---|---|---
$n$ | # |  Optimal |  Time (sec) |  Gap (%) |  Optimal |  Time (sec) |  Gap (%)
10 | 10 |  10 |  0.1 |  – |  10 |  0.1 |  –
15 | 10 |  10 |  4.6 |  – |  10 |  0.6 |  –
20 | 10 |  10 |  2.2 |  – |  10 |  1.5 |  –
25 | 10 |  10 |  12.4 |  – |  10 |  8.6 |  –
30 | 10 |  9 |  204.5 |  1.66 |  10 |  48.2 |  –
35 | 10 |  6 |  152.8 |  0.89 |  9 |  51.5 |  0.79
40 | 10 |  2 |  1860.0 |  1.16 |  10 |  342.3 |  –
45 | 10 |  0 |  – |  2.74 |  10 |  361.3 |  –
50 | 10 |  0 |  – |  4.26 |  10 |  696.0 |  –
All | 90 |  57 |  117.0 |  2.56 |  89 |  169.1 |  0.79
Processor |  | Intel i7 3.5 GHz | Xeon E5-2687W 3.1GHz

We now turn our attention to the discrete TWAVRP. Table 2 compares the
numerical performance of our algorithm with the results published in [32]
across all 80 instances of the discrete TWAVRP. The columns in this table have
the same meaning as in Table 1. The two algorithms are also compared in the
performance profiles of Figure 8(b). Our algorithm is able to solve 54 out of
80 benchmark instance to optimality, utilizing an average computation time of
274 seconds; of these, 22 instances were unsolved by the best previous method,
while the remaining unsolved instances were determined to be within 1.2% of
optimality, on average. As in the continuous setting, our algorithm strongly
outperforms the existing method: it solves more instances and achieves the
fastest computation time in all of them.

Table 2: Computational comparison of the proposed algorithm against the
existing state-of-the-art algorithm (SD15) [32] on all 80 benchmark instances
of the discrete TWAVRP.

|  | SD15 | Proposed algorithm
---|---|---|---
$n$ | # |  Optimal |  Time (sec) |  Gap (%) |  Optimal |  Time (sec) |  Gap (%)
10 | 10 |  10 |  3.9 |  – |  10 |  0.1 |  –
15 | 10 |  10 |  185.9 |  – |  10 |  15.2 |  –
20 | 10 |  9 |  1247.6 |  0.06 |  10 |  33.8 |  –
25 | 10 |  3 |  504.4 |  n/a† |  9 |  248.8 |  0.43
30 | 10 |  0 |  – |  n/a† |  9 |  581.7 |  0.13
40 | 10 |  0 |  – |  n/a† |  5 |  1263.5 |  1.22
50 | 10 |  0 |  – |  n/a† |  1 |  533.6 |  1.31
60 | 10 |  0 |  – |  n/a† |  0 |  – |  1.30
All | 80 |  32 |  457.5 |  n/a† |  54 |  274.4 |  1.21
Processor |  | Intel Core i5-2450M 2.5 GHz | Xeon E5-2687W 3.1GHz

  * •

Note. The reported results for “SD15” are the best entries of Tables 3 and 4
from that publication [32].

  * •

† The optimality gaps for the unsolved instances have not been reported in the
publication.

(a) Continuous

(b) Discrete

Figure 8: Log-scaled performance profiles across all benchmark instances. The
left graph compares the performance in the continuous setting, in which “DS18”
shows the performance of the algorithm of [10], while the right graph compares
the performance in the discrete setting, in which “SD15” shows the performance
of the algorithm of [32]. In both graphs, “This work” shows the performance of
our proposed algorithm. For each curve (i.e., algorithm), the value at $t=0$
gives the fraction of benchmark instances for which it is fastest, while the
limiting value at $t\to\infty$ gives the fraction of instances which it could
solve within the time limit of one hour.

###  6.3 Detailed Discussion of Results

A comparison of Tables 1 and 2 shows that the discrete TWAVRP instances take
longer to solve than the continuous ones. This can be partly explained by the
fact that, in the continuous setting, the separation problem (7) is a linear
program, while in the discrete setting, it is a mixed-integer linear program.
Consequently, the algorithm spends a greater fraction of the total time in
solving the separation problem in the latter case (see Table 3).

Table 3: Percentage of computing time spent in various parts of the algorithm (averaged across instances solved to optimality in one hour). “Solving VRPTW” and “Separation problem” refer to Steps 3′ and 5 of the algorithm from Section 4.2 respectively, while “Upper bounding” refers to the steps in Section 4.3. | Continuous | Discrete
---|---|---
$n$ |  Solving VRPTW |  Separation problem |  Upper bounding |  Solving VRPTW |  Separation problem |  Upper bounding
25 |  84.6 |  0.0 |  15.0 |  81.0 |  15.4 |  3.3
30 |  86.8 |  0.1 |  12.9 |  75.7 |  21.4 |  2.7
40 |  89.2 |  0.2 |  10.4 |  90.6 |  8.2 |  1.2
50 |  93.1 |  0.0 |  6.8 |  95.7 |  1.6 |  2.8
$\geq 25$ |  89.0 |  0.1 |  10.8 |  81.6 |  15.6 |  2.6

To show the efficacy of the path-based disjunctions and the associated
branching rules (see Section 4.2), we disable them and run only the basic
version of the algorithm from Section 4.1. Table 4 compares the performance of
this basic version with the one incorporating the path-based disjunctions.
They are also compared in the performance profiles of Figure 9. The results
indicate that the path-based branching rules are important to improve the
tractability of the overall algorithm. In particular, they are essential in
reducing the total number of nodes that are processed (i.e., the number of
times Steps 3 and 3′ are executed in the overall search trees) from $2,400$ to
$100$, on average. We remark, however, that this reduction comes at a price:
it requires modifying the underlying VRPTW solver (see Section 5.1).
Nevertheless, even without the path-based branching rules, the basic version
of our algorithm outperforms the existing ones (see Figure 9), while having
the advantage of being able to utilize any VRPTW solver in a modular fashion.

Table 4: Computational comparison of the algorithm with and without the path-based disjunctions on all 170 benchmark instances of the continuous and discrete TWAVRP. |  | Without path-based disjunctions | With path-based disjunctions
---|---|---|---
$n$ | # |  Optimal |  Nodes |  Time (sec) |  Gap (%) |  Optimal |  Nodes |  Time (sec) |  Gap (%)
$\left[10,15\right]$ | 40 |  40 |  5,887 |  91.5 |  – |  40 |  26 |  4.0 |  –
$\left[20,25\right]$ | 40 |  36 |  933 |  147.7 |  0.17 |  39 |  158 |  68.7 |  0.43
$\left[30,35\right]$ | 30 |  22 |  949 |  327.6 |  0.16 |  28 |  132 |  220.7 |  0.46
$\left[40,45\right]$ | 30 |  19 |  469 |  210.7 |  0.66 |  25 |  139 |  534.1 |  1.22
$\left[50,60\right]$ | 30 |  8 |  557 |  1,061.0 |  1.12 |  11 |  105 |  681.2 |  1.30
All | 170 |  125 |  2,427 |  229.4 |  0.75 |  143 |  108 |  208.9 |  1.19

(a) Continuous

(b) Discrete

Figure 9: Log-scaled performance profiles across all benchmark instances. The
profile “This work (no path)” refers to our proposed algorithm without the
path-based disjunctions (see Section 4.2). The other profiles as well as the
axes have the same meaning as in Figure 8.

Tables 5 and 6 present detailed results on all benchmark instances of the
continuous and discrete TWAVRP, respectively. In these tables, if an instance
could be solved to optimality within one hour, then Opt [UB] reports the
corresponding optimal objective value, while Time (sec) [LB] reports the time
to solve the instance to optimality. Otherwise, the columns respectively
report in brackets the best upper and lower bounds found within the time limit
of one hour.

Table 5: Detailed performance of the proposed algorithm on all 90 benchmark instances of the continuous TWAVRP. Instance | $n$ |  Opt [UB] |  Time (sec) [LB] | Instance | $n$ |  Opt [UB] |  Time (sec) [LB] | Instance | $n$ |  Opt [UB] |  Time (sec) [LB]
---|---|---|---|---|---|---|---|---|---|---|---
1 | 10 |  17.65 |  0.0 | 31 | 25 |  31.43 |  2.6 | 61 | 40 |  46.13* |  9.7
2 | 10 |  15.56 |  0.2 | 32 | 25 |  30.71 |  3.0 | 62 | 40 |  48.35 |  27.8
3 | 10 |  17.42 |  0.0 | 33 | 25 |  33.71 |  9.9 | 63 | 40 |  44.48* |  18.2
4 | 10 |  18.51 |  0.2 | 34 | 25 |  33.34 |  3.4 | 64 | 40 |  43.75 |  192.2
5 | 10 |  16.07 |  0.1 | 35 | 25 |  29.05 |  3.1 | 65 | 40 |  43.39* |  20.2
6 | 10 |  18.00 |  0.0 | 36 | 25 |  30.50 |  3.9 | 66 | 40 |  44.68* |  62.7
7 | 10 |  17.02 |  0.0 | 37 | 25 |  28.68 |  17.2 | 67 | 40 |  46.88* |  110.2
8 | 10 |  23.89 |  0.1 | 38 | 25 |  35.69 |  38.5 | 68 | 40 |  44.96* |  2,506.0
9 | 10 |  20.31 |  0.1 | 39 | 25 |  32.55 |  2.5 | 69 | 40 |  43.07* |  31.1
10 | 10 |  16.31 |  0.1 | 40 | 25 |  32.14 |  2.0 | 70 | 40 |  43.00* |  445.3
11 | 15 |  17.78 |  0.6 | 41 | 30 |  36.38 |  2.7 | 71 | 45 |  50.65* |  20.4
12 | 15 |  27.10 |  1.3 | 42 | 30 |  34.69* |  9.7 | 72 | 45 |  51.74* |  69.4
13 | 15 |  29.37 |  0.4 | 43 | 30 |  35.48 |  285.4 | 73 | 45 |  41.70* |  39.6
14 | 15 |  23.18 |  1.9 | 44 | 30 |  35.88 |  19.9 | 74 | 45 |  47.77* |  228.6
15 | 15 |  24.15 |  0.3 | 45 | 30 |  35.55 |  26.1 | 75 | 45 |  49.39* |  582.0
16 | 15 |  21.03 |  0.3 | 46 | 30 |  37.47 |  11.8 | 76 | 45 |  49.83* |  1,082.5
17 | 15 |  22.04 |  0.3 | 47 | 30 |  32.54 |  5.7 | 77 | 45 |  51.09* |  1,241.8
18 | 15 |  22.30 |  0.4 | 48 | 30 |  36.32 |  7.0 | 78 | 45 |  53.33* |  102.2
19 | 15 |  26.52 |  0.4 | 49 | 30 |  35.30 |  67.6 | 79 | 45 |  48.09* |  99.4
20 | 15 |  22.11 |  0.3 | 50 | 30 |  40.27 |  46.2 | 80 | 45 |  50.26* |  146.8
21 | 20 |  28.08 |  0.8 | 51 | 35 |  43.46 |  102.9 | 81 | 50 |  58.11* |  559.8
22 | 20 |  29.80 |  0.7 | 52 | 35 |  41.84 |  25.5 | 82 | 50 |  52.61* |  211.0
23 | 20 |  30.30 |  1.0 | 53 | 35 |  45.03* |  39.3 | 83 | 50 |  58.58* |  2,826.9
24 | 20 |  24.16 |  1.3 | 54 | 35 |  41.54* |  42.8 | 84 | 50 |  53.92* |  115.7
25 | 20 |  29.84 |  7.8 | 55 | 35 |  37.92 |  12.7 | 85 | 50 |  54.96* |  1,113.0
26 | 20 |  29.72 |  0.8 | 56 | 35 |  44.49* |  17.9 | 86 | 50 |  52.83* |  306.1
27 | 20 |  26.48 |  1.0 | 57 | 35 |  [41.04] |  [40.72] | 87 | 50 |  53.71* |  93.3
28 | 20 |  26.14 |  0.8 | 58 | 35 |  41.22 |  64.1 | 88 | 50 |  56.12* |  203.3
29 | 20 |  26.61 |  0.6 | 59 | 35 |  43.43 |  14.8 | 89 | 50 |  60.23* |  1,299.0
30 | 20 |  26.36 |  0.6 | 60 | 35 |  42.27 |  143.1 | 90 | 50 |  58.93* |  231.6
*Instances solved to optimality for the first time are indicated with an asterisk.

Table 6: Detailed performance of the proposed algorithm on all 80 benchmark instances of the discrete TWAVRP. Instance | $n$ |  Opt [UB] |  Time (sec) [LB] | Instance | $n$ |  Opt [UB] |  Time (sec) [LB] | Instance | $n$ |  Opt [UB] |  Time (sec) [LB]
---|---|---|---|---|---|---|---|---|---|---|---
1 | 10 |  12.83 |  0.1 | 31 | 25 |  35.47 |  24.3 | 61 | 50 |  [52.23] |  [51.57]
2 | 10 |  16.84 |  0.2 | 32 | 25 |  32.66* |  16.0 | 62 | 50 |  [55.61] |  [55.19]
3 | 10 |  16.60 |  0.1 | 33 | 25 |  [31.75] |  [31.61] | 63 | 50 |  [50.75] |  [50.09]
4 | 10 |  15.96 |  0.2 | 34 | 25 |  34.14* |  236.9 | 64 | 50 |  51.17* |  533.6
5 | 10 |  19.65 |  0.2 | 35 | 25 |  30.29* |  617.7 | 65 | 50 |  [54.11] |  [53.44]
6 | 10 |  18.13 |  0.3 | 36 | 25 |  32.54* |  611.6 | 66 | 50 |  [57.52] |  [56.57]
7 | 10 |  12.17 |  0.1 | 37 | 25 |  27.48* |  568.5 | 67 | 50 |  [58.14] |  [57.54]
8 | 10 |  17.09 |  0.2 | 38 | 25 |  34.83 |  15.9 | 68 | 50 |  [56.37] |  [55.27]
9 | 10 |  20.14 |  0.1 | 39 | 25 |  34.39* |  123.7 | 69 | 50 |  [53.85] |  [53.40]
10 | 10 |  17.17 |  0.1 | 40 | 25 |  30.73 |  25.0 | 70 | 50 |  [57.37] |  [56.37]
11 | 15 |  23.04 |  7.4 | 41 | 30 |  36.39* |  37.9 | 71 | 60 |  [64.83] |  [63.60]
12 | 15 |  25.27 |  1.0 | 42 | 30 |  40.59* |  245.4 | 72 | 60 |  [62.60] |  [61.57]
13 | 15 |  22.12 |  2.8 | 43 | 30 |  37.18* |  88.7 | 73 | 60 |  [64.92] |  [63.91]
14 | 15 |  18.46 |  0.7 | 44 | 30 |  [38.02] |  [37.97] | 74 | 60 |  [69.14] |  [68.59]
15 | 15 |  24.87 |  129.9 | 45 | 30 |  36.72* |  311.2 | 75 | 60 |  [63.61] |  [63.15]
16 | 15 |  19.82 |  2.4 | 46 | 30 |  34.76* |  237.8 | 76 | 60 |  [64.49] |  [63.77]
17 | 15 |  21.96 |  4.7 | 47 | 30 |  42.24* |  133.6 | 77 | 60 |  [61.24] |  [60.82]
18 | 15 |  22.93 |  0.7 | 48 | 30 |  37.04* |  3,501.3 | 78 | 60 |  [64.77] |  [63.06]
19 | 15 |  23.14 |  1.3 | 49 | 30 |  40.47* |  202.6 | 79 | 60 |  [65.48] |  [64.91]
20 | 15 |  18.84 |  0.7 | 50 | 30 |  39.89* |  477.2 | 80 | 60 |  [64.42] |  [63.76]
21 | 20 |  27.99 |  1.2 | 51 | 40 |  [41.96] |  [41.44] |  |  |  |
22 | 20 |  25.63 |  16.8 | 52 | 40 |  [47.43] |  [47.37] |  |  |  |
23 | 20 |  26.53 |  199.2 | 53 | 40 |  41.76* |  1,829.8 |  |  |  |
24 | 20 |  32.36 |  3.9 | 54 | 40 |  45.96* |  1,385.2 |  |  |  |
25 | 20 |  28.84 |  4.7 | 55 | 40 |  [48.68] |  [48.10] |  |  |  |
26 | 20 |  26.99 |  6.3 | 56 | 40 |  [44.88] |  [43.99] |  |  |  |
27 | 20 |  27.55* |  80.0 | 57 | 40 |  43.90* |  619.3 |  |  |  |
28 | 20 |  26.53 |  14.4 | 58 | 40 |  43.09* |  918.1 |  |  |  |
29 | 20 |  29.49 |  8.6 | 59 | 40 |  [49.27] |  [48.50] |  |  |  |
30 | 20 |  23.55 |  2.5 | 60 | 40 |  47.13* |  1,565.2 |  |  |  |
*Instances solved to optimality for the first time are highlighted with an asterisk.

###  6.4 Instances containing a Large Number of Scenarios

We now turn our attention to benchmark instances containing a large number of
scenarios. Our goals are two-fold. First, we aim to understand how our
algorithm performs as the number of considered scenarios ($S$) increases.
Second, we aim to understand the cost benefits of considering more scenarios
during strategic time window assignment. In pursuit of these goals, we
consider the 80 benchmark instances consisting of 15 demand scenarios each
(see Section 6.1). For each of these instances, we obtain time window
assignments using our algorithm by considering the following sample average
approximations: (i) the original instance with all 15 scenarios, and (ii) the
original instance with only the first $S$ scenarios, where
$S\in\\{1,3,5,10\\}$. For each approximation, we implement a parallel version
of our algorithm in which Step 3′ and the upper bounding step in Section 4.3
are each parallelized using up to $\min\\{S,10\\}$ threads.

For $S=1$ and $S=3$, all 80 instances were solved to optimality, while 73 out
of 80 instances were solved to optimality for $S=5$ (the average gap for the 7
unsolved instances being less than 0.2%). Therefore, at the interest of
brevity, we shall not show tabulated results for these cases. On the other
hand, for the cases with the higher number of scenarios, namely $S=10$ and
$S=15$, Table 7 presents a summary of the performance. The columns in this
table have the same interpretation as in Tables 1 and 2; the only difference
is that the column Time (sec) is now broken into two parts: Wall denotes the
average wall clock time, while CPU denotes the average CPU time. The ratio
between these quantities is a good measure of how well the algorithm scales
across multiple threads (i.e., how much it benefits from parallelism), which
is also plotted as a function of $S$ in Figure 10.

Table 7: Summary of computational performance of the parallelized algorithm on 80 benchmark instances of the (continuous and discrete) TWAVRP, each containing $S$ scenarios. |  | $S=10$ | $S=15$
---|---|---|---
|  |  | Time (sec) |  |  | Time (sec) |
$n$ | # |  Optimal |  Wall |  CPU |  Gap % |  Optimal |  Wall |  CPU |  Gap %
10 | 20 |  20 |  0.7 |  3.8 |  – |  20 |  1.1 |  7.0 |  –
15 | 20 |  18 |  127.7 |  474.6 |  0.21 |  18 |  441.7 |  1,639.2 |  0.64
20 | 20 |  16 |  266.4 |  779.8 |  0.46 |  11 |  426.2 |  1,403.3 |  0.33
25 | 20 |  5 |  465.3 |  841.6 |  0.60 |  2 |  2,189.9 |  3,777.1 |  0.70
All | 80 |  59 |  150.9 |  428.9 |  0.53 |  51 |  334.2 |  1,032.1 |  0.58
Figure 10: The ratio of CPU time to wall clock time with increasing number of
scenarios (for instances whose solution took at least one wall clock second
and at most one wall clock hour). For each $S$, the lower and upper most
dashes denote the minimum and maximum values of the ratio, the lower and upper
edges of the box denote the first and third quartile, while the middle bar (in
red color) denotes the median ratio.

Table 7 shows that we can consistently solve all benchmark instances with 15
scenarios to an optimality gap of less than 1% within a time limit of one
hour. Figure 10 shows that the speedup in wall clock time is sublinear with
respect to the number of parallel threads. The reasons for deviating from a
perfect linear speedup are two-fold. First, each node of our search tree does
not necessarily require the solution of $S$ VRPTW instances (see concluding
paragraph of Section 4.1). Indeed, for the considered benchmark instances, the
average number of VRPTW instances solved in a typical node is smaller than
$S/2$. Second, the variance in solution times across the VRPTW instances
solved in a node is typically large because the feasible route sets in a
particular scenario might be drastically different compared to other scenarios
in the same node (because of different demand realizations). Consequently, the
different CPU threads are not necessarily balanced, i.e., do not perform equal
amount of work. Nevertheless, as Figure 10 shows, on average, our algorithm
performs four times as many computations for instances containing 15 scenarios
by utilizing up to ten threads as compared to using just one.

Finally, Table 8 shows the cost savings in routing that are to be expected by
considering more than just one scenario during time window assignment. We
calculate these savings as follows. First, for a particular instance (with a
given $S$), we obtain the optimal time window assignment $\tau^{\star}$ from
our algorithm. Next, we calculate the “in-sample” expected costs by (i) fixing
the time windows to $\tau^{\star}$, (ii) solving a VRPTW instance with time
windows fixed to $\tau^{\star}$ for each of the 15 postulated demand scenarios
in the original benchmark instance, and (iii) averaging the costs. The “out-
of-sample” expected costs are obtained in exactly the same manner except that
in step (ii), the costs are evaluated over 100 independently generated demand
scenarios, which are randomly drawn using the procedure described in Section
6.1. Table 8 shows that, on average, we expect to do better than the
deterministic solution by about $2.3\%$ when considering up to $3$ scenarios
and by about $3.2\%$ when considering up to $15$ scenarios (based on out-of-
sample evaluations). If we consider only the discrete instances, the out-of-
sample savings increase to about $3.0\%$ for $S=3$ and $3.7\%$ for $S=15$. In
either case, we observe that the marginal benefits diminish as the number of
scenarios grows. Although we expect this trend to generally hold, we also
expect the actual magnitudes of the cost savings to depend on the problem
parameters as well as the dimensionality of the uncertainty (the number of
customers in this case). Indeed, for high-dimensional problems with narrow and
discrete time windows, we expect greater cost benefits by considering more
scenarios of the uncertainty.

Table 8: Expected cost savings from considering $S$ scenarios during strategic time window assignment relative to considering only one scenario (averaged across all 60 instances with $n\geq 15$). |  $S=1$ |  $S=3$ |  $S=5$ |  $S=10$ |  $S=15$
---|---|---|---|---|---
In-sample |  0.00% |  2.26% |  2.88% |  3.13% |  3.22%
Out-of-sample |  0.00% |  2.36% |  2.89% |  3.18% |  3.21%

##  7 Conclusions

While there has been significant research done on efficient algorithms for
vehicle routing problems, rapidly changing business models and advances in
technology are creating new challenges that are yet unsolved. This paper
attempts to address the challenge of dealing with operational uncertainty
during strategic decision-making in the context of vehicle routing. In
particular, we studied problems in which the decisions correspond to an
allocation of long-term delivery time windows to customers. These problems are
motivated from several real-world distribution operations, and are
particularly common in retail. We proposed a novel algorithm to solve this
problem that is highly competitive with existing methods. It draws on existing
algorithms for deterministic vehicle routing problems (both exact and
heuristic) and can use any vehicle routing solver as a subroutine, thus
facilitating its deployment in practice. From a modeling viewpoint, it allows
the user to postulate potential scenarios of future uncertainty corresponding
to different routing-specific parameters, as well as incorporate modular
changes in routing-specific constraints. Our business insights, aided via
numerical experiments, are that long-term costs are expected to decrease by
modeling more scenarios of future uncertainty but the marginal benefits
rapidly diminish as a function of the number of postulated scenarios. In other
words, a few scenarios are sufficient to obtain time window assignments that
would incur lower long-term costs compared to assignments that would be
obtained by completely ignoring operational uncertainty.

Our results open multiple avenues for future research. First, the scope of the
proposed algorithm can be broadened further by incorporating features such as
time-dependent travel times (which capture predictable variations in traffic
conditions) as well as soft time windows (which allow violation of the
assigned time windows in a controlled manner). Second, the time window
assignment model can be generalized to allow customers to express preferences
in the time window allocation process. More ambitiously, it would be
instructive to explore how we can reinterpret these concepts in the context of
consistent vehicle routing problems.

## Acknowledgments

We acknowledge support from the United States National Science Foundation,
award number CMMI-1434682. Anirudh Subramanyam and Akang Wang also acknowledge
support from, respectively, the John and Claire Bertucci Graduate Fellowship
program and the James C. Meade Fellowship program, at Carnegie Mellon
University.

## References

  * [1] Niels Agatz, Ann Campbell, Moritz Fleischmann, and Martin Savelsbergh.  Time Slot Management in Attended Home Delivery.  Transportation Science, 45(3):435–449, 2011.
  * [2] Niels Agatz, Ann M. Campbell, Moritz Fleischmann, and Martin Savelsbergh.  Challenges and opportunities in attended home delivery.  In Bruce Golden, S. Raghavan, and Edward Wasil, editors, The Vehicle Routing Problem: Latest Advances and New Challenges, pages 379–396. Springer US, Boston, MA, 2008.
  * [3] Norbert Ascheuer, Matteo Fischetti, and Martin Grötschel.  Solving the Asymmetric Travelling Salesman Problem with time windows by branch-and-cut.  Mathematical Programming, 90(3):475–506, 2001.
  * [4] Roberto Baldacci, Nicos Christofides, and Aristide Mingozzi.  An exact algorithm for the vehicle routing problem based on the set partitioning formulation with additional cuts.  Mathematical Programming, 115(2):351–385, 2008.
  * [5] Roberto Baldacci, Aristide Mingozzi, and Roberto Roberti.  New route relaxation and pricing strategies for the vehicle routing problem.  Operations Research, 59(5):1269–1283, 2011.
  * [6] Olli Bräysy.  A Reactive Variable Neighborhood Search for the Vehicle-Routing Problem with Time Windows.  INFORMS Journal on Computing, 15(4):347–368, 2003.
  * [7] Olli Bräysy and Michel Gendreau.  Vehicle routing problem with time windows, Part I: Route construction and local search algorithms.  Transportation Science, 39(1):104–118, 2005.
  * [8] Bruno P. Bruck, Jean-François Cordeau, and Manuel Iori.  A practical time slot management and routing problem for attended home services.  Omega, 81:208 – 219, 2018.
  * [9] Carlos F. Daganzo.  The length of tours in zones of different shapes.  Transportation Research Part B: Methodological, 18(2):135–145, 1984\.
  * [10] Kevin Dalmeijer and Remy Spliet.  A branch-and-cut algorithm for the Time Window Assignment Vehicle Routing Problem.  Computers & Operations Research, 89:140–152, 2018.
  * [11] Guy Desaulniers, Oli B.G. Madsen, and Stefan Ropke.  Chapter 5: The Vehicle Routing Problem with Time Windows.  In Paolo Toth and Daniele Vigo, editors, Vehicle Routing: Problems, Methods and Applications, pages 119–159. Society for Industrial and Applied Mathematics, Philadelphia, PA, 2014.
  * [12] Elizabeth D. Dolan and Jorge J. Moré.  Benchmarking optimization software with performance profiles.  Mathematical Programming, 91(2):201–213, 2002.
  * [13] Moshe Dror, Gilbert Laporte, and Pierre Trudeau.  Vehicle routing with stochastic demands: Properties and solution frameworks.  Transportation science, 23(3):166–176, 1989.
  * [14] Miguel Andres Figliozzi.  Planning approximations to the average length of vehicle routing problems with time window constraints.  Transportation Research Part B: Methodological, 43(4):438–447, 2009\.
  * [15] Ricardo Fukasawa, Humberto Longo, Jens Lysgaard, Marcus Poggi de Aragão, Marcelo Reis, Eduardo Uchoa, and Renato F. Werneck.  Robust branch-and-cut-and-price for the capacitated vehicle routing problem.  Mathematical Programming, 106(3):491–511, 2006.
  * [16] Charles Gauvin, Guy Desaulniers, and Michel Gendreau.  A branch-cut-and-price algorithm for the vehicle routing problem with stochastic demands.  Computers & Operations Research, 50:141–153, 2014.
  * [17] Michel Gendreau, Ola Jabali, and Walter Rei.  Chapter 8: Stochastic Vehicle Routing Problems.  In Paolo Toth and Daniele Vigo, editors, Vehicle Routing: Problems, Methods and Applications, volume 18, pages 213–239. Society for Industrial and Applied Mathematics, Philadelphia, PA, 2014.
  * [18] Chris Groër, Bruce Golden, and Edward Wasil.  The Consistent Vehicle Routing Problem.  Manufacturing & Service Operations Management, 11(4):630–643, 2009.
  * [19] Florent Hernandez, Michel Gendreau, and Jean-Yves Potvin.  Heuristics for tactical time slot management: a periodic vehicle routing problem view.  International Transactions in Operational Research, 24(6):1233–1252, 2017.
  * [20] Ola Jabali, Roel Leus, Tom van Woensel, and Ton de Kok.  Self-imposed time windows in vehicle routing problems.  OR Spectrum, 37(2):331–352, 2015.
  * [21] Brian Kallehauge, Natashia Boland, and Oli B.G. Madsen.  Path inequalities for the vehicle routing problem with time windows.  Networks, 49(4):273–293, 2007.
  * [22] Attila A. Kovacs, Bruce L. Golden, Richard F. Hartl, and Sophie N. Parragh.  Vehicle routing problems in which consistency considerations are important: A survey.  Networks, 64(3):192–213, 2014.
  * [23] Attila A. Kovacs, Sophie N. Parragh, and Richard F. Hartl.  A template-based adaptive large neighborhood search for the consistent vehicle routing problem.  Networks, 63(1):60–81, 2014.
  * [24] Gilbert Laporte, François Louveaux, and Hélène Mercure.  Models and exact solutions for a class of stochastic location-routing problems.  European Journal of Operational Research, 39:71–78, 1989.
  * [25] Kunlei Lian.  Service Consistency in Vehicle Routing.  PhD thesis, University of Arkansas, Fayetteville, 2017.
  * [26] Jens Lysgaard, Adam N. Letchford, and Richard W. Eglese.  A new branch-and-cut algorithm for the capacitated vehicle routing problem.  Mathematical Programming, 100(2):423–445, 2004.
  * [27] Diego Pecin, Artur Pessoa, Marcus Poggi, and Eduardo Uchoa.  Improved branch-cut-and-price for capacitated vehicle routing.  Mathematical Programming Computation, 9(1):61–100, 2017.
  * [28] Artur Pessoa, Marcus Poggi De Aragão, and Eduardo Uchoa.  Robust branch-cut-and-price algorithms for vehicle routing problems.  The vehicle routing problem: Latest advances and new challenges, 43:297–325, 2008.
  * [29] Nicola Secomandi.  A rollout policy for the vehicle routing problem with stochastic demands.  Operations Research, 49:796–802, 2001.
  * [30] Marius M. Solomon.  Algorithms for the vehicle routing and scheduling problems with time window constraints.  Operations Research, 35(2):254 – 265, 1987.
  * [31] Remy Spliet, Said Dabia, and Tom van Woensel.  The time window assignment vehicle routing problem with time-dependent travel times.  Transportation Science, 52(2):261–276, 2018.
  * [32] Remy Spliet and Guy Desaulniers.  The discrete time window assignment vehicle routing problem.  European Journal of Operational Research, 244(2):379–391, 2015\.
  * [33] Remy Spliet and Adriana F. Gabor.  The Time Window Assignment Vehicle Routing Problem.  Transportation Science, 49(4):721–731, 2015.
  * [34] Anirudh Subramanyam and Chrysanthos E. Gounaris.  A branch-and-cut framework for the consistent traveling salesman problem.  European Journal of Operational Research, 248(2):384–395, 2016\.
  * [35] Anirudh Subramanyam and Chrysanthos E. Gounaris.  A Decomposition Algorithm for the Consistent Traveling Salesman Problem with Vehicle Idling.  Transportation Science, 52(2):386–401, 2018.
  * [36] C.D. Tarantilis, F. Stavropoulou, and P.P. Repoussis.  A template-based Tabu Search algorithm for the Consistent Vehicle Routing Problem.  Expert Systems with Applications, 39(4):4233–4239, 2012.
  * [37] Marlin W. Ulmer and Barrett W. Thomas.  Enough waiting for the cable guy - estimating arrival times for service vehicle routing, 2017.  Available online at http://web.winforms.phil.tu-bs.de/paper/ulmer/Ulmer_SID.pdf.
  * [38] Anastasios D. Vareias, Panagiotis P. Repoussis, and Christos D. Tarantilis.  Assessing Customer Service Reliability in Route Planning with Self-Imposed Time Windows and Stochastic Travel Times.  Transportation Science, Articles in Advance, 2017.
  * [39] B.F.A. Vercammen.  Improving the planning accuracy and route efficiency at a freight distribution company: A case study at Van Opzeeland.  Master’s thesis, Department of Industrial Engineering & Innovation Sciences, Eindhoven University of Technology, 2016.
  * [40] Chengliang Zhang, George Nemhauser, Joel Sokol, Myun-Seok Cheon, and Dimitri Papageorgiou.  Robust Inventory Routing with Flexible Time Window Allocation, 2015\.  Available at Optimization Online. http://www.optimization-online.org/DB_HTML/2015/01/4744.html.