Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
The robust vehicle routing problem
(USC Thesis Other)
The robust vehicle routing problem
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
THE ROBUST VEHICLE ROUTING PROBLEM
by
Ilgaz Sungur
A Dissertation Presented to the
FACULTY OF THE GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(INDUSTRIAL AND SYSTEMS ENGINEERING)
August 2007
Copyright 2007 Ilgaz Sungur
Dedication
To the Burning Man Community
ii
Acknowledgements
I would like to first thank my external committee member Sven Koenig for being an invaluable
asset in my PhD career. His open-minded approach, understanding, and positive attitude as well
as his financial support were all far beyond my expectations.
Iwouldlike tothank mydepartmentchairJamesMooreforhis kindbehaviorandconstructive
approach to the suggestions I have made on the department. Not only did he listen but also took
certain supportive actions which improved the life of many others as well in their journey to PhD.
I would like to thank my co-adviser Maged Dessouky for being an intermediary between me
and my adviser and for being able to handle certain managerial decisions.
I would like to thank my adviser Fernando Ord´ o˜ nez for being able to adequately address all
the technical questions I have raised during the course of PhD.
I would like to thank my friends Zikri Altun, Caner Ta¸ skın and Evrim Dalkıran for their
financial support which made it possible for me to movefrom Istanbul to Los Angeles and to start
my PhD career smoothly.
I would like to also thank Burning Man for giving me a special admission to this year’s event.
As I promised, a copy of this dissertation will be delivered to the Temple to be shared with the
community before it is burnt.
Lastly, this researchis supported by NSF under grantCMS-0409887and by METRANS under
grant 06-11.
iii
Table of Contents
Dedication ii
Acknowledgements iii
List of Tables vi
List of Figures viii
Abstract ix
1 Introduction 1
1.1 Motivation 1
1.2 Proposed Methodology 3
1.2.1 Robust Optimization 4
1.2.2 Recourse Models 6
1.3 Structure of the Study 7
2 Literature Review 9
2.1 Deterministic Vehicle Routing Problem 9
2.1.1 Exact Solution Methods 10
2.1.2 Approximate Solution Methods 12
2.2 Stochastic Vehicle Routing Problem 13
2.2.1 Stochastic Optimization 13
2.2.2 Probabilistic Analysis 16
2.2.3 Real-time Routing 17
2.3 Robust Optimization 17
2.4 Courier Delivery Problem 20
3 Robustness of VRP Formulations 23
3.1 Methodology 23
3.2 VRP Formulations 26
3.2.1 Polynomial Vehicle Flow Models 27
3.2.2 Commodity Flow Models 29
3.2.3 Exponential Vehicle Flow Models 30
3.3 Experimental Results 31
3.3.1 Problem Size 31
3.3.2 Effect of C and K on Solution Time 32
3.3.3 Effect of the Travel Time Matrix on Solution Time 32
3.3.4 Understanding the Effect of CV
c
on Solution Time 35
iv
3.3.5 Effect of Geographical Distribution of Customers 39
3.4 Summary 40
4 A Priori Robust Optimization for VRP 42
4.1 Methodology 42
4.2 Robust Counterpart of VRP 44
4.2.1 Identifying the VRP Formulation 44
4.2.2 Deriving a General RVRP Formulation 47
4.3 Uncertainty in Demand 48
4.3.1 Derivation of the RVRP 48
4.3.2 Performance Measures and Bounds 52
4.4 Uncertainty in Demand and Travel Times 54
4.4.1 Derivation of the RVRP 54
4.4.2 Performance Measures and Bounds 56
4.5 Experimental Results 57
4.5.1 Robust versus Deterministic on Standard Problems 58
4.5.2 Robust versus Deterministic on Family of Clustered Instances 62
4.5.3 Robust versus Distributions of Excess Vehicle Capacity 69
4.6 Summary 72
5 Robust Optimization with Recourse for VRP 74
5.1 Methodology 74
5.2 Deterministic VRP Formulation 76
5.3 A Priori Robust VRP Formulation 78
5.3.1 Uncertainty in Service Times 79
5.3.2 Uncertainty in Time Windows 80
5.3.3 Probabilistic Customers 81
5.4 Recourse Model for Robust Counterpart of VRP 81
5.4.1 Recourse Action of Omitting Customers 82
5.4.2 Recourse Action of Rescheduling Routes 86
5.5 Real Life Problem Data and Uncertainty 90
5.6 Two-Phase Heuristic Solution Procedure 94
5.7 Experimental Results 99
5.7.1 Two-Phase versus Independent Daily Insertion 100
5.7.2 Two-Phase versus Real Life Solution 105
5.8 Summary 107
6 Conclusions and Future Research Directions 109
References 112
v
List of Tables
4.1 Augerat et al. (1995) sets 59
5.1 Comparison of two-phase with independent daily insertion 103
5.2 Comparison of two-phase with the real life solution 106
vi
List of Figures
3.1 LEFT: Average CPU time as a function of the vehicle capacity. RIGHT: Median,
1
st
and 3
rd
quartiles. With μ
c
=4.5, σ
c
= 3.375, K = 5, 50 random instances 33
3.2 Average CPU time as a function of mean travel time μ
c
. With C = 6, K = 5,
σ
c
= 3.375, 50 random instances 34
3.3 Average CPU time as a function of standard deviation of travel time σ
c
. With
C = 6, K = 5, μ
c
= 4.5, 50 random instances 34
3.4 AverageCPUtimeasafunctionofthecoefficientofvariationoftraveltimeσ
c
/μ
c
.
With C = 6, K =5, 50 random instances 35
3.5 Number of triangular inequalities violated and asymmetry as a function of the
coefficient of variation CV
c
. With C = 6, K = 5, 50 random instances 36
3.6 Initial LP gap as a function of triangular inequalities violated and solution time
as a function of initial LP gap for VRP1 and VRP3. With C = 6, K = 5 37
3.7 Average CPU time as a function of the asymmetry in the matrix of travel times.
With C = 6, K =5 38
3.8 Average CPU time as a function of the coefficient of variation CV
c
for symmetric
instances. With C =6, K = 5, 50 random instances 39
3.9 AverageCPUtime asafunction ofthe clusteringratior
c
forsymmetricinstances.
With C = 6, K =5, 50 random instances 40
4.1 Effect of lifting on solution time 46
4.2 Deterministic and robustsolutions for Augeratet al. (1995)set B instance 3. No.
of vehicles 5, capacity 100, value next to route is the deterministic demand served 61
4.3 Deterministic and robustsolutions for Augeratet al. (1995)set B instance 1. No.
of vehicles 5, capacity 100, value next to route is the deterministic demand served 62
4.4 Cluster generation in random sets 63
4.5 Comparison of deterministic and robust solutions for random set 1 64
4.6 Comparison of deterministic and robust solutions for random set 2 64
vii
4.7 Comparison of deterministic and robust solutions for random set 3 65
4.8 Deterministic and robust solutions for a bad case for robust from random set 1.
No. of vehicles 4, capacity 1500, value next to routes is the deterministic demand
served 66
4.9 Deterministic and robust solutions for a good case for robust from random set 3.
No. of vehicles 4, capacity 1500, value next to routes is the deterministic demand
served 67
4.10 Comparison of three uncertainty models for random set 3 68
4.11 Comparison of uniform buffer capacity and robust solutions for random set 1 70
4.12 Comparison of non-uniform buffer capacity and robust solutions for random set 1 71
5.1 Location of potential customers 91
5.2 A typical daily occurrence data of customers 91
5.3 Actual and fitted service time distributions 92
5.4 Actual and modified probability distributions of the occurrences of customers 92
viii
Abstract
We consider a Vehicle Routing Problem (VRP) under uncertainty and seek to find efficient solu-
tions in practice. Routing problems contain a fairly high degree of uncertainty in real life; thus
pre-planned optimal routes are no longer of practical use. Current methods to address the un-
certainty in VRP either make strong assumptions or increase the complexity of the deterministic
model significantly. In this study, we propose a routing model that uses robust optimization to
represent uncertainty in demand and travel times, we adapt this model for a real life courier de-
livery problem, we develop a heuristic for this large scale problem, and we evaluate the quality of
this solution procedure in practice.
To develop the routing model, we first experimentally compare the robustness of VRP formu-
lations. We then adapt the most suitable VRP formulation and develop a robust counterpart of
VRP (RVRP) with demand and travel time uncertainty that can be solved efficiently. The RVRP
is not more difficult to solve than VRP; in fact the RVRP with demand uncertainty is simply a
single instance of VRP with modified demand data. We show experimentally that the solution of
this RVRP when compared to the solution of the deterministic VRP protects against the unmet
demand at the expense of slightly increased cost, especially when the network is clustered.
We combine robust optimization with scenario based stochastic optimization to address a real
life courier delivery problem, a variant of VRP with time windows. We use robust optimization
for uncertainty in service times and time windows, and we represent uncertainty in demand using
stochastic optimization with recourse. Our model maximizes number of customers serviced and
similarity of daily routes while minimizing route lengths and total penalty. We develop a two-
phaseinsertionbasedheuristicbalancingtheseobjectives. Ourexperimentsshowthatourheuristic
improvessimilarityofroutesandtotalpenaltyattheexpenseofincreasingroutelengthscompared
toindependentdeterministicsolutionsforeachdemandoutcome. Oursolutionforreallifeproblem
is also better than the current practice.
ix
Chapter 1
Introduction
1.1 Motivation
Many industrial applications deal with the problem of routing a fleet of vehicles from a depot to
service a set of customers that are geographically dispersed. This type of problem is faced daily
by courier services (e.g., Federal Express, United Parcel Service, and Overnight United States
Postal Service), local trucking companies, and demand responsive transportation services, just to
name a few. These types of services have experienced tremendous growth in recent years. For
example, both United Parcel Service (UPS) and Federal Express have annual revenue of well over
$10 billion, and the dial-a-ride service for the disabled and handicapped is today a $1.2 billion
industry (Palmer, Dessouky, and Abdelmaguid 2004). However, congestion and variability in
demand and travel times affect these industries on three major service dimensions: travel time,
reliability, and cost (Meyer 1996). Therefore, there is a need to develop routing and scheduling
tools that directly account for the uncertainty in demand and varying congestion in the road
network for this vital industry.
Problems where a given set of vehicles with finite capacity have to be routed to satisfy a
geographicallydisperseddemandatminimumcostareknownasVehicleRoutingProblems(VRP).
ThisclassofproblemswasintroducedbyDantzigandRamser(1959)andgeneralsurveysofvehicle
routingresearchcanbefoundinTothandVigo(2002a),Fisher(1995),LaporteandOsman(1995),
1
and Golden and Assad (1988). During the last two decades, the VRP and a wealth of extensions
have been the subject of a wide body of research (Bodin, Golden, Assad, and Ball 1983; Laporte
1992; Fisher 1995; Toth and Vigo 2002a). Most of the solution methods developed for the VRP
assumethatthe inputdataisknown. Forexample, the customerdemands, travelcosts, and travel
times are deterministic and known a priori. The built-in assumption of these approaches is that
there will be only small deviations on the realization of the demand and travel times from the
plan so that the pre-determined routes form a basis for either the pickup or delivery schedule.
In the real world, however, operations in any traffic network contain a fairly high degree of
uncertainty, which includes arrival of new orders, cancellation of existing orders, variable waiting
timesandvariabletraveltimesdue totrafficcongestion. Thus, inahighlystochasticenvironment,
the pre-planned optimal routes are no longer of practical use and there is a need for models
addressing this uncertainty.
There exists a vast literature on models and algorithms that take into account uncertainty
in routing problems. The most studied areas in the stochastic vehicle routing problem (SVRP)
literature have been the VRP with stochastic demands, and with stochastic customers (Bertsimas
1992; Gendreau, Laporte, and Seguin 1995; Gendreau, Laporte, and Seguin 1996; Bertsimas and
Simchi-Levi 1996; Secomandi 2001). Despite its importance and practical application, especially
in major cities with traffic congestion, research efforts on the VRP with stochastic travel times
have been limited (Carraway,Morin, and Moskowitz1989; Laporte, Louveaux, and Mercure 1992;
Lambert, Laporte, and Louveaux 1993).
Generally speaking, current methods to address the uncertainty in SVRP consider one of the
following strategies (or both): they either make strong assumptions regarding the distribution of
theuncertainparameters,orthemethodsapproximatetheSVRPwithamuchlargerdeterministic
modelinwhichtheuncertaintyisrepresentedthroughscenarios(BertsimasandSimchi-Levi1996).
The objective of these methods is to obtain a solution that minimizes the expected value, or
to analyze the performance of a routing policy in expected value or worst case. The resulting
solution is potentially sensitive to the actual data that occurs in the problem. Considering that
2
the VRP solution in a given application will only face a single realization of the uncertainty, a
morereasonablegoalmightbe to obtain a “robust”solution, i.e. a “good”solution forall possible
data uncertainty.
In this study, we address the uncertainty in routing vehicles by following the robust opti-
mization methodology. We propose to obtain a robust solution for all possible realizations of the
uncertainty. We introduce uncertainty in several different problem parameters. We first develop
an a priori robust optimization model and we also analyze the sensitivity of the model and effi-
ciency of the robust solution by using exact solution procedures. Lastly, we develop a recourse
robust optimization model for a specific large scale real life problem. We also propose a heuristic
solution procedure for this courier delivery problem in industry.
1.2 Proposed Methodology
We propose to directly obtain a robust solution for the VRP, where by robust solution we mean
a solution that is fixed and efficient for all realizations of the problem uncertainty. We propose
to achieve this by introducing a new problem formulation, the Robust Vehicle Routing Problem
(RVRP), whose goal is a solution that optimizes the worst case value over all data uncertainty.
Robust solutions have the potential to be viable solutions in practice, since they tend not to be
far from the optimal solution and significantly outperform the optimal solution in the worst case
(Goldfarb and Iyengar 2003).
A well-knownmethod of addressinguncertainty through scenariosis recoursein which some of
the variables are fixed for all possible realizations of the uncertainty (first-stage) and some others
take values based on a particular realization of the uncertainty (second-stage) whereas in a priori
optimization there is only one stage. Although this method gives a flexible and adaptive solution,
it usually requires additional assumptions and it increases the complexity of the formulation.
Therefore,generallyspeakinganapriorirobustmodelbenefitsfrom“simplicity”intheformulation
and the setof assumptions at the expense of losing this “flexibility”of recourseactions. In certain
3
applications, however recourse actions are necessary; and thus we also extend our a priori robust
model to problems with recourse by incorporating scenario based recourse models when applying
the methodology to a real life courier delivery problem.
In this study, we will develop efficient routing strategies for routing problems with uncertainty
and validate them on realistic instances. We will first adapt the robust optimization approach to
consider the uncertainty in several problem parameters. This will require identifying an efficient
VRP formulation, which is robust, and deriving the robust counterpart of the VRP, namely the
RVRP. Then we will combine the RVRP model with scenario based stochastic models to consider
recourseactionsforalargescalecourierdeliveryprobleminindustryforwhichwewillalsodevelop
a heuristic solution procedure. In the next subsection, we elaborate more on robust optimization
and in the last subsection more on the scenario based stochastic optimization with recourse.
1.2.1 Robust Optimization
The term “robust” has been used in the VRP literature in two different contexts. First, research
has aimed at obtaining solution methodologies that are robust, which means that the method
obtains optimal or close to optimal solutions for many different problem instances of the VRP
(Haughton 1998;Wassanand Osman2002). The second useofrobustnessin the VRP literatureis
closer to the one that concerns us. Within the SVRP, the solution obtained through an algorithm
or policy is termed robust if its worst-case value over the uncertainties considered, is not far-off
the optimal value (Bertsimas and Simchi-Levi 1996). Our definition of robustness differs from
these earlier cases since it directly takes into consideration all possible realizations of the problem
uncertaintyindeterminingthesolutioninsteadofsimplyconsideringthesensitivityofthesolution
approach to various problem instances.
ToformulateandsolvetherobustversionoftheVRP,orRVRP,wewillfollowthecurrentliter-
ature in mathematical programmingwhich defines the robustcounterpartfor convexoptimization
problems (Ben-Tal and Nemirovski 1998; Ben-Tal and Nemirovski 1999; Ben-Tal, El-Ghaoui, and
Nemirovski 2000) and has been recently extended to integer programming problems (Bertsimas
4
and Sim 2004; Bertsimas and Sim 2003). The general approach of robust optimization is to op-
timize against the worst instance that might arise due to data uncertainty by using a min-max
objective. The resulting solution from the robust counterpart problem is insensitive to the data
uncertainty, as it is the one that minimizes the worst case, and therefore is “immunized” against
this uncertainty. The robust optimization methodology assumes the uncertain parameters belong
to a bounded uncertainty set. Clearly the size of the uncertainty set influences the deviation
from optimality of the robust solution. For fairly general uncertainty sets, the resulting robust
counterpartproblem is modestly largerthan the originalproblem, and thereforehas a comparable
complexity.
Inourrobustoptimizationmethodology,weareinterestedinintroducinguncertaintyinseveral
parameters of VRP (i.e. demand, travel time, service time, time windows) and formulating a
tractable RVRP which can be solved efficiently. From a practitioner’s point of view, the challenge
of developing and solving RVRP is threefold:
1. identifying a VRP formulation which is robust with respect to the problem parameters and
which is not inherently difficult to solve,
2. formulatingtherobustcounterpartofthisVRPformulationbasedonmeaningfuluncertainty
sets, which should yield an amenable RVRP formulation,
3. and developing an efficient solution procedure for the RVRP formulation.
Since the mathematical derivation based on the uncertainty set will not reduce the complexity
of the VRP, the resulting RVRP will be at least as complex as the VRP. Therefore the first
challenge is to adopt a “suitable” VRP formulation which will not further complicate the RVRP
formulation, and which will yield an amenable RVRP for computation. Here by “suitable”, we
refer to the complexity of the VRP formulation. The running time complexity of a problem is
affected by three aspects: the problem data, the problem formulation and the solution procedure.
Different problem formulations lead to different complexities, and we refer to this as the inherent
difficulty of the problem formulation. We are interested in a VRP formulation which makes the
5
problem easier to solve in the following sense: the running time complexity of the problem due to
the formulation itself should be reasonable.
The second challenge is to construct “meaningful” uncertainty sets for the VRP to derive
RVRP accordingly. Here by “meaningful”, we refer to two different aspects: the definition of
the uncertainty set should be closely related to real life applications so that a practitioner could
construct this uncertainty set with ease from real life data, and the uncertainty set should yield a
tractable robust counterpart.
Finally, when the RVRP is derived from a suitable VRP by using meaningful uncertainty sets,
thecomputationalaspectofsolvingtheRVRPwillbethethirdchallenge. Eventhoughatractable
RVRP might be obtained, the problem will still be NP-Hard and an efficient solution procedure
will be necessary to develop in order to reduce the running time so that real life applications can
be solved.
1.2.2 Recourse Models
In the stochastic optimization literature, when the uncertainty is represented through scenarios,
a well-known method to address this uncertainty is to consider recourse actions for each scenario.
Generallyspeaking,thistypeofmodelsincludesapriorivariableswhosevaluesarefixedinthefirst
stage and adjusted variables whose values are determined in the second stage based on the values
of a priori variables and the outcome of the uncertainty. Therefore the model can be considered
as an extension of the a priori model which now includes adjusted variables as a secondary stage.
The purpose of recourse models is to provide specific solutions for each realization of the un-
certainty via adjusted variables for each scenario. How these second stage variables are adjusted
is determined by arecourseactionwhich is usually problemdomain specific. This definition ofthe
recourse action should be properly translated into the formulation which necessitates additional
modeling efforts and/or probabilistic calculations of the expected outcome of the recurse action.
Overall, the flexibility of a recourse model in obtaining a specific solution for each scenario comes
6
with the cost of increased complexity in the formulation compared to an a priori model. In par-
ticular, as the number of scenarios increases, the difficulty of solving the recourse model becomes
more pronounced.
For a given problem domain, the four main challenges a practitioner should address in order
to formulate a recourse model are therefore:
1. identifying an efficient a priori model,
2. identifying a realistic and practical recourse action for the problem domain,
3. formulating an efficient recourse model based on the recourse action and the a priori model,
4. and developing an efficient solution procedure for the recourse model.
In this study, after developing an a priori RVRP formulation for a real life courier delivery
problem which is a variant of VRP with time windows, we will be interested in incorporating
scenario based stochastic optimization with recourse into the formulation due to the probabilis-
tic nature of the customers, which requires adaptable routes for each possible realization of the
uncertainty in customer occurrences. For this, we first have to come up with an efficient a priori
RVRP model. Second, we must identify a problem domain specific recourse action which should
be realistic and which can be implemented with ease both in practice and in the formulation. The
third challenge is to incorporate the recourse action into the a priori model to obtain an efficient
recourse model. Lastly, we develop an efficient solution procedure to be able to solve a realistic
courier delivery problem whose complexity is further increased due to the recourse formulation.
1.3 Structure of the Study
The literaturereview in Chapter 2providesa broaderdiscussionofVRP, uncertaintyin VRP, and
robustoptimization. In particular,we exploreefficientexactandapproximatesolutionprocedures
for VRP, we contrasta priori and recourse models for uncertainty in VRP, and discuss the courier
delivery problem.
7
Chapter 3 addresses the question of identifying an efficient VRP formulation for robust opti-
mization. Our experimental analysis suggests that the VRP formulation with polynomial number
of subtour elimination constraints exhibits a more robust behavior compared to the others when
solved by a general integer programming solver.
In Chapter 4, taking into consideration the nature of the uncertainty in cost of travel times
and demands, we formulate a RVRP amenable for computation. The possible conservative nature
of a robust solution that has to be efficient for all uncertainties, and thus its utility in practice, is
explored experimentally by investigating its dependence on the uncertainty sets when it is solved
by a competitive VRP solver in the literature. Our results for the case of demand uncertainty
reveals that the robust solution is preferable when the network has clustered zones far from the
depot and a dense random zone near the depot. We also show that the robust optimization
provides a better strategy in distributing the excess vehicle capacity compared to uniform and
non-uniform distribution strategies.
In Chapter 5, we consider a real life courier delivery problem in industry and develop a RVRP
with uncertain service times and time windows. We also consider the demand locations to be
probabilistic for which we incorporate scenario based recourse formulations. Lastly, we propose a
two-phase insertion based heuristic to solve such a large scale real life problem. Our experiments
show the sensitivity of the heuristic on uncertain parameters and contrast the two-phase solution
with the solution obtained when all the scenarios are solved independently. We also show that
our two-phase solution provides improvements over the reallife solution of the particular problem
considered in this study.
Lastly, Chapter 6 concludes the work and indicates future research directions.
8
Chapter 2
Literature Review
We break down the literature review in four sections. We discuss first the different formulations
and solution methods for the deterministic VRP, second approaches to treat uncertainty in the
VRP, third the literatureon robustoptimization, and fourththe courierdeliveryproblem. In each
section we point out where our work fits in and what deficiencies in the literature it addresses.
2.1 Deterministic Vehicle Routing Problem
Due to the wide rangeof applicationsthat the VRP has addressed, there are a number ofdifferent
variants of the VRP, obtained by adding specific features or constraints to the basic capacitated
vehicle routing problem (CVRP). Some of the most important variants are: VRP with time
windows (VRPTW), VRP with multiple depots (MDVRP), VRP with back-hauls (VRPB), and
VRP with pick-up and deliveries (VRPPD). We will limit our study to the CVRP in the next two
chapters and then extend it to VRPTW in the last chapter.
The research on solution procedures in the deterministic VRP can be separated into exact
algorithms and heuristics. Although the optimization solution techniques can be used to find
optimalsolutionsonlyforsmalltomediumsizeproblems,theycanalsobethebasisfordetermining
efficient approximate solutions for larger problems (Koskosidis, Powell, and Solomon 1992). In
the next subsection, we reviewthe literature on exactsolution approachesand then briefly discuss
9
the various approximate solution approaches in the last subsection. Even though the size of the
problem and the nature of the objective function is different for the RVRP, we will still be able to
draw from the existing body of work on finding both optimal and approximate solutions for the
VRP, to our robust problem.
Since the VRP polynomially reduces to the Traveling Salesman Problem (TSP), the com-
plexity of the VRP is the same as the complexity of the TSP. Therefore for general d
ij
, the
distance between locations i and j, the problem is NP-hard and there is no known approximation
scheme with a guaranteed worst case bound (Orponen and Mannila 1987). If d
ij
is a metric, then
there exists a polynomial time algorithm with a guaranteed worst case performance ratio of
3
2
(Christofides 1976). Finally if d
ij
is the Euclidean distance, then the problem has a polynomial
time approximation scheme (Arora 1996).
It is well-known that the VRP is NP-Hard and the robust counterpart of the VRP will be at
leastasdifficult tosolve. However,wewill benefit from the currentbody ofresearchaboutsolving
the VRP when it comes to solving the RVRP. The nature ofthe formulation of the RVRP and the
related definition of the uncertainty sets will play an important role in the complexity of solving
the RVRP, either by an exact or approximate solution procedure.
2.1.1 Exact Solution Methods
The exact algorithms can be further classified into dynamic programming (DP) and integer pro-
gramming (IP) methods; see Laporte (1992). Since there has been recent advances in solving
largescaleIPproblems, mostofthe recentexactsolution approacheshavefocused onIP methods.
Furthermore, IP models can be the basis for developing good approximate solutions. Hence, our
research approach will be based on an IP formulation (see next chapter) and is the focus of our
review in this section. In terms of the DP approaches, Eilon et al. (1971) were the first to apply
dynamic programming to the VRP. An efficient way to reduce the state space was introduced by
Christofides et al. (1981b). Psaraftis (1980, 1983) applied the DP approach to the dial-a-ride
problem.
10
Therearetwogeneralclassesofintegerprogrammingsolutionmethods, whichaddressdifferent
formulations of the VRP. Classical branch and bound methods address the more natural, or edge
basedformulation,wherethereareintegervariablesassigningavehicletoeachedgeinthenetwork.
The work by Christofides et al. (1981a) uses lower bounds obtained from spanning trees of the
networkgraph. Other competitivebranchand bound methods aredue to Laporteetal. (1985)for
the symmetric distance case and Laporte et al. (1992)and Achuthan et al. (1996)for asymmetric
distance. Recently,adistributedbranch-and-cutsearchmethod,whichemployedanLPbasedTSP
code and specialized cutting planes, solved a 120 customer, 4 vehicle VRP, and proved optimality
ofthesolutionsettlingtheWhizzkids’96mathematicschallenge(Applegate,Cook,Dash,andRohe
2002). Some other recent improvements of branch-and-cut algorithm are Lysgaard et al. (2004),
Baldacci et al. (2004), Ralphs et al. (2003), and Ralphs (2003) for a parallel implementation. A
novel approach, branch-and-cut-and-price algorithm, is due to Fukasawa et al. (2006). Among
those, the implementation of Ralphs et al. (2003) is competitive and the open source solver can
be found in the SYMPHONY Library at the webpage http://branchandcut.org/VRP.
AnalternateformulationoftheVRPisthepathbasedformulation. Thisformulationconsiders
a binary variable for each vehicle and each possible route, thus generating an exponential number
of integer variables. The natural method to attack this formulation is column generation. This
formulation for the VRP is due to Balinski and Quandt (1964). Column generation techniques
were applied to this problem by Rao and Zionts (1968), Foster and Ryan (1976), Orloff (1976),
and Desrosiers et al. (1984).
Therehasbeen furtherworkin both solutionapproaches,butithasconcentratedonparticular
extensions to the VRP that are important in applications. The two principal extensions are the
VRP with time windows and the Pickup and Delivery problem. The resulting algorithms exploit
the structure of these specific problems, and in some cases do not readily extend to the general
VRP.
11
2.1.2 Approximate Solution Methods
Cordeau et al. (2002) and Laporte et al. (2000) provide an excellent survey of the different
approximatesolutionapproachestotheVRPandtheVRPwithtimewindows. Thereareavariety
of different heuristics that have been applied to the VRP. Briefly, they are a savings algorithm
by Clarke and Wright (1964), sweep algorithm (Wren 1971; Wren and Holliday 1972; Gillett
and Miller 1974), two-phase algorithm (Christofides, Mingozzi, and Toth 1979), partitioning and
partial solution (Fisher and Jaikumar 1981), simulated annealing (Alfa, Heragu, and Chen 1991;
Osman 1993), genetic algorithm (Potvin and Bengio. 1996), evolutionary algorithms (Homberger
and Gehring 1999), tabu search(Willard 1989; Gendreau, Hertz, and Laporte 1991; Taillard1993;
Taillard, Badeau, Gendreau, Guertin, and Potvin 1997),and constraintprogramming(De Backer,
Furnon, Kilby, Prosser, and Shaw 2000; Caseau and Laburthe 1999).
An alternative approach is to develop an optimization-based heuristic. They are also called
mathematical programming-based heuristics that include algorithms that are directly based on
the mathematical programming formulation of the underlying routing problem. These heuristics
require more computational effort, but produce high-quality solutions. Sexton and Choi (1986)
minimized the total vehicle operating time and total customer penalty due to missing any of the
time windows for the Pickup and Delivery Problem by applying a Bender’s decomposition proce-
dure to a mixed binary nonlinear formulation. They used a heuristic to solve the computationally
demanding routing problem. Koskosidiset al. (1992)formulated the vehicle routing problem with
soft time windows as a mixed integer program and heuristically decomposed the problem into an
assignment/clustering subproblem and a series of routing and scheduling subproblems. Savels-
bergh and Sol (1998) used an approximation algorithm based on the branch-and-price algorithm
proposed by Dumas et al. (1991) to solve a large-scale Pickup and Delivery Problem with hard
time windows. The dynamic programfor column generation in Dumas et al. (1991)is substituted
with a heuristic construction and insertion algorithm.
12
On the other hand, optimization-based heuristics are not viable solutions for large scale prob-
lems due to computation time, for which fast and efficient heuristics are preferred. One pop-
ular heuristic is insertion and recent examples are Quadrifoglio and Dessouky (2007), and Lu
and Dessouky (2006). Solomon (1987) compares different classical insertion methods with other
heuristics on well-knownbenchmark problems. In a more recent study, Cordeau et al. (2002)uses
asimilarmethod tocompareseveralheuristicsonlargeprobleminstancesbasedon fourmeasures:
accuracy, speed, simplicity and flexibility. Insertion is still one of the widely used and promising
heuristics for such large scale problems that we will also benefit in our research when solving a
large scale real life problem.
2.2 Stochastic Vehicle Routing Problem
There are many problem parameters of the VRP that in practice have a significant degree of
uncertainty. For example in the case of pickups, load sizes and travel times are typically unknown
with any certainty prior to arriving at the stop. The stochastic vehicle routing problem attempts
to construct a routing policy (as opposed to a fixed route) that minimizes expected cost, given
this inherent randomness. Stochastic Vehicle Routing Problems (SVRP) arise whenever some or
all elementsofthe VRP problemarerandom. Solution approachesforthe SVRP canbe separated
into three categories: stochastic optimization, probabilistic analysis, and real-time routing.
2.2.1 Stochastic Optimization
Stochastic programmingwas introduced by Dantzig (1955), and good referencesare Ermolievand
Wets (1988) and Mulvey et al. (1995). Extensions to stochastic integer programming appear in
Schultz et al. (1998), and Takriti and Birge (2000).
The most studied areas in the stochastic vehicle routing problem literature have been the
VRP with stochastic demands (VRPSD), and with stochastic customers (VRPSC). A major con-
tribution to the study of the VRPSD comes from Bertsimas (1992). This work illustrates the a
13
priori method with different recourse policies (re-optimization is allowed) to solve the VRPSD
and derives several bounds, asymptotic results and other theoretical properties. Bertsimas and
Simchi-Levi (1996) surveys the development in the VRPSD with an emphasis on the insights
gained and on the algorithms proposed. Besides the conventional stochastic programming frame-
work, a Markov decision process for single stage and multistage stochastic models are introduced
to investigate the VRPSD in Dror et al. (1989)and Dror (1993). More recently, a re-optimization
type routing policy for the VRPSD is introduced by Secomandi (2001).
The VRPSC, in which customers with deterministic demands and a probability p
i
of being
present, and the VRP with stochastic customers and demands (VRPSCD), which combines the
VRPSCandVRPSD,firstappearintheliteratureofJ´ ez´ equel(1985),Jaillet(1987)andJailletand
Odoni (1988). Bertsimas (1988) gives a more systematic analysis and presents several properties,
bounds and heuristics. Gendreau et al. (1995) proposes the first exact solution, an L-shaped
method, and a meta-heuristic, a tabu search for the VRPSCD.
The biggest drawback of the stochastic optimization solution approach is the fact that rep-
resenting the distribution of the uncertainty is typically done through discrete scenarios, which
increasesthe dimension of the problem tremendously. Additionally the exactdistribution is rarely
known in practiceand mostsolutionmethods forthe VRP with stochastictraveltimes requirethe
knowledgeof the distribution of the sum ofthe traveland service times along the routes (Laporte,
Louveaux,and Mercure1992; Gendreau, Laporte, and Seguin 1996). However,sucha distribution
cannot be easily computed in the presence of hard time windows.
On the other hand, recourse models considering variants of VRP including the single-vehicle
VRP(TSP)withstochasticdemandand/orstochasticcustomershavebeeninvestigatedin several
studies; a good summary is provided by Gendreau et al. (1996). Most of the work has been
dedicated to direct probabilistic computation of the expected cost of recourse action. TSP with
probabilistic customers is first introduced in Jaillet (1988). The idea of a priori optimization
based on the recourse action of omitting non-occurring customers is defined and a probabilistic
method of calculating the expected tour length is described. Later, Bertsimas et al. (1990)
14
extends this work to other applications of TSP, and provides insights for VRP. They define two
recourse actions for VRP: follow the same sequence of customers in the a priori route but omit
non-occurring customers; follow exactly the a priori route until the capacity is depleted in which
case go back to depot to replenish and resume the route. Waters (1989) gives models for several
variantsofVRPandsuggeststhreerecourseactions: followexactlytheaprioriroute(fixed-route),
follow the same sequence of customers in the a priori route but omit non-occurring customers
(semi-fixed route), and first omit all the non-occurring customers then completely reschedule
the recourse routes (variable route). This work also provides an experimental study to compare
these recourse actions and suggests a break-even point for rescheduling the routes with respect
to the number of non-occurring customers. Swihart and Papastavrou (1999) develops a queuing
theory approach for single vehicle pick-up and delivery problem where customers appear in the
system dynamically with stochastic locations. They propose three routing policies for which
they provide analytical and simulation results. Laporte et al. (2002) is a recent study of VRP
with stochastic demand using the recourse action of resuming the route after replenishing the
vehicle at the depot. The probabilistic computation is used to approximate the expected tour
length and an exact solution procedure based on Integer L-Shaped Method is devised. Bent
and Hentenryck (2004) combines VRP with stochastic and dynamic customers in which some of
the customers arrive dynamically and there is available information about the sampling of the
customers. They develop a solution procedure to maximize the total number of customers to be
serviced in which they maintain a pool of feasible solutions which service current known requests
and accommodate possible future requests. Morales (2006) uses robust optimization approach
for the VRP with stochastic demand and adopts the recourse action of resuming the route after
replenishing the vehicle at the depot. The study considers three types of uncertainty: no advance
demand information, total demand visibility, and one customer ahead demand visibility. For a
given route, the determination of the worst-casevalue for the recourse action is reduced to finding
the longest path on an augmented network. A tabu search heuristic is developed and the robust
solution is compared to the traditional expectation minimization approaches.
15
In most of the recourse literature, either the original deterministic model is significantly aug-
mentedwhileformulatingtherecourseactionthroughscenariosortheexpectedcostoftherecourse
actioniscalculatedbasedonalargesetofassumptionsabouttheuncertainty. Inbothofthecases,
the flexibility of having an adaptive recourse solution comes with an increased complexity in the
formulation. Our a priori robust optimization approach gives away this flexibility by providing
a fixed “good” solution for all possible realizations of the uncertainty to benefit from the “sim-
plicity” of the formulation and the set of assumptions. However, we will combine later robust
optimization and scenario based stochastic optimization to develop “simple and flexible enough”
models with recourse.
Compared with stochastic customers and demands, the VRP with stochastic service time and
travel time (VRPSSTT) has received less attention. Laporte et al. (1992) proposes three models
forthe VRPSSTT: chanceconstrainedmodel, 3-indexrecoursemodeland 2-indexrecoursemodel,
andpresentsageneralbranch-and-cutalgorithmforall3models. TheVRPSSTTmodelisapplied
to a banking context and an adaptation of the savings algorithm is used in the work of Lambert
et al. (1993). Jula et al. (2006) develops a procedure to estimate the arrival time to the nodes in
the presence of hard time windows. In addition, they use these estimates embedded in a dynamic
programming algorithm to determine the optimal routes.
2.2.2 Probabilistic Analysis
This solution strategy for the SVRP provides analysis of different algorithms or policies for the
VRP in expected value and worst case for given assumptions on the distributions of uncertain
parameters and the behavior of the system (e.g. heavy traffic). References for this method are
the work by Federgruen and Simchi-Levi (1995), Powell et al. (1995), and the survey paper by
Bertsimas and Simchi-Levi (1996).
This solution approachhowever does not provide the means to obtain the best solution for the
SVRP, simply the means to analyze the effectiveness of a proposed solution. This effectiveness is
undoubtedly tied with the assumptions on distributions and network behavior.
16
2.2.3 Real-time Routing
A drastically different method of addressing the uncertainty present in applications of the VRP
is to develop efficient strategies to update a solution to meet changes in demand, customers, and
travel times.
Real-time routing has only received limited attention in the literature and a general overview
can be found in Psaraftis (1995). Until recently the research on real-time routing was limited
to the analytical demand-responsive research (Daganzo 1978; Daganzo 1984; Stein 1978). A
novel direction of research is the work by Bertsimas and Van Ryzin (1991, 1993), Bertsimas and
Simchi-Levi(1996),andKamounandHall(1996)whichprovideapproximationandboundsonthe
performance of the updating strategy for various traffic, network, and uncertainty assumptions.
To date, no one has developed a viable routing product that can optimally (or near optimally)
route vehicles in real-time (Hall and Partyka 1997).
It is important to note that research on real-time routing does not preclude from attempting
to obtain an a priori solution that is close to optimal, regardless of the resulting uncertainty. The
idea being that if we start from a robust solution, it would provide the conditions to guarantee
better bounds on the performance of a real-time routing algorithm.
2.3 Robust Optimization
Thecurrentrobustoptimizationresearchiscloselyrelatedtorobustcontroltheory;seeforexample
(Zhou, Doyle, and Glover 1996). Robust optimization, or minimizing the worst case of a set of
parametrized functions, can be traced to the work by Soyster (1973), but the recent activity
in the field started with Ben-Tal and Nemirovski (1998, 1999, 2000a), Ben-Tal et al. (2000),
and El-Ghaoui et al. (1998). Renewed interest in the robust optimization methodology can be
attributed to recent developments in interior-point methods, that allow the efficient solution of
the robust problem. For instance, for a linear programming problem and reasonable conditions
on the uncertainty, the robust solution is obtained by solving a Second Order Cone Program
17
(SOCP) (Ben-Tal and Nemirovski 1998), which can be solved to the desired accuracy efficiently
using interior-point methods.
The general approach of robust optimization is to optimize against the worst instance that
might arise due to data uncertainty by using a min-max objective. The resulting solution from
therobustcounterpartproblemisinsensitivetothedatauncertainty,asitistheonethatminimizes
the worst case and that is feasible under all outcomes of the data uncertainty, and therefore is
immunized against this uncertainty. Robust solutions have the potential to be viable solutions
in practice, since they tend not to be far from the optimal solution of the deterministic problem
and significantly outperform the deterministic optimal solution in the worst case realization of
the uncertainty (Goldfarb and Iyengar 2003; Bertsimas and Sim 2004). Here, robust optimization
methodology assumes that the problem remains feasible for any outcome of the data. If in a
particular application, a given realization results in infeasibility, one should address this problem
with a different methodology and a different model.
Robust optimization has been applied to different problem domains. For example, El-Ghaoui
and Lebret (1997) and Pinar (2002) address least-square problems and data fitting. In Ben-Tal
and Nemirovski (2000b), the authors introduce the robust optimization version of the structural
truss design problem, which turns out to be a semi-definite program. There has also been recent
workapplyingrobustoptimizationtotheportfolioselectionproblems(GoldfarbandIyengar2003;
Ben-Tal, Margalit, and Nemirovski 2000; El-Ghaoui, Oks, and Oustry 2003; Lutgens and Sturm
2002).
An important conclusion of the robust optimization methodology is that, for an appropriate
uncertainty set, the robust counterpart problem is a tractable problem, however it can have an
increasedcomplexity. For example the robustcounterpartproblemof a convexquadraticprogram
is a SOCP if the uncertainty set is ellipsoidal, but it becomes NP-hard if the uncertainty set is
defined by interval bounds; see Ben-Tal and Nemirovski (1998).
TherecentworkbyBertsimasand Sim (2001,2002)expandsthe robustoptimizationapproach
tointegerprogrammingproblems. Theyshowthatiftheuncertaintyisinthecostcoefficientsorin
18
the matrixcoefficients, then the resultingrobustcounterpartproblem is only modestly larger,and
consequently of the same complexity than the original problem. Another significant contribution
is the work by Ben-Tal et al. (2003) and Guslitser (2002), which extends robust optimization to
problems which have decision variables that adapt to the uncertainty (problems with recourse).
Since a robust solution optimizes the worst case over the uncertainty set, it is bound to have a
worseobjective value than an optimal solution for the deterministic data. This difference between
the robust solution and the optimal solution provides a measure of the conservativeness of the
robust solution, and therefore of its significance for the problem in practice, for the additional
level of robustness. The first computational study of the difference between robust and optimal
solutions appears in Ben-Tal and Nemirovski (2000a). This work considers linear programs and
assumes interval uncertainty around all the “ugly” coefficients (coefficients with long decimal
places). The conclusion of this study is that the resulting robust solution typically did not cost
much in terms of optimality. The same sort of conclusion is observed in Goldfarb and Iyengar
(2003),whereforaportfoliooptimizationproblemtheauthorsobtainrobustsolutionsthatachieve
80% of the optimal solution while being 200% better than the optimal solution in the worse case.
Efficient robust solutions are also obtained for integer programming in Bertsimas and Sim (2004).
This work in addition defines a budget of uncertainty, with which the user can set how much
uncertainty to consider in the problem, and thus enables the study of how a solution deteriorates
from optimality as it becomes more robust. These works suggest that, in many settings, a robust
solution with objective value close to optimal exists.
Other alternatives to obtain robust solutions, in the sense that they are resistant to problem
data variation, for integerprogrammingproblems havebeen explored. In Kouvelis and Yu (1997),
the procedure that considers the robust problem causes an increase in the complexity of the
problem. Other methods such as Averbakh and Berman (2000), Averbakh (2001), and Yaman
et al. (2001) are combinatorially oriented and do not seem to be easily generalizable to other
integer optimization problems.
19
Weconcludethattherobustoptimizationmethodology,atleastintheory,couldleadtoarobust
counterpart of the VRP that is tractable since for reasonable assumptions on the uncertainty,
the robust counterpart is only modestly larger than the original problem. A robust version of
the VRP has not been fully explored in the literature, and the potential of obtaining robust
solutions, which will perform well for any realization of the uncertainty could lead to better
practical implementations for the VRP.
2.4 Courier Delivery Problem
The courier delivery problem which concerns routing of daily courier services is frequently faced
in the industry. For example, both UPS and Federal Express have annual revenue of well over
$10 billion. This problem is a variant of VRPTW and has particular characteristics compared
to a general VRP. In this industry, customers are usually located in an urban area and travel
distances are fairly small in the city and therefore have less significance. On the other hand, the
number of potential customers who can make a delivery request each day is very large. Therefore
the couriers are to make frequent and short stops every day. Each delivery to a customer site
is associated with a service time and it must be done within a specific interval of time given
by time windows which can be either hard constraints or soft constraints allowing penalty for
deviations. The problem could be either capacitated or uncapacitated depending on the type of
the goods delivered. Uncertainty arises most of the time due to the unknown service times and
to the probabilistic nature of the customers, i.e. daily delivery requests from potential customers
are not known beforehand but they usually become available in the morning. The challenge is to
construct a daily schedule of couriers to cover all the delivery requests at minimum cost. Mostly,
it is also desired to have similarity between daily schedules to eliminate significant re-planning in
each day and to increase the familiarity of drivers with daily routes.
The particular real life courier delivery problem we will consider later in this study is an unca-
pacitatedVRPTWwith uncertaintyinservicetimes andtimewindowstogetherwith probabilistic
20
customers. The delivery requests become available each morning and it is desired to construct a
master plan for the planning horizon which is to be used as a basis for daily schedules to minimize
the cost and to maximize the coverage of customers as well as the similarity of daily routes with
the master plan.
Someofthe recentfindings in the literature onthe courierdeliveryproblemand on its variants
follow, mainly from the review in Zhong et al. (2007). For the deterministic pickup and delivery
problem systems, Hall (1996) studied the optimal route designs for overnight carriers using con-
tinuous space approximation models and demonstrated how these constraints affect the designs.
Haughton and Stenger (1998) modeled customer service for fixed route delivery systems under
stochastic demand. Later, Haughton (2000)developed a frameworkfor quantifying the benefits of
route re-optimization, again under stochastic customer demands. Recent research that considers
the design of large scale logistic systems for uncertain environments comes from Erera(2000). His
dissertationadaptsacontinuumapproximationmethodologydevelopedfordeterministicproblems
for analysis of large scale vehicle dispatching problems. The methodology is to find expected cost
approximations that allow near-optimal configuration of such systems under stochastic customer
locationsanddemand. ArecentworkbyZhongetal. (2007)proposesthe conceptsof“cell”, “core
area”, and “flex-zone” to develop a two-stage model which uses the capacity that is not assigned
in the first stage to adapt to the demand uncertainty in the second stage. They also improve
familiarity of drivers with the daily routes over the planning horizon by increasing the similarity
of each day’s schedule.
When the customer locations and demand are revealed at the beginning of each day, a simple
solution of routing for a planning horizon is to repeat the routing solution procedure each day
based on that day’s requirements, which will yield variable routes as opposed to fixed routes
which remain constantand which are usually correctedby a recourseprocedureto deal with route
failures. SavelsberghandGoetschalckx(1995)studiedtheviabilityoffixedroutesasanalternative
to variable routes. In their two-phase method, a generalized assignment heuristic constructs an
initial solution based on mean customer demand during the first phase; and in the second phase,
21
a local search procedure improves the current solution by minimizing the objective function of a
stochastic program with recourse. Beasley and Christofides (1997) studied the delivery operation
ofalargemailordercompany,whichusedasetoffixeddeliveryroutesineachday. Byaggregating
the customers according to postcode, the effective number of customers was reduced significantly.
22
Chapter 3
Robustness of VRP Formulations
3.1 Methodology
The first challenge in formulating a RVRP is to identify a suitable VRP formulation to derive
a tractable RVRP from uncertainty sets which are closely related to the application, and which
can yield a tractable robust counterpart. There exists a number of different VRP formulations
and since the RVRP formulation is related to the VRP formulation from which it is derived, it
is important to identify a VRP formulation which will yield a RVRP amenable for computation.
In this chapter, we identify the deterministic VRP formulation which will be used to derive the
RVRP formulation. The intuitive motivation is the following: if the VRP formulation is not
“inherently” difficult to solve, then the robust counterpart of it should not be more difficult to
solve. Therefore, in the search of a computationally amenable RVRP, it is important to identify
first a computationally amenable VRP.
Whendealingwiththecomplexityandpracticaldifficultyofsolvingaproblem,threeaspectsof
the problem become relevant: the problem size, how the problem is formulated, and the algorithm
used to solve the problem. In fact, classic complexity theory studies the asymptotic behavior of
algorithms as the problem size increases, and on the practical side, a number of algorithms are
developed to address specific formulations of the problem. In particular, the VRP is known to be
NP-hard (Lenstra and Rinnooy Kan 1981) and it is well-known that highly constrained VRPs,
23
such as those with tight capacity or time window constraints, are more amenable to a column
generation solution approach since the constraints lead to pruning up front a large number of
infeasible paths (Lu and Dessouky 2004). Therefore, there is clearly a strong relationship between
specific algorithms and problem formulation, and different problem formulations lead to different
complexities. We refer to this as the “inherent” difficulty of the problem formulation. How
the formulation of the problem influences this intrinsic difficulty of the problem and hence the
solution run time of any algorithm for VRPs of a fixed size in practice is not fully explored
in the literature. We are interested in identifying a VRP formulation which is computationally
robust with respect to the severalproblem parameterswhen the solution procedure is fixed. Since
our goal here is to discover the inherent difficulty of the VRP, this chapter concentrates on the
most fundamental version of VRP, the capacitated VRP. All the other variants of VRP which are
obtainedfromCVRPbyaugmentingtheformulation(byaddingside-constraintsand/oradditional
variables) are not fundamentally easier to solve than CVRP since they “inherit” the complexity
of CVRP while augmenting the formulation. Therefore, this chapter is dedicated to explore a
“suitable”’ formulation of CVRP which is to be later used in a robust optimization framework
when introducing uncertainty in problem parameters.
We aimto study experimentallythe effectofproblemformulationon solutionrun times forthe
CVRP when solved optimally by a commercial IP solver. The goal of this work is then similar to
Jones et al. (1993), which investigates the effect of problem formulation of multicommodity flow
problems on decomposition algorithms. However, by considering a general purpose solver we aim
to identify problem features that can indicate whether a problem formulation is fundamentally
easier to solve or not. Therefore, special solution procedures favoring particular formulations are
avoided. We note that there is a vast literature of experimental analysis on well-known hard
VRP instances combining formulation and solution procedure. Since such studies develop special
solution procedures for a given formulation, they do not identify the sole effect of the problem
formulation independent from its special solution procedure.
24
The effect of problem formulation is investigated based on a priori performance measures.
We are interested in identifying measures based on problem data to describe the difficulty of the
problem instance due to problem formulation, without solving the problem. There are other per-
formancemeasuresproposedintheVRPliteraturetoinvestigatetheeffectofproblemformulation.
However, these measures are mainly posterior and typically depend on the optimal solution. For
example the averagenumber of customers per route and the number of arcsselected in all optimal
solutions also inform how difficult it is to solve problems. Here, we experimentally investigate the
relation between problem data and problem formulation on the difficulty faced by a general IP
solver in solving a CVRP.
We denote by performance measures the measures of a problem that provide insight into the
difficulty of solving a problem. Therefore, problem size is by definition the essential performance
measure for the worst case analysis of algorithms for combinatorial problems. However, we are
interested in understanding how the problem formulation impacts solution times. Hence, we are
interested in data based performance measures not related to problem size. We also focus on
the observed solution run times and not on the worst case complexity of the algorithm. There
exists research on performance measures for NP-hard optimization problems that aim to explain
the observed variation in solution times. It is shown that these problems tend to exhibit sharp
changes in computational complexity, with instances going from being easy to being hard to solve
(or vice-versa) as a certain performance measure changes. These sharp frontiers are referred to
as phase transitions. In particular the work on the traveling salesman problem (TSP) (Gent and
Walsh 1996; Cheeseman et al. 1991; Zhang 2004)is related to performance measures for the VRP
due to the similarities between the problems. Phase transitions for the TSP have been identified
based on asymptotic properties of the optimal tour length as the problem size increases and the
existence of solution backbones, i.e. variables that are present in all optimal solutions (Gent and
Walsh 1996; Zhang 2004). Our approachis closer to Cheesemanet al. (1991)and Zhangand Korf
(1996), which identify phase transitions that depend on properties of the distance matrix of the
TSP, such as the range and standard deviation and not on asymptotic properties of the problem.
25
We consider the capacitated asymmetric vehicle routing problem with unit demand at each
demand location. Our study concentrates on three different arc-based formulations of the VRP
and uses CPLEX 8.1 as our state-of-the-art commercial IP solver to obtain exact solutions. We
limit our study to arc-based formulations because path-based formulations of the VRP have an
additional difficulty due to the variable explosion and thus require specialized column generation
algorithms. The experiments conducted explore whether vehicle capacity, number of vehicles,
and properties of the matrix of travel times, such as mean travel time or number of triangular
inequalities violated, as well as geographical distribution of the customers influence solution run
times. We present plots comparing the effect on the run times of the different formulations by
modifying each of these problem parameters.
In the next section we discuss different formulations of the VRP and introduce the problems
consideredin this chapter. Then we presentour experimental results in the following section. The
last section finishes with summary, conclusions, and general remarks.
3.2 VRP Formulations
The CVRP problem considers the problem of routing at minimum cost a uniform fleet of K vehi-
cles,eachwithcapacityC,toservicegeographicallydispersedcustomers,eachwithadeterministic
unit demand that must be serviced by a single vehicle. Let V be the set of n demand nodes and
a single depot, denoted as node 0. Let d
i
be the demand at each node i. Since we assume unit
demand, then d
i
= 1 for all i∈ V \{0}. We consider the fully connected network, and denote the
deterministic travel time between node i and node j by c
ij
. We assume asymmetric deterministic
travel times (i.e. c
ij
6=c
ji
in general) which might in addition violate the triangular inequality.
We now describe different formulations for the CVRP. Our presentation here also follows the
review in the book Toth and Vigo (2002b). In broad terms, formulations for the CVRP can be
separated into two main categories: path-based formulations and arc-based formulations. Path-
based formulations, alsoreferred to as set-partitioningformulations, considera binary variable for
26
each path (or vehicle route) which indicates whether that route is part of the routing solution or
not. These formulations therefore consider a number of integer variables that grow exponentially
with n and are typically tackled via column generation schemes. Our work does not consider
path-based formulations because of its dependence on specialized algorithms. We present now
different arc-based formulations, which we group into polynomial vehicle flow models, commodity
flow models, and exponential vehicle flow models.
3.2.1 Polynomial Vehicle Flow Models
Arc-based models consider integer variablesx
ij
which indicate whether a vehicle goes from node i
tonodej ornot. Problem(1)belowincludesadditionalcontinuousvariablesu
i
foreveryi∈V\{0}
that represent the flow in the vehicle after it visits customer i.
(VRP1) min
X
i∈V
X
j∈V
c
ij
x
ij
(1.1)
s.t.
X
i∈V
x
ij
=1 j ∈V \{0} (1.2)
X
j∈V
x
ij
= 1 i∈ V \{0} (1.3)
X
i∈V
x
i0
=K (1.4)
X
j∈V
x
0j
=K (1.5)
u
i
−u
j
+Cx
ij
≤C−d
j
i,j ∈ V \{0},i6=j (1.6)
d
i
≤u
i
≤C i∈ V \{0} (1.7)
x
ij
∈{0,1} i,j ∈ V (1.8) .
(1)
Problem VRP1 minimizes the total travel time while ensuring a feasible route. The first two
constraints, (1.2)and (1.3), force that a single vehicle goes into and out of every node; constraints
(1.4) and (1.5) ensure that a total of K vehicles leave the depot and then return to it; constraints
(1.6)and (1.7)aresubtour elimination constraintsimposing both the capacity and connectivity of
the feasible routes. They are also known as Miller-Tucker-Zemlin (MTZ) constraints. Note that
27
isolatedsubtoursviolatethe subtour eliminationconstraintsand thatinourcasetheseconstraints
simplify slightlysince the demand ateachnode isd
i
= 1. Although this formulationis polynomial
in size, it is well known that it can lead to weak LP relaxations. In fact, Desrochers and Laporte
(1991) propose tightening constraints for VRP1.
(VRP2) min
X
i∈V
X
j∈V
c
ij
K
X
k=1
x
ijk
(2.1)
s.t.
K
X
k=1
y
ik
= 1 i∈V \{0} (2.2)
K
X
k=1
y
0k
=K (2.3)
X
j∈V
x
ijk
=
X
j∈V
x
jik
=y
ik
i∈V, k ∈1...K (2.4)
X
i∈V
d
i
y
ik
≤C k ∈1...K (2.5)
u
ik
−u
jk
+Cx
ijk
≤C−d
j
i,j ∈ V \{0}, i6=j,
k ∈1...K
(2.6)
d
i
≤u
ik
≤C i∈V \{0}, k ∈ 1...K (2.7)
x
ijk
∈{0,1} i,j ∈V, k ∈1...K (2.8)
y
ik
∈{0,1} i,∈V, k ∈1...K (2.9) .
(2)
Our second vehicle flow model considered is known as the three-index formulation. This
formulation is derived from VRP1 by adding an index k to the integer variables x
ijk
to identify
which vehicle k traverses arc (i,j). In addition this formulation considers an integer variable y
ik
which indicates if vehicle k services customer i or not. Hence the three-index formulation we
consider is VRP2.
Constraints (2.2) and (2.3) ensure that exactly one vehicle services every customer i∈ V \{0}
and that K vehicles leave the depot. Constraints (2.4) force vehicle k to arrive and leave from
node i only if it services that node. Constraints (2.5) enforce the capacity constraint for every
vehicle k, while constraints (2.6) and (2.7) are the subtour elimination constraints.
28
3.2.2 Commodity Flow Models
Commodity flow models were recently introduced for the symmetric CVRP in the paper by Bal-
dacci et al. (2004). Below we describe an asymmetric version of this model. Commodity flow
models consider, in addition to the routing variables x
ij
, two sets of continuous non-negative flow
variables y
ij
and z
ij
for each arc (i,j). If a vehicle travels from i to j in the network, it carries a
load y
ij
and a residual capacityz
ji
so that y
ij
+z
ji
=C. The model also requires a dummy depot
node n+1 to the network. Let V
0
= V ∪{n+1} and also let D be the set of demand nodes, so
D∪{0,n+1}=V
0
. The vehicles in the model now depart node 0 and finish at n+1.
(VRP3) min
X
i∈V
0
X
j∈V
0
c
ij
x
ij
(3.1)
s.t.
X
j∈V
0
(y
ji
−y
ij
)=d
i
i∈D (3.2)
X
j∈V
0
(z
ji
−z
ij
)=d
i
i∈D (3.3)
X
j∈D
y
0j
=d(D) (3.4)
X
j∈D
z
j0
=KC−d(D) (3.5)
X
j∈D
z
n+1j
=KC (3.6)
y
ij
+z
ji
=Cx
ij
i,j ∈V
0
(3.7)
X
j∈V
0
(x
ij
+x
ji
) =2 i∈D (3.8)
y
ij
≥ 0 i,j ∈V
0
(3.9)
z
ij
≥ 0 i,j ∈V
0
(3.10)
x
ij
∈{0,1} i,j ∈V
0
(3.11) .
(3)
Constraints (3.2) and (3.3) enforce that the flow and residual capacity variables are adjusted
according to the demand at every node i∈D; constraints (3.4) send the total demand out of the
depot; constraints (3.5) and (3.6) respectively enforce total amounts of residual capacity into the
depot and out of the dummy node; constraints (3.7) ensure the flow and residual capacity add
29
up to C for every arc that is selected, and set these variables to zero on arcs that are not used;
finally constraints (3.8) are the routing constraints enforcing that exactly one vehicle arrives to
and departs from every node i.
3.2.3 Exponential Vehicle Flow Models
Different formulations of the VRP are obtained by replacing the subtour elimination constraints
(1.6) and (1.7) in problem (1) with either the capacity-cut constraints (1.9) or the generalized
subtour elimination constraints (1.10) below:
X
i6∈S
X
j∈S
x
ij
≥γ(S) S ⊆V \{0}, S 6=∅ (1.9) ,
X
i∈S
X
j∈S
x
ij
≤|S|−γ(S) S ⊆V \{0}, S 6=∅ (1.10) .
Constraint (1.9) enforces that the number of vehicles entering any subset S must be enough to
satisfy the demand at S. Constraint (1.10) implies that the number of vehicles that leave S is
enough to satisfy the demand at S. Here the function γ(·) is defined for every S ⊆ V, such that
γ(S) equals the optimal value for the bin packing problem of demand d(S) =
P
i∈S
d
i
(= |S|)
with bins of size C. The function γ(S) can be equivalently replaced by the bin packing problem
lower bound given by dd(S)/Ce. In fact, in the case of unit demand this lower bound is tight.
Duetotheexponentialexplosioninthenumberofconstraints,theabovemodelsbecomesignif-
icantly large as the problem size increases. Therefore in practice, they are not solved completely,
but they are used in a constraint generation procedure. Since we are interested in this study on
the performance of formulations on a general IP solver, we do not present the results for these
exponential models. However, it is known that these models are mostly insensitive to problem
parameters since they result in tight LP relaxations.
30
3.3 Experimental Results
In this experimental section we first identify, for each formulation, the problem size for which an
exact solution can be obtained in a reasonable amount of time to conduct a parameter sensitivity
analysis. We then explore the effect of vehicle capacity C and number of vehicles K on problem
solutiontime. Weinvestigatehowdifferentaspectsofthe traveltime matrixc
ij
affectthe problem
solution time. Our last computational experiment explores the effect of geographical distribution
of customers on the solution time.
For each problem formulation, and each setting of problem parameters (demand nodes n,
vehicle capacity C, number of vehicles K, and mean μ
c
and standard deviation σ
c
of the travel
time matrix c
ij
) we create random instances by randomly generating the n(n+ 1) off-diagonal
values of the travel time matrix following a lognormal distribution with mean μ
c
and standard
deviation σ
c
. The lognormal distribution assumption for actual travel time distances is standard
in the transportation literature, see Dessouky et al. (1999). Our computational analysis considers
a base case scenario with the parameter values C =6, K = 5, μ
c
= 4.5, and σ
c
=3.375. For each
scenario, 50 instances are randomly generated and here we define the problem solution time as
the average solution time of these instances.
All experiments were carried out on a Dell Dimension 8200 computer with 2.4 GHz Pentium
4 Processor and 1GB RDRAM running Red Hat Linux 7.3. Exact solutions for the three VRP
formulations considered were obtained using the mixed integer programming solver from CPLEX
8.1 with default settings.
3.3.1 Problem Size
Our goal in selecting a problem size is to obtain runtimes that allow a study of the sensitivity of
the solution time with respect to changes in problem parameters. Therefore, we select different
problemsizesforeachformulation. AftersomepreliminarycomputationsweobservedthatVRP2,
the three-indexformulation,ismoredifficulttosolvethan the otherpolynomialformulationssince
31
VRP2 has more than K times the number of integer variables than either VRP1 or VRP2. In
the remainder of the paper we use n = 26 for VRP1 and VRP3, and n = 20 for VRP2. With
these sizes, the mean solution time with base case parameters are between 1 to 2 minutes for all
formulations. These solution times are small enough so that the experiments can be executed in
a reasonable time and also significant enough so that changes to solution time can be noticed.
3.3.2 Effect ofC and K on Solution Time
We plot in the left side of Figure 3.1 the mean solution time for each of the 3 formulations as
we vary the vehicle capacity C. Note that in these experiments, μ
c
and σ
c
are kept at their base
case values of 4.5 and 3.375, respectively and we fix the number of vehicles K = 5. We note that
as the vehicle capacity increases all three formulations become easier to solve essentially because
the problem is becoming less constrained and optimal solutions are simpler. In particular, if the
vehicle capacity increases, then it is more likely that efficient routing solutions will have enough
capacity to satisfy the demand. We obtained very similar results when we fix the vehicle capacity
C and vary the number of vehicles K because as the number of vehicles K increases, each vehicle
can perform shorter routes in the optimal solution.
IntherightsideofFigure3.1,wealsoplotinformationaboutthedistributionofthe50instances
for each data point. For all three formulations, around each median value, we plot vertical lines
between the first and the third quartiles of 50 instances. Despite some fluctuations, we believe
the averages are good enough representations of the underlying distribution of the 50 instances.
We observed similar trends in all of the experimental results and therefore we only present the
averages in the reminder of this chapter to make the figures clearer.
3.3.3 Effect of the Travel Time Matrix on Solution Time
To investigate the effect of the travel time matrix on solution time we present in Figures 3.2, 3.3,
and 3.4 the mean solution times as we respectively vary only the mean μ
c
, standard deviation σ
c
,
32
5 6 7 8 9 10
0
10
20
30
40
50
60
Vehicle Capacity
Average CPU Time (Sec)
VRP1 n26
VRP2 n20
VRP3 n26
5 6 7 8 9 10
0
10
20
30
40
50
60
Vehicle Capacity
Distribution of CPU Time (Sec)
VRP1 n26
VRP2 n20
VRP3 n26
Figure 3.1: LEFT: Average CPU time as a function of the vehicle capacity. RIGHT: Median, 1
st
and 3
rd
quartiles. With μ
c
= 4.5, σ
c
= 3.375, K = 5, 50 random instances
and coefficient of variation CV
c
=
σ
c
μ
c
, maintaining C = 6 and K = 5. Since the coefficient of
variation is defined in terms of both μ
c
and σ
c
, there are many combinations that yield the range
explored in the experiment. In our settings we consider:
μ
c
2.5 3.5 4.5 5.5 6.5 7.5 8.5 9.5 10.5 11.5
σ
c
0.625 1.75 3.375 5.5 9.75 15 25.5 38 52.5 69
CV
c
0.25 0.5 0.75 1 1.5 2 3 4 5 6
For VRP1 and VRP3 the mean does not seem to impact solution times. However, there is a
largedependencyonσ
c
,whichshowsanincreaseofmorethan25secondsandthenadecrease. The
resultswithrespecttoCV
c
arelessclear,asitseemstoshowastrongincreaseforsmallcoefficients
andtobealmostindifferenttolargercoefficients. FormulationVRP2exhibitsadifferentbehavior.
It decreases significantly with an increase in μ
c
and it increases significantly with an increase in
σ
c
. However for CV
c
, the effect is mixed.
33
30
40
50
60
Average CPU Time (Sec)
VRP1 n26
VRP3 n26
2 2.5 3.5 4.5 5.5 6.5 7.5 8
80
90
100
110
Mean, μ
c
Average CPU Time (Sec)
VRP2 n20
Figure 3.2: Average CPU time as a function of mean travel time μ
c
. With C = 6, K = 5,
σ
c
= 3.375, 50 random instances
30
40
50
60
70
80
Average CPU Time (Sec)
VRP1 n26
VRP3 n26
0.625 1.75 3.375 5.5 9.75 15 16
50
100
150
200
250
Standard Deviation, σ
c
Average CPU Time (Sec)
VRP2 n20
Figure 3.3: Average CPU time as a function of standard deviation of travel time σ
c
. With C =6,
K =5, μ
c
= 4.5, 50 random instances
34
30
40
50
60
70
Average CPU Time (Sec)
VRP1 n26
VRP3 n26
0.25 0.75 1.5 2 3 4 5 6 7
50
100
150
200
250
Coefficient of Variation
Average CPU Time (Sec)
VRP2 n20
Figure 3.4: Average CPU time as a function of the coefficient of variation of travel time σ
c
/μ
c
.
With C =6, K = 5, 50 random instances
3.3.4 Understanding the Effect of CV
c
on Solution Time
Taking a closer look at the results with respect to varying CV
c
in Figure 3.4, we can explain the
split behavior in VRP1, VRP2, and VRP3 (first increase followed by relative independence from
CV
c
)bytheinterplayoftwocompetingforces: anincreaseinthenumberoftriangularinequalities
violated and an increase in the asymmetry of the travel time matrix. Before discussing these
two competing forces in detail, we comment on how the random problems are generated for the
experiments.
For each random instance we generate n(n + 1) random lognormal numbers with mean μ
c
and standard deviation σ
c
. This instance is used for all three formulations. However since each
formulationconsidersa different number ofnodes, each formulationuses as manyof the generated
random numbers as appropriate. We generate 26∗27random lognormalnumbers, all of which are
used for VRP1 and VRP3, while VRP2 uses the first 20∗21 numbers.
The fact that the number of triangular inequalities violated and the asymmetry of the travel
time matrix increase as CV
c
increases is not surprising. For instance, for a fixed mean, CV
c
increasesby increasingthe standard deviation. Ifthe traveltime matrix is generatedwith a larger
35
variability then it is reasonableto have more triangularinequality violations and a largerdistance
from symmetric instances. We also show this relationship experimentally in Figure 3.5. The
number of triangular inequalities violated is obtained by checking all inequalities for every subset
of three arcs in the matrix. We define the asymmetry of a square matrix A by
1
√
2
kA−A
T
k
F
=
s
X
i>j
(A
ij
−A
ji
)
2
.
0.25 0.75 1.5 2 3 4 5 6 7
0
2000
4000
6000
Triangular Inequalities Violated
VRP1 n26
VRP2 n20
VRP3 n26
0.25 0.75 1.5 2 3 4 5 6 7
0
500
1000
1500
Coefficient of Variation
Asymmetry
VRP1 n26
VRP2 n20
VRP3 n26
Figure 3.5: Number of triangular inequalities violated and asymmetry as a function of the coeffi-
cient of variation CV
c
. With C =6, K = 5, 50 random instances
We note that the number of triangular inequalities violated increases sharply for small values
of CV
c
and slowerfor largeCV
c
, while the rate of increase of the asymmetry is very linear. Hence
for small CV
c
the number of triangular inequalities violations will be the dominating factor due
to its sharp increase, whereas for greater values of CV
c
the asymmetry of the matrix becomes
more important. Note also that since VRP1 and VRP3 use the same data, both the number of
triangular inequalities violated and the asymmetry are the same for those formulations.
The fact that the number of violated triangular inequalities is positively correlated with so-
lution times is suggested from the existence of better worst case complexity results for the TSP
when the distances satisfy triangularinequalities (Christofides 1976; Orponen and Mannila 1987).
In particular, we observe experimentally that the number of triangular inequality violations is
36
positively correlated with the initial LP gap, which is of course related to the solution time. Both
these trends appear in Figure 3.6 for VRP1 and VRP3, and also hold for the other formulation
VRP2. Notethattheseplotsconsiderrangesforthevaluesonthex-axis,asitisextremelydifficult
to construct random examples with a specific number of violated triangular inequalities or with
a given initial LP gap. Therefore, in these plots, the y values are the mean of a certain number
of problems that fall within each range, or bucket. The actual number of problems used varies
for each formulation since it is determined by the bucket with fewest problems. For each problem
formulation we indicate the number of problems used to compute the mean in the graph legend,
next to the formulation name and number of demand points.
[0−1000) [1000−2000) [2000−3000) [3000−4000) [4000−5000+)
0
2
4
6
8
10
Triangular Inequalities Violated
% LP Gap
VRP1 n26 b45
VRP3 n26 b30
[0−2) [2−4) [4−6) [6−8) [8−10+)
20
40
60
80
100
120
140
160
% LP Gap
Average CPU Time (Sec)
VRP1 n26 b50
VRP3 n26 b50
Figure 3.6: Initial LP gap as a function of triangular inequalities violated and solution time as a
function of initial LP gap for VRP1 and VRP3. With C = 6, K =5
In conclusion, for small CV
c
, the main effect of increasing CV
c
is an increase in the number
of violated triangular inequalities, which implies an increase in solution time, a phenomenon that
is observed in Figure 3.4 for small CV
c
.
For large CV
c
we expect the dominating factor to be the asymmetry of the travel time matrix.
We observe the effect of this factor on solution time in Figure 3.7, and note that for VRP1 it
37
becomes easier to solve as asymmetry increases and for VRP2 and VRP3 the effect is mixed
showing an increase and then a decrease.
[0−250) [250−500) [500−750) [750−1000) [1000−1250+)
30
35
40
45
50
55
60
65
70
75
Asymmetry
Average CPU Time (Sec)
VRP1 n26 b50
VRP3 n26 b50
[0−125) [125−250) [250−375) [375−500) [500−625+)
50
100
150
200
250
300
350
Asymmetry
Average CPU Time (Sec)
VRP2 n20 b50
Figure 3.7: AverageCPUtime as a function of the asymmetryin the matrix oftraveltimes. With
C = 6, K =5
Finally, to validate the fact that there are two competing factors influencing solution time as
we increase CV
c
, we plot in Figure 3.8 how the solution times are affected by an increase in CV
c
for symmetric data instances. We note that for formulation VRP1, and to some extent VRP2, the
solutiontimes increaseasthecoefficientofvariationincreases. Thereforethe factthattheincrease
stops for large CV
c
in Figure 3.4 is most likely due to the asymmetry. A second observation from
Figure 3.8 is that instances with symmetric data are harder to solve, forcing a reduction in the
size of problems considered for formulations VRP2 (from 20 to 13), and VRP1 and VRP3 (from
26 to 18). We note that these symmetric instances were solved using the asymmetric formulations
presentedearlier,hencethedifficultystemsfromhavingmanyrouteswithsimilartotalcosts,since
now traversing a route clockwise or counterclockwise are equivalent.
38
0
200
400
600
Average CPU Time (Sec)
VRP1 n18
0.25 0.75 1.5 2 3 4 5 6 7
0
20
40
60
Coefficient of Variation
Average CPU Time (Sec)
VRP2 n13
VRP3 n18
Figure 3.8: Average CPU time as a function of the coefficient of variation CV
c
for symmetric
instances. With C = 6, K =5, 50 random instances
3.3.5 Effect of Geographical Distribution of Customers
Lastly, we investigate the effect of geographical distribution of customers on the solution time.
For this, we randomly generate customers lying on a 2-dimensional Euclidean plane and then
construct the cost matrix. Problems generated with this procedure will have symmetric cost
matrix which also obeys the triangular inequalities. We randomly generate customers on the
plane with different degrees of clustering. For that, we consider a depot at location (0,0) and
randomly generate customer locations in 5 different clusters (one for each vehicle) with identical
radius of r = 20. We center each cluster at a random location a distance R from the depot. By
increasing the value of R, we change the geographical distribution of customers from uniformly
distributed to clustered. We measure how much the customers are clustered with the following
clustering ratio: r
c
=
R
r
. Figure 3.9 shows averages of 50 random instances for each clustering
ratio. Clearly, VRP3 benefits the most from the effect of clustering while VRP1 benefits less. For
those two formulations clustering helps decomposing the problem naturally and makes it easier to
solvesince eachvehicleisassignedtoa singlecluster. Onthe otherhand, this hasanegativeeffect
on VRP2. Since VRP2 distinguishes between the vehicles, there are several different assignments
39
of vehicles to the clusters resulting in the same total cost and not being able to prune those
solutions increases the solution time of VRP2.
0 2 4 6 8
0
50
100
150
200
250
Clustering Ratio
Average CPU Time (Sec)
VRP1 n18
VRP2 n13
VRP3 n18
Figure 3.9: Average CPU time as a function of the clustering ratio r
c
for symmetric instances.
With C =6, K = 5, 50 random instances
3.4 Summary
Inthis chapterweinvestigatetheeffectofproblemformulationfordifferentcapacitatedVRPsand
of data parameters on the solution run time of a state-of-the-art general purpose IP solver. By
concentrating on the performance of a general IP solver we aim to represent an intrinsic difficulty
oftheproblemformulationandintheprocessidentifythe bestproblemformulationsforasolution
method that has become more relevant both for exact and heuristic solutions.
We find that the initial integrality gap of polynomial arc-based formulations makes the perfor-
mance sensitive to the problem parameters. We now summarize our findings of the experiments
presented.
• All models show a significant decrease in solution time if we increase the capacity or the
number of vehicles.
40
• The effect of the travel time matrix on solution time is more obscure. VRP2 shows a strong
dependency on the standard deviation, while VRP1 and VRP3 show first an increase in
solution time, then oscillatory behavior as CV
c
increases.
• This mixed dependency on CV
c
is explained by two factors. As CV
c
increases we have that
the number of violated triangular inequalities increases and the asymmetry of the matrix
increases. The increase in triangular inequalities violations is most significant for smaller
values of CV
c
and leads to an increase in solution time, while the asymmetry increase
becomes the dominating factor for large CV
c
.
• Problem VRP1 becomes easier to solve as the asymmetry increases, while problems VRP2
and VRP3, first show an increase then a decrease in solutions times with asymmetry.
• Forsymmetricproblems,anincreaseincoefficientofvariationleadstoanincreaseinsolution
times.
• For symmetric problems, when the degree of clustering of the geographical distribution
of customers increases, VRP1 and VRP3 become easier to solve whereas VRP2 shows an
increase in solution times.
We observed similar performance trends when considering symmetric problems, and have no
reason to suspect a deviation from this behavior for the symmetric versions of the formulations
considered. However,specialized algorithmsthat favora specific formulation, can lead to different
results. The question of which formulation has a specialized algorithm that performs best is at
the heart of a number of VRP challenges and competitions, focuses on specific problem instances
and is not investigated here. For instance, using a path-based formulation could be a competitive
solution approach for problems with hard time windows, which helps reducing the number of
feasible solutions.
41
Chapter 4
A Priori Robust Optimization for VRP
4.1 Methodology
The first challenge in developing RVRP is to identify a suitable VRP. Based on the findings in
the previous chapter, we identify in this chapter this suitable VRP formulation and address the
challengeofdevelopingrobustcounterpartsforthecaseofuncertaintyindemandandtraveltimes.
We first develop RVRP with demand uncertainty and gradually add uncertainty in travel times.
The proposed RVRP in this chapter identifies an optimal “a priori” route that is feasible for
every realization of the uncertainty and does not consider recourse actions. In the case of demand
uncertainty, such solutions are possible if the uncertainty is small or the vehicle fleet has sufficient
capacity. A different model is needed if there is no route that can meet all demand realizations
and we saythat the RVRP is infeasible. In problems with recoursethere is an additional difficulty
of determining optimal recourse actions. This leads either to augmenting the original determin-
istic model or requiring strong uncertainty assumptions to evaluate the recourse action. In both
cases, the additional flexibility of the recourse solution comes with an increased complexity in
the formulation and solution procedure. In this chapter we propose a model that uses a simple
representation of the uncertainty and does not make the problem more difficult to solve. We con-
centrate on studying when this a priori RVRP can be beneficial, providing solutions that address
42
well the uncertainty without being overly conservative, and leave to the next chapter extensions
of this model to recourse actions.
A natural method to address demand uncertainty is to reserve vehicle capacity to be able to
adapt to cases when the realized demand is greater than the expected demand. In fact, if there
is abundant vehicle capacity, such as in the uncapacitated VRP, the optimal routing solution
can easily accommodate changes in the demand levels. However, in capacitated cases where the
vehicle capacity is slightly greater than the expected demand, the interesting problem is how to
manage the extra vehicle capacity to distribute slack among routes to better address the demand
uncertainty. The RVRP distributes this slack by finding a routing solution at minimum cost that
satisfies all possible demand realizations.
A few methods developed for the deterministic VRP have focused explicitly on how to dis-
tribute the vehicle capacity among routes. For instance, Daganzo (1988) proposes the use of a
consolidation center as a strategy to better manage vehicle capacity. Charikar et al. (2001) intro-
duces a constantratio approximationalgorithm which assigns each vehicle half of its capacity and
uses the remaining capacity to improve routes with a matching algorithm. Branke et al. (2005)
shows that managing the slack by waiting at strategic locations can increase the probability of
meeting additional demand. Our work departs from these prior results as we consider a different
problem domain: a standard CVRP with no transshipment nodes and with a small capacity to
demand ratio. For the stochastic VRP, Zhong et al. (2007) developed a two-stage model that
uses the capacity that is not assigned in the first stage to adapt to the demand uncertainty in the
second stage. The first stage creates “core areas” to be serviced and after the demand is realized
the recourseactions decide how to route these areas allowingfor exchanging demand nodes on the
“flex-zones” between core areas. They show that keeping unassigned customers near the depot is
a good strategy for balancing the workload due to daily demand variations. Although our work
develops an a priori routing strategy with no recourse, we also notice that having customers near
the depot facilitates the creation of efficient robust routes.
43
When it comes to the challenge of solving RVRP, we benefit from an open source solver as
an efficient exact solution procedure for VRP. We propose performance measures to contrast the
deterministic and the robust solution and provide experimental results for the case of demand
uncertainty for which we also explore the effect of distributing the excess vehicle capacityover the
network. These include a comparison of the robust and deterministic solution on a well-known
suite of CVRP problems (Augerat et al. 1995), a comparison on a family of clustered instances,
and verifying that the robust solution better addresses the demand uncertainty than a uniform
and a non-uniform distributions of unused vehicle capacity.
In the next section, we develop a general RVRP. We gradually add uncertainty in demand and
travel times in the following two sections. We present our computational results in Section 4.4
and finish this chapter with some conclusions in Section 4.5.
4.2 Robust Counterpart of VRP
In this section, we first identify the VRP formulation to be used in the derivation of the robust
counterpart and we also justify the choice of the formulation. Then we mathematically derive a
general RVRP formulation with uncertainty in both demand and travel times.
4.2.1 Identifying the VRP Formulation
Wehaveinvestigatedinthepreviouschapterthequestionofidentifyinga“suitable”VRPformula-
tionbyconductinganexperimentalcomparisonofdifferentCVRPformulations. Ourexperimental
analysis on a general IP solver revealed that the CVRP formulation with MTZ subtour elimina-
tion constraints which are polynomially many in the size of the number of nodes, outperformed
computationally the other edge-based VRP formulations.
Anotherimportantcriterioninidentifyinga“suitable”formulationforourrobustoptimization
framework is the nature of the formulation with respect to uncertain parameters. In this chapter,
we are interested in introducing uncertainty into both demand and travel times. In particular,
44
when we consider the parts of the formulation related to demand, the MTZ formulation has
constraints in the form of inequalities. In the robust optimization methodology, the inequalities
involving uncertain parameters are preferable to the equalities involving uncertain parameters.
Ben-Tal et al. (2003) shows the latter has a further increased complexity since an equality con-
straint can only be satisfied for multiple values by a problem with recourse which allows some
variables to be decided after the uncertainty is realized. Since the MTZ formulation results in
reasonable solution times and the nature of the formulation is suitable to introduce uncertainty
in the robust optimization framework, we adopt the MTZ formulation as our deterministic CVRP
formulation and use it to derive the robust counterparts.
On the other hand, although the nature of the MTZ formulation is the preferred one with
respect to uncertain parametersfor our robustoptimization framework, this formulation is known
to be very sensitive to the problem parameterswhich is most likely due to the largeinitial LP gap
oftheformulation. However,itispossibletoimprovetheinherentweakLPgapofthisformulation
by adopting the lifting techniques proposed by Desrochers and Laporte (1991).
The idea of lifting a constraint relies on introducing new variables into that constraint with
dummy coefficients and then analytically determining the maximum value these dummy coeffi-
cients can take. Meaningful choice of additional decision variables may yield better values for
the dummy coefficients to be fixed, which will make the constraint tighter. The resulting lifted
constraints can be added to the formulation as cuts to decrease the LP gap. In our case, this
lifting technique will further augmentthe CVRP formulation by adding the constraints(1.11)and
(1.12) which are the first and second order lifted versions of the MTZ constraints respectively; we
also need to replace constraint (1.7) by (1.13) and (1.14).
45
u
i
−u
j
+Cx
ij
+(C−d
i
−d
j
)x
ji
≤C−d
j
i,j ∈ V \{0} (1.11)
u
i
−u
j
+Cx
ij
+(C−d
i
−d
j
)x
ji
+d
k
(x
ik
+x
kj
)≤C−d
j
i,j,k∈ V \{0} (1.12)
d
i
+
X
j∈V, j6=i
x
ij
d
j
≤u
i
i∈ V \{0} (1.13)
u
i
≤C−(C−d
i
− max
j∈V, j6=i
d
j
)x
0i
−
X
j∈V, j6=i
x
ij
d
j
i∈ V \{0} (1.14)
23 24 25 26 27 28 29
0
50
100
150
200
250
300
Average CPU Time (Sec)
Number of Nodes
MTZ
Lifted MTZ
Figure 4.1: Effect of lifting on solution time
Theliftingmethodkeepsthetotalnumberofconstraintspolynomialontheorderofthenumber
of nodes, without further increasing the complexity of the CVRP formulation. The gain in the
solution times however is significant due to improved LP gap. Figure 4.1 shows solution times of
the CVRP with MTZ constraints versus the CVRP with lifted MTZ constraints when solved by a
generalIPsolver. Eachdatapointinthefigureisanaverageof30randominstances. However,the
46
effectoftheLPgapofthe formulationand the improvementofthe lifting technique aredependent
on the solution procedure used. The lifted MTZ constraints have not been fully explored in the
literature when it comes to developing a solution procedure for the VRP.
4.2.2 Deriving a General RVRP Formulation
In this subsection, we apply robust optimization methodology to obtain the robust counterpart
of the CVRP formulation identified in the previous subsection. Our goal is to obtain a compu-
tationally tractable RVRP with uncertainty in both demands and travel times. For the former,
we will slightly modify the formulation. Notice that the demand d
i
appears by itself and only on
constraints (1.6) and (1.7). However, the lower bound on constraint (1.7) is implied from (1.6),
the fact that every node is visited, and that u
i
≥ 0 for all i ∈ V \{0}. We will therefore only
consider the uncertainty in (1.6) and replace all d
i
with 0 in constraint (1.7).
For the purposes of our discussion here, we will consider that the vector of demands d and
matrixoftraveltimescareuncertainandbelongtosomeboundeduncertaintysets,i.e. d∈ U
D
and
c∈ U
C
respectively. Note that we do not make any assumptions on the probability distributions
of d and c other than they are bounded. The robust solution aims at obtaining the solution that
will perform well for all the uncertainty scenarios, and is obtained by minimizing the worst case.
Toillustratethe robustcounterpartproblem, orthe problemwhosesolutionisthe robustsolution,
consider the slightly modified version of the CVRP formulation in Equation 1 as explained. Let
P(c,d) denote the feasible set, and let
z(x,u;c,d)=
X
(i,j)∈A
c
ij
x
ij
+Υ
P(c,d)
(x,u) ,
whereΥ
S
(x,u)istheindicatorfunctionofsetS,thatisΥ
S
(x,u)= 0ifx,u∈S,andΥ
S
(x,u)=∞
otherwise. In this notation, the robust counterpart problem is given by the following mix integer
programming problem
z
∗
= min
x,u
max
c∈UC, d∈UD
z(x,u;c,d) .
47
As an alternative method, instead of minimizing the maximum occurrence of uncertainty,
the worst case, we can alter the problem in a way to minimize over all possible realizations of
uncertainty. The two methods are equivalent and result in the same solution. We will follow the
latter definition of robustness. Then the general form of the robust CVRP with uncertainty in
demand and travel times to minimize over all possible realizations of uncertainty is given below
by the RVRP formulation.
(RVRP) min γ
s.t.
X
i∈V
X
j∈V
c
ij
x
ij
≤γ ∀c∈U
C
u
j
−u
i
+C(1−x
ij
)≥d
j
i,j ∈V \{0}, i6=j ∀d∈ U
D
(1.2−1.5), (1.7−1.8) .
4.3 Uncertainty in Demand
In this section, we consider only the demand parameter d to be uncertain and therefore we fix
the definition of uncertainty for travel times as U
C
= {c} in the RVRP formulation to make it
deterministic. Recall that we consider the problem only with uncertainty in constraint (1.6) with
constraint (1.7) equal to 0 ≤ u ≤ C. We first propose the robust counterpart problem RVRP
for CVRP with demand belonging to an uncertainty set U
D
. Then we discuss the performance
measures to be used in the experiments together with some bounds on the measures.
4.3.1 Derivation of the RVRP
Weconsideruncertaintysetswhichareconstructedasdeviationsaroundanexpecteddemandvalue
d
0
. Thepossibledeviationdirectionsfromthesenominalvaluesarefixedandidentifiedbyscenario
vectors, d
k
, each of size n which is the number of nodes. The scenario vectors are allowed to have
negative deviation values. For a given number of scenario vectors, s, the general uncertainty set
U
D
is a linear combinationofthe scenariovectors. These weights, unknown combinationvariables
48
y
k
, then constitute a vector y of size s. The only restriction on the uncertainty set U
D
is that the
combination variables y
k
associated with each fixed scenario vector belongs to a bounded set Y.
U
D
={d|d
0
+
X
k
y
k
d
k
, k = 1,...,s, y ∈Y}
The set U
D
yields a tractable RVRP for different types of bounded sets Y. In particular, we
are interested in the following three sets for Y: convex hull, box, and ellipsoidal sets.
If set Y is a convex hull as in Y
1
, the resulting uncertainty set for demand is U
D1
.
Y
1
={y | y
k
≥ 0,
X
k
y
k
≤ 1, k =1,...,s}
If set Y is a box as in Y
2
, the resulting uncertainty set for demand is U
D2
.
Y
2
={y |−1≤y
k
≤ 1, k =1,...,s}
If set Y is an ellipsoid as in Y
3
, the resulting uncertainty set for demand is U
D3
.
Y
3
={y |
X
k
y
2
k
≤1, k = 1,...,s}
We now propose the robust counterparts, RVRPs for problems with demand uncertainty. By
replacing constraint (1.6) in CVRP with the constraint below (1.15) we obtain the RVRP with
respect to the general definition of the uncertainty set U
D
.
u
j
−u
i
+C(1−x
ij
)≥d
j
i,j ∈ V \{0}, i6=j ∀d∈U
D
(1.15)
49
If we substitute in the definition of the uncertainty set U
D
, we can write the robust constraint
(1.15) as the following inequality.
u
j
−u
i
+C(1−x
ij
)−d
0
j
≥
s
X
k=1
y
k
d
k
j
∀y ∈Y i,j ∈ V \{0}, i6=j (1.16)
For given decision variables x and u we refer to the left hand side of the above inequality as
φ
ij
(x,u) = u
j
−u
i
+C(1−x
ij
)−d
0
j
for i,j ∈ V \{0}, i 6= j. Then, to enforce that the above
inequalityholdsforally ∈ Y itsufficestoenforceitfor sup
y∈Y
s
X
k=1
y
k
d
k
j
= sup
y∈Y
y
T
D
j•
. Herewedenote
byD = [d
1
...d
s
]∈<
n×s
the matrix ofscenariovectorsand we denote byD
j•
= (d
1
j
,...,d
s
j
)
T
the
j-th rowof D as a column vector. Let us alsodenote e as the column vectorof all 1 of appropriate
dimension.
Proposition 1. Under uncertainty set U
D1
, the robust counterpart is obtained by replacing con-
straint(1.6) in CVRPwiththe constraintbelow (1.17). We refer tothe resultingRVRPas RVRP
1
.
u
j
−u
i
+C(1−x
ij
)≥d
0
j
+max{max
k
d
k
j
,0} i,j ∈ V \{0}, i6=j (1.17)
Proof: Using the definition of Y
1
we can write sup
y∈Y1
y
T
D
j•
and its dual as the following pair of
LPs:
(Primal) max y
T
D
j•
(Dual) min θ
s.t. e
T
y ≤ 1 s.t. θe≥D
j•
y ≥ 0 θ ≥0 .
From weak duality, the condition φ
ij
(x,u) ≥ sup
y∈Y1
y
T
D
j•
is equivalent to φ
ij
(x,u) ≥ θ for some
dual feasible θ. This means that φ
ij
(x,u)≥0 and φ
ij
(x,u)≥d
k
j
for k = 1,...,s. Combining these
conditions together for all φ
ij
(x,u) gives (1.17).
50
Proposition 2. Under uncertainty set U
D2
, the robust counterpart is obtained by replacing con-
straint(1.6) in CVRPwiththe constraintbelow (1.18). We refer tothe resultingRVRPas RVRP
2
.
u
j
−u
i
+C(1−x
ij
)≥d
0
j
+
X
k
|d
k
j
| i,j ∈ V \{0}, i6=j (1.18)
Proof: Using the definition of Y
2
we can write sup
y∈Y2
y
T
D
j•
and its dual as the following pair of
LPs:
(Primal) max y
T
D
j•
(Dual) min e
T
(α+β)
s.t. y ≤e s.t. α−β =D
j•
y ≥−e α,β ≥ 0 .
It is simple to verify that the optimal solution to the dual problem will satisfy the following
for every k = 1,...,s: α
∗
k
+β
∗
k
= |d
k
j
|. Therefore the dual optimal objective value is
s
X
k=1
|d
k
j
|.
Enforcing the robust feasibility condition onφ
ij
(x,u) with the above optimal dual objective value
we obtain (1.19).
Proposition 3. Under uncertainty set U
D3
, the robust counterpart is obtained by replacing con-
straint(1.6) in CVRPwiththe constraintbelow (1.19). We refer tothe resultingRVRPas RVRP
3
.
u
j
−u
i
+C(1−x
ij
)≥d
0
j
+
q
D
T
j•
D
j•
i,j ∈ V \{0}, i6=j (1.19)
Proof: Using the definition of Y
3
we have that sup
y∈Y3
y
T
D
j•
= maxy
T
D
j•
: y
T
y ≤ 1. From the
KKToptimalityconditionswehavethattheoptimalsolutiontothisproblemisy
∗
=
1
√
D
T
j•
Dj•
D
j•
.
When we plug this optimal solution into the robust feasibility condition φ
ij
(x,u)≥(y
∗
)
T
D
j•
, we
obtain (1.19).
For the three RVRPs with demand uncertainty studied, the only change from the original
CVRPformulationisanincreaseinthe demandsthatappearin(1.6). Since the deviationvectors,
51
d
k
, are fixed, each of the RVRPs is an instance of the CVRP. We can therefore make use of the
efficient exact algorithmsin the literature to solve the robustproblems. We note that the demand
parameters used in the RVRPs are at least as big as the deterministic demand parameters, thus
the RVRPsaretypicallymorecapacityconstrainedthan thecorrespondingdeterministic problem.
We note this because it has been observed in practice that solving CVRP becomes harder as
the problem is more capacity constrained. Thus, although the RVRPs are instances of CVRP,
in practice solving RVRPs is likely to be more difficult than solving the deterministic versions.
Lastly note that, depending on the nature of the scenario vectors, RVRPs may result in infeasible
problems even though the deterministic CVRP is feasible.
For different types of demand uncertainty sets d ∈ U
D
, the key step in the derivation of the
RVRP is to compute sup
d∈UD
d
j
and substitute this value into the right hand side of equation (1.15).
This can be done for different uncertainty sets than considered here. We do not pursue these
formulations here for simplicity, since many require additional constraints and variables making
the resulting robust problem not a CVRP.
4.3.2 Performance Measures and Bounds
We now propose two performance measures to contrast the robust and deterministic CVRP. The
first performance measure, the ratio r
rd
, quantifies the relative extra cost of the robust solution
with respect to the cost of the deterministic solution. It is given by r
rd
=
z
r
−z
d
z
d
where z
d
is
the optimum solution of the deterministic CVRP and z
r
is the optimum solution of the robust
counterpart. This ratio gives information on how much extra cost we will incur if we want to
implement the robust solution to protect against the worst case realization of the uncertainty,
instead of implementing the deterministic solution. Note that the calculation of the ratio requires
solving only two instances of CVRP.
The second performance measure considers the effect of the solutions on the demand when
it is subject to demand uncertainty. The ratio r
ud
is the relative unsatisfied demand for the
deterministic solution when it is implemented under the worst case realization of uncertainty. It
52
is given byr
ud
=
z
ud
P
i∈V
d
0
i
where the nominatorz
ud
is the maximum unsatisfied demand that can
occur when implementing the policy given by the deterministic solution, under the worst possible
realization of the uncertainty with respect to this policy. The denominator is the total demand of
the deterministic case and it is assumed that the deterministic problem is feasible. To obtain z
ud
,
we first solve the deterministic CVRP where the uncertainty set is fixed to the nominal values.
Thentheoptimalvalueofthedecisionvariablesxconstitutesapolicytoimplement. Now,wewant
to test the performanceofthis policy by fixingx, and by makingthe realizationofthe uncertainty
set as bad as possible for this particular policy. This provides the value of z
ud
. For that, we just
need to feed the worst possible demand values from the robust counterpart formulation to the
particular policy suggested by z
d
. Note that the calculation of the ratio requires solving only one
instance of the CVRP.
The robust counterpart gives a solution which has zero unmet demand in every possible re-
alization of the uncertainty and has the best cost in the worst case, whereas the deterministic
counterpart may have unmet demand depending on the particular scenario but it has the best
cost in the expected case. In other words, robust solution guarantees to have zero unmet demand
under uncertainty with an increase in the cost with respect to the deterministic solution. Since
the two dynamics competing here between robust and deterministic are the cost of the solution
and the amount of unmet demand, we will be using the above performance measures to compare
robust and deterministic latter in the experimental section.
It is possible to derive an upper bound on r
ud
. Observe that the maximum amount of unmet
demand which can be incurred by the deterministic solution is the worst case possible increase in
the nominal demand values due to uncertainty. If we assume that the deterministic solution is not
able to meet any of this additional demand, then the upper bound based on the type of demand
uncertainty is given by
r
ud
≤
X
i∈V
UD
i
X
i∈V
d
0
i
53
where for the case of convex hull uncertainty UD
i
= max{max
k
d
k
i
,0}, for the case of box uncer-
tainty UD
i
=
X
k
|d
k
i
|, and for the case of ellipsoidal uncertainty UD
i
=
s
X
k
(d
k
i
)
2
.
4.4 Uncertainty in Demand and Travel Times
In this section, in addition to the uncertainty in demand, we consider uncertainty in travel times.
Since the robust formulations with uncertainty in demand, RVRP
1
, RVRP
2
, and RVRP
3
are
instances of CVRP, it is possible to introduce uncertainty in travel time in addition to the un-
certainty in demand by using the approach proposed by Bertsimas and Sim (2003) for integer
programs with uncertain cost coefficients. They show that the optimum solution of the robust
counterpart can be obtained by solving a polynomial number of nominal problems with modified
cost coefficients. Since the RVRPs for our proposed uncertainty sets are instances of a general
integer program, formulating and solving the robust counterparts with independent uncertainty
in both demand and travel time is an application of this methodology. Similar to the previous
section, we also provide performance measures and bounds at the end of this section.
4.4.1 Derivation of the RVRP
Addinguncertaintyintraveltimesinadditiontotheuncertaintyindemandisequivalenttoadding
uncertainty in travel times to the CVRP formulation directly since our RVRPs under demand
uncertainty lead to CVRP instances. Here, we follow the approach introduced by Bertsimas and
Sim (2003) which solves integer problems with cost uncertainty by only a polynomial increase in
the time complexityofsolvingthe deterministic problem. This specialtypeofboxuncertaintyset,
U
C
,allowseachtraveltimec
ij
todeviatefromitscurrentvalue, ¯ c
ij
,byincuringanadditionalcost,
η
ij
, but it restricts the maximum number of arcs in the network to be uncertain, m. Without loss
54
of generality, it is assumed that the additional costs are indexed in non-increasing lexicographical
order: η
00
≥η
01
≥η
ij
≥η
nn
. Let’s further define η
n+1n
= 0.
U
C
={c| ¯ c
ij
≤c
ij
≤ ¯ c
ij
+y
ij
η
ij
, y
ij
∈{0,1}, i,j ∈V,
X
i∈V
X
j∈V
y
ij
≤m}
As Bertsimas and Sim (2003) showed, for a given value of m, the robust counterpart problem
hasanoptimalsolution,z
∗
,whichistheminimumofatmost2n+3nominaloptimizationproblems.
The nominal problems have an altered objective function with respect to the parameter m and
the additional costs. That is, for i = 0,...,n+1 and j =0,...,n the nominal problem N
ij
is given
by
(N
ij
) min
x,u
¯ cx+
i
X
k=0
j
X
l=0
(η
kl
−η
ij
)x
kl
+mη
ij
s.t. (1.2−1.8).
Then the corresponding optimal solution of the robust counterpart is z
∗
=min
i,j
N
ij
.
If f is the number of distinct values among the additional deviation costs, η
ij
, then the above
procedure solves only f + 1 distinct nominal optimization problems. If T is the time to solve
one nominal optimization problem, the procedure solves the robust counterpart in T(f +1) time,
thus preserving the polynomial solvability of the nominal optimization problem. Also note that
setting m to a bigger value than the total number of arcs in the solution, K +n, will result in
the same solution as setting m = K +n since allowing more arcs in the network to deviate than
the maximum number of arcs which could be selected in the solution can not make the robust
solution worse which reaches its worst case when those two quantities are equal.
In our framework, by adding uncertainty in demand, the CVRP will be altered to RVRP
1
,
RVRP
2
orRVRP
3
. In addition to this changein demand data, we can apply the aboveprocedure
to add uncertainty also in travel times which basically requires a change in the cost data of
f+1 nominal problems. These manipulations of data together with the above solution procedure
55
increase the total complexity of solving the RVRP with uncertainty in demand and travel times
by only a polynomial factor of the complexity of the original deterministic CVRP.
4.4.2 Performance Measures and Bounds
To compare the costof the solutions of the robust and the deterministic CVRP when the problem
includes uncertainty in cost, we can adapt two performance measures proposed by Ord´ o˜ nez and
Zhao (2004). We will both define the measures and explain how to compute them. The first
performance measure is the ratio r
wc
which is the relative improvement of the robust solution
in the worst case realization of uncertainty set. It is given by r
wc
=
z
wc
−z
r
z
r
where z
wc
is the
objective function value ofthe optimum solution of the deterministic CVRP when it is used under
the worstcase realizationofuncertainty andz
r
is the optimum solution of the robustcounterpart.
Since the latter is an optimization problem with respect to worst case realization, its solution will
be always at least as good as any other solutions implemented under the worst case. Therefore
this ratio is always non-negative.
The value z
r
is obtained by the procedure explained in the previous section. To get z
wc
, we
first solve the deterministic CVRP where the uncertainty set is fixed to the nominal values. Then
the optimal value of the decision variables x will constitute a policy to implement. Now, we want
to test the performanceofthis policy by fixingx, and by makingthe realizationofthe uncertainty
set as bad as possible for this particular policy. This will give the value of z
wc
. For that, we
need to solve the below IP formulation. Note that this is a very easy to solve IP with only one
constraint and (n+1)
2
binary variables.
max
X
i∈V
X
j∈V
( ¯ c
ij
+η
ij
y
ij
)x
ij
s.t.
X
i∈V
X
j∈V
y
ij
≤m
y
ij
∈{0,1} i,j ∈ V
56
It is possible to derive an upper bound on r
wc
by following the methodology shown by O˜ guz
(2000). If there is a constant σ such that
η
ij
¯ c
ij
≤σ holds for all (i,j)∈A, then the bound is
r
wc
≤
2σ
1+σ
.
The second performance measure is the ratio r
nc
which is the relative loss of optimality of
the robust solution when implemented under the nominal case realization of uncertainty set. It is
given by r
nc
=
z
nc
−z
d
z
d
where z
nc
is the objective function value of the optimum solution of the
robust counterpart problem when it is used under the nominal case realization of uncertainty and
z
d
is the optimum solution of the deterministic CVRP. Following a similar reasoning we did for
the ratio r
wc
, it is easy to see that the above ratio r
nc
is always non-negative.
The value z
d
is obtained just by solving the deterministic CVRP. To get z
nc
, we first solve
the robust problem optimally and fix the value of the decision variables suggested by the policy.
Now, we want to test the performance of this policy when the realization of the uncertainty set
is exactly equal to the deterministic case. Then, to compute z
nc
, we just need to recalculate the
objective function value in (1.1) using the policy suggested by the robust solution.
4.5 Experimental Results
In this section, we focus only on the demand uncertainty and consider the related performance
measures that will be used to compare robust and deterministic solutions. We then analyze the
trade-offs of robust solutions on instances from the literature and families of clustered instances.
We finish the section verifying that the robust solution provides better protection to the demand
uncertainty than a uniform and a non-uniform distributions of the unused capacity among the
vehicles. We present computational results only for the convex hull uncertainty set. Similar
results are obtained for box and ellipsoidal uncertainty sets.
57
Our experiments require the solution of the deterministic and robust versions of the routing
problem. Since all are instances of CVRP we use an efficient exact solution procedure for generic
CVRP. We use the branch-and-cutbased VRP solverin the open sourceSYMPHONY librarydue
to Ralphs et al. (2003), available on-line at http://branchandcut.org/VRP.All experiments are
carried out with a runtime limit of one hour on a Dell Precision 670 computer with a 3.2 GHz
Intel Xeon Processor and 2 GB RAM running Red Hat Linux 9.0.
4.5.1 Robust versus Deterministic on Standard Problems
Our first set of experiments addresses problem set A (Random Instances), set B (Clustered In-
stances), and set P (Modified Instances from the literature) of the CVRP suite of problems by
Augeratetal. (1995). Theinstancesrangefrom15to 100customers. Wemodified these instances
to include demand uncertainty. We allow eachdemand parameterto further increaseup to a fixed
percentage of the deterministic value. We randomly generate a total of 5 scenarios within the
allowed percent deviation for the demand uncertainty set. More specifically, we use the following
values of percent deviation in demand parameters: 5, 10, 15, and 20.
Table 4.1 shows the results based on the performance measures r
ud
for the percent unmet
demand ratio and r
rd
for the percent cost ratio of the solutions, where “No” indicates the number
of the instance, “T” indicates the percent tightness ratio of the instance which is defined as the
ratio of the total expected demand to total vehicle capacity, “IN” indicates infeasible instance,
and “NA” indicates that an optimal solution could not be found within the 1 hour runtime limit.
The first observation is that since the original instances are already tight (they have a percent
tightness ratio between 81% and 99%), the robust counterparts run quickly into infeasibility as
the percent deviation of uncertainty increases even though the deterministic CVRPs are feasible
and could be solved to optimality. In almost all of the instances with a tightness ratio greater
than 90%, the cost ratio could not be calculated for the percent deviation values 15% and 20%
since either the robust counterpart became infeasible or the runtime limit was reached. Recall
that the robust instances are more capacity constrained due to increased demand in the data
58
Table 4.1: Augerat et al. (1995) sets
Set A Set B Set P
Unmet Cost Unmet Cost Unmet Cost
Demand Ratio Demand Ratio Demand Ratio
No T 5 10 15 20 5 10 15 20 T 5 10 15 20 5 10 15 20 T 5 10 15 20 5 10 15 20
1 82 0.5 2.5 4.8 6.9 0 1.5 4.2 9.3 82 0 0.7 2.6 4.5 0 1.3 1.6 3.3 87 0.5 1.2 2.2 3.3 0 2.2 5.1 IN
2 89 0.2 1.7 3.5 5.9 0 2.1 6.5 IN 91 1 2.5 4.4 6.5 0 0.1 5.5 IN 96 0.5 2.3 5 8.2 0 IN IN IN
3 90 0.3 1.8 4 6.5 0 2.7 7.1 IN 87 0.7 2.3 4.4 6.7 0.1 2.8 3.2 4.4 96 0.8 2.2 4.8 7.5 0 IN IN IN
4 92 0 0.3 1.6 3.2 0 1.5 4.6 IN 85 0.4 2 4.2 7.3 0 0.1 1.7 1.9 93 0 0 0.9 3.4 0 0 0 IN
5 88 0.8 2.1 4 6.3 0 1.8 4.5 6.3 88 1.2 2.7 4.6 6.7 0 2.2 2.6 6.4 96 0.9 2.4 3.8 5.9 0 3.2 IN IN
6 81 0.2 1.6 3.3 4.8 0 2.7 4.8 5.4 94 0.9 2.4 4.3 6.7 0 4.2 IN IN 93 0.7 1.9 3.3 4.9 14.9 IN IN IN
7 95 0.7 2.4 4.5 6.8 0.1 NA IN IN 86 1.1 3.5 6.4 9.4 0 1.5 4.9 5.5 97 1.2 3.4 5.8 8.2 2.8 IN IN IN
8 96 1.4 3.8 6.4 9 0 IN IN IN 91 1.2 3.6 6 8.5 0 4.6 8.5 IN 88 0.3 1.2 2 3.1 0 0.7 2.4 IN
9 95 1 3.1 5.2 7.9 0 1.5 IN IN 97 0.9 3.1 5.3 8.1 0 IN IN IN 92 1.2 3.3 5.5 7.6 0.6 2 IN IN
10 87 1.3 4.2 7.3 10 0 0.4 2.3 5.1 98 1.8 4.4 6.9 9.3 0 IN IN IN 90 0.5 2 4.1 6.1 0.5 2 NA IN
11 95 0.9 2.2 3.8 5.8 0 1.5 IN IN 87 0.6 1.5 2.8 4.2 0 0.5 1.3 1.9 99 1.5 3.9 6.7 9.3 IN IN IN IN
12 98 1 3.7 6.4 9.1 0 IN IN IN 91 1.2 3 5.3 7.4 NA 2.7 NA IN 95 0.7 2.1 4.3 6.7 NA IN IN IN
13 90 0.8 2.4 4 6 NA NA NA IN 97 1.2 3.2 5.8 8.7 0 IN IN IN 97 1.1 3.2 5.3 7.7 NA IN IN IN
14 86 1.4 3.4 5.7 8.1 0 4.2 6.8 8 86 0.7 2 3.7 5.7 0 0.7 1.5 3.2 87 0.4 2 4.3 6.9 NA NA NA NA
15 89 1 2.6 4.5 6.6 0 NA NA IN 88 1.1 3 4.9 6.9 0 3 3.4 5.8 81 0.5 1.9 3.6 5.3 NA NA NA NA
16 94 1.1 3.3 5.7 8.2 0.5 3.7 IN IN 99 2 5 8.1 11 IN IN IN IN 90 0.8 2.9 5.1 7.7 NA NA NA IN
17 95 0.4 1.9 4.5 7 NA NA IN IN 89 1.2 3.1 5.3 7.5 NA NA NA NA 99 1.6 4.4 7.2 10 NA IN IN IN
18 93 1.3 3.4 5.8 8.5 0.1 2.8 IN IN 92 1.1 3 5.1 7.4 NA NA NA IN 94 0.5 1.7 3.7 5.9 NA NA IN IN
19 92 0.9 2.1 3.8 5.8 NA NA NA IN 97 1.2 3.3 5.6 8 0 IN IN IN 94 0.8 2.3 4 5.9 NA NA IN IN
20 98 1.7 4.4 7.2 10 NA IN IN IN 95 1.4 3.6 5.9 8.4 NA IN IN IN 93 0.6 1.8 3.3 5.3 NA NA IN IN
21 91 1.3 3.6 6.1 8.6 NA 4.7 NA IN 90 1 2.5 4.4 6.5 NA NA NA IN 97 0.5 2.3 4.4 6.6 NA IN IN IN
22 97 1.2 3.1 5.4 7.7 NA IN IN IN 93 1.6 3.9 6.4 9.1 NA NA IN IN 97 0.8 2.4 5 7.7 0.5 IN IN IN
23 93 1.1 3.3 5.6 8.1 NA NA IN IN 93 1.4 3.7 6.2 8.9 NA NA IN IN 97 0.9 2.5 5.1 7.9 NA IN IN IN
24 94 1.2 3.5 5.9 8.4 NA NA IN IN 91 0 0.9 2.9 5 0 0.4 0.6 IN
25 97 1.3 3.8 6.3 8.9 NA IN IN IN
26 93 0.8 2.5 4.5 6.8 NA NA IN IN
27 94 0.7 2.5 4.4 6.6 NA NA IN IN
59
and empirical observations have shown that CVRPs become more difficult to solve in practice
as problems are more capacity constrained. In the interpretation of Table 4.1, we remove “IN”
instances since by definition our robust optimization model does not account for infeasibility of a
given scenario and these instances require a different model as discussed before; we also remove
“NA” instances even though in several of these cases a best feasible solution could be calculated
for robust, since comparing a robust feasible solution with a deterministic optimal solution would
be bias. Therefore, we only interpret cases where both the robust and deterministic could be
solved optimally on feasible instances.
To compare the ratios in each problem set to analyze the average behavior of the robust and
thedeterministic,wecountthenumberoftimesthecostratioisnotworsethantheunmetdemand
ratio. When the increase in the unmet demand on the worst case scenario for the deterministic
solution is greater than the increase in cost for the robust solution on the nominal case, we
say the robust solution is preferable. We understand that the magnitude for comparison of the
ratios is application and scenario dependent. For simplicity, we weight both ratios equally in our
comparisonsin this study. In setA, robustispreferablein 27casesoutof42possiblecomparisons;
in setB,itis 39casesoutof45; andinsetP,itis16outof23. We seethatthe robustoutperforms
the deterministic the most in set B where the instances are clustered; that is the additional cost
of implementing the robust is no worse than the amount of unmet demand incurred in the worst
case if the deterministic was implemented.
A close analysis of the solutions revealed that the success of the robust for the clustered
instances is directly related to the distribution of the unused vehicle capacity slack over the
network. When the realized demand is greater than the expected demand, the additional slack
helps address the increase in demand. If there is not enough slack in the vehicle for the additional
demand, then the robustwill routethe vehiclesdifferently. In this case, if there is a vehiclenearby
with enough slack, then this demand can be serviced by this latter vehicle with only a slight cost
increase. However, if a vehicle with enough slack is far from this extra demand, then the robust
solution that can satisfy this demand will be significantly more expensive.
60
Consider for example the deterministic and robust optimal solutions for instance number 3
from set B depicted in Figure 4.2. The numbers next to the routes indicate the total expected
demand, d
0
serviced on each route. The slack of each route is the vehicle capacity, 100, minus
this total expected demand. We label the customers according to the deterministic route that
services them, keeping the label for the robust route to emphasize the change in solution. We
see that five units of slack from the vehicle which services 88 units of expected demand in the
deterministic is distributed in the robust to the other vehicles without significantly increasing the
total cost, even though itaffected three other routes. Also note thatthe routeofthe vehiclewhich
services 56 units of total demand in the deterministic stays the same in the robust since there
is already enough slack in this vehicle to accommodate the increase in the demand due to worst
case realization of the uncertainty. Overall, for this particular case, the increase in the cost due
to robust is less significant than the total unmet demand that occurs in the worst case realization
of the uncertainty for the deterministic.
The Deterministic
100
88
56
97
96
The Robust
96
93
96
56
96
Figure 4.2: Deterministic and robust solutions for Augerat et al. (1995) set B instance 3. No. of
vehicles 5, capacity 100, value next to route is the deterministic demand served
On the other hand, when the network structure does not allow an easy distribution of slack,
then the robust may result in a poor performance as depicted in Figure 4.3, instance number 1
from the set B. The vehicles have again the identical capacity of 100. In this case, the slack in the
vehicle which services 38 units of total demand in the deterministic is costly to distribute in the
61
robust to the other vehicles since the route of that vehicle is relatively far from the routes of the
others.
The Deterministic
97
86
96
38
95
The Robust
82
91
96
87
56
Figure 4.3: Deterministic and robust solutions for Augerat et al. (1995) set B instance 1. No. of
vehicles 5, capacity 100, value next to route is the deterministic demand served
Theseexamples showthatthe structureofthe networkplaysakey rolein determining optimal
routes and thus the distribution of the unused capacity and its impact on the success of a robust
solution. The instances in the Augerat et al. (1995) suite of problems suggest that when the
network is clustered optimal solutions will have some vehicles close which could share unused
capacity at a low cost favoring a robust solution.
4.5.2 Robust versus Deterministic on Family of Clustered Instances
To validate our findings and to generalize them with respect to the structure of the network,
we randomly generate instances with 4 vehicles of capacity 1500 and 49 customers with uniform
demand of 100, in three different problem sets. In each set, there are 4 clusters of customers.
For the clustering, we generalize the clustering idea in Chapter 3. First of all, we consider points
which are on the circle of a given radius R, centered at a depot, and we randomly select a point
on that circle to be the center of a cluster (see Figure 4.4). Then we generate customers for that
cluster within the circle with a given radius r. We also use as before the measure for clustering r
c
which is given by: r
c
=
R
r
. We fix the value of r = 20 and consider the following values for R: 0,
2r, 4r, 6r, and 8r.
62
R
r
Depot
Center of Cluster
Figure 4.4: Cluster generation in random sets
When r
c
= 0 all the clusters are centered at the depot and the instance becomes random with
no clustering effect; however as r
c
increases, clusters separate from each other. In problem sets 1
and 2, the three clusters have 13 customers and the fourth one has 10. Note that a vehicle can
service up to 15 customers. The reasonfor this selection is that, in clustered instances where each
cluster will be serviced by only one vehicle, there will be one vehicle with relatively more slack,
namely the one servicing the fourth cluster. The only difference between sets 1 and 2 is that in
the latter as we increase r
c
, we always keep the fourth cluster centered at the depot. This serves
the purpose of having a random zone around the depot and some clusters far from the depot. In
set 3, we make this random zone around the depot denser by increasing the number of customers
in the fourth cluster to 25 and decreasing the one for the others to 8. Figures 4.5, 4.6, and 4.7
display the results of the three sets for percent unmet demand ratio r
ud
and percent cost ratio
r
rd
as a function of percent clustering ratio r
c
for different values of percent deviations of the
uncertainty set. Each data point on the figures is an average of 30 instances. For all these three
sets, we observedthatthe nature ofrobustsolution is significantlydifferentthan the deterministic
solution in most of the cases when we calculated the average ratio of the number of different arcs
in both solutions to total number of arcs. Depending on the percent deviation in demand, this
difference varies from 3% to 33% for set 1, from 20% to 34% for set 2, and from 24% to 33% for
set 3.
63
0 2 4 6 8
0
2
4
6
8
10
12
14
Average Unmet Demand Ratio
Clustering Ratio
Uncertainty
Range
5%
10%
15%
20%
0 2 4 6 8
0
2
4
6
8
10
12
14
Average Cost Ratio
Clustering Ratio
Uncertainty
Range
5%
10%
15%
20%
Figure 4.5: Comparison of deterministic and robust solutions for random set 1
0 2 4 6 8
0
5
10
15
20
25
30
Average Unmet Demand Ratio
Clustering Ratio
Uncertainty
Range
5%
10%
15%
20%
0 2 4 6 8
0
5
10
15
20
25
30
Average Cost Ratio
Clustering Ratio
Uncertainty
Range
5%
10%
15%
20%
Figure 4.6: Comparison of deterministic and robust solutions for random set 2
64
0 2 4 6 8
0
1
2
3
4
5
6
7
8
9
10
11
Average Unmet Demand Ratio
Clustering Ratio
Uncertainty
Range
5%
10%
15%
20%
0 2 4 6 8
0
1
2
3
4
5
6
7
8
9
10
11
Average Cost Ratio
Clustering Ratio
Uncertainty
Range
5%
10%
15%
20%
Figure 4.7: Comparison of deterministic and robust solutions for random set 3
The results of set 1 suggestthat both the deterministic and the robustbenefit from clustering.
For percent deviation of 5% and 10%, both results are comparable, for 15% the robust is better,
andfor20%the robustis worsewhenr
c
≥ 2. Infactastheuncertaintyincreases,wewouldexpect
the robust to outperform the deterministic. The reason for this odd behavior is the distribution
of the slack in the network. When the instances are clustered for the bigger values of r
c
, each
cluster is serviced only by one vehicle in the deterministic. In case of high uncertainty such as
20% deviation, if the total demand of a cluster exceeds the vehicle capacity then another vehicle
has to be routed to this cluster by the robust. When these vehicles are not close, the robust
results in a large travel cost. The network structures with pure clusters as in set 1 therefore do
not allow a good distribution of slack on the average and the robust is not convenient for high
uncertainties. An illustrative example of such a phenomenon is depicted in Figure 4.8. We see
that the slack in the vehicle which services 1300 units of total demand in the deterministic is
distributed in the robust to the clusters which are relatively far from the route of that vehicle.
65
Thus, in this particular case, the percent cost ratio due to the robust, 12.81%, is more significant
than the unmet demand ratio due to the deterministic, 4.2%.
The Deterministic
1000
1100
1500
1300
The Robust
1300
1100
1200 1300
Figure 4.8: Deterministic and robust solutions for a bad case for robust from random set 1. No.
of vehicles 4, capacity 1500, value next to routes is the deterministic demand served
When we look at the results of set 2, as before we see the same phenomenon in the increase
of the cost ratio for the robust with 20% deviation. However, clustering helps only after r
c
> 2.
The reason is due to the random zone around the depot. When r
c
≤ 2, the circles of clusters
intersect and the vehicles do not necessarily service only customers for the same cluster. This
interaction of customers keeps the instance as random until r
c
>2 since from that point onwards
the three clusters become more distinct than the fourth one around the depot and the effect of
clustering gets more pronounced in the instance. Increasing r
c
until 2 only makes the size of the
network enclosing all the customers bigger, and therefore the cost of robust increases on these
bigger random instances. When it comes to the amount of unmet demand of the deterministic,
the effect of the random zone is more drastic. No matter how much the network is clustered,
the unmet demand is always constant and much worse compared to set 1. The vehicles in the
deterministic service customers in the random zone on their way to the clusters and usually 3 out
of4vehiclesarefilledtocapacity,whichisnotthecaseinthe deterministicofset1. Thesevehicles
with full capacity are the minimum cost solutions but they have a very big potential of incurring
unmet demand under uncertainty. The network structures with a scattered random demand zone
around the depot as in set 2 therefore have a very negative effect on the deterministic.
66
When the random zone is denser around the depot, the results of set 3 are similar to the
results of set 2. The deterministic results in high unmet demand values and is outperformed by
the robust in almost all the cases. Having more customers in the random zone helps the robust
even further. The reason why the phenomenon with 20% uncertainty disappearsis due to the fact
that the vehicles are close now since the number of customers they service in the random zone
on their way to the clusters is significantly larger compared to set 2. Therefore when the slack
in one vehicle needs to be distributed to the network, this can be achieved through customers in
the random zone. The network structures with clusters and dense random zone around the depot
as in set 3 therefore allow a good distribution of slack on the average and the robust solution
benefits from this with little extra cost. An illustrative example of such a phenomenon is depicted
in Figure 4.9. We see that the slack in the vehicle which services 700 units of total demand in the
deterministic is easily distributed in the robust to the other vehicles by exchanging the customers
in the random zone near the depot. Thus, in this particular case, the percent cost ratio due to the
robust, 0.78%, is less significant than the unmet demand ratio due to the deterministic, 6.76%.
The Deterministic
1500
1200
700
1500
The Robust
1300
1000
1300
1300
Figure 4.9: Deterministic and robust solutions for a good case for robust from random set 3. No.
of vehicles 4, capacity 1500, value next to routes is the deterministic demand served
Lastly, Figure 4.10 depicts results for the three uncertainty models (convex hull, box, and
ellipsoidal) on the random set 3 with 5% deviation in demand. As expected, box results in the
mostunmetdemandandconvexhulltheleast. Basically,ellipsoidalandboxare“shifted”versions
67
of convex hull uncertainty. Therefore the general trends observed for convex hull is carried out to
the other two uncertainty sets whose results are omitted in general due to space limitations. In
this particular case, the deterministic solution is not attractive at all under box and ellipsoidal
whereas for convex hull it is acceptable for values of clustering ratio r
c
< 4. This is due to the
fact that the increase in the nominal demand is much more drastic for box and ellipsoidal in the
worst case making the performance of the deterministic solution very poor.
0 2 4 6 8
0
1
2
3
4
5
6
7
8
9
10
Average Unmet Demand Ratio
Clustering Ratio
Uncertainty
Model
CH
BOX
ELL
0 2 4 6 8
0
1
2
3
4
5
6
7
8
9
10
Average Cost Ratio
Clustering Ratio
Uncertainty
Model
CH
BOX
ELL
Figure 4.10: Comparison of three uncertainty models for random set 3
We conclude this experimental subsection by emphasizing that our findings in the instances
by Augerat et al. (1995) are confirmed by a larger class of random instances from the population
of instances with the same characteristics. In particular, we showed that both the existence of
enough slack in the solution and its distribution over the network are very important factors
affecting the quality of the robust. Our experiments reveal that clustered network structures with
a dense random zone around the depot favor the robust. For this scenario, we showed that the
deterministic could result in a large amount of unmet demand and the extra cost of the robust is
relatively small.
68
4.5.3 Robust versus Distributions of Excess Vehicle Capacity
The robust solution distributes the excess vehicle capacity in the expected demand case aiming to
obtain routes at minimum cost that satisfy all demand outcomes from the uncertainty set. In this
section we explore how this compares to two simple strategies of distributing this excess capacity
among all the vehicles: uniformly and non-uniformly.
Werandomlygenerateinstanceswith4vehiclesofcapacity2100and68customerswithuniform
expected demand of 100. Therefore there is a total of 1600 units of excess vehicle capacity to be
used to address the demand uncertainty. We generate these instances according to the three sets
as before. The uniform distribution of excess capacity will reserve a buffer of 0, 100, 200, 300,
and 400 units of capacity in each vehicle. The distribution of the remaining amount of buffer
(except for the case of 400 units which distributes all the available excess capacity uniformly) is
automaticallydeterminedbythe solutionproceduretominimizethe totalcostofthe deterministic
solution based on the particular instance. That is, the uniform buffer amount (excess capacity)
is removed from consideration to compute the optimal deterministic solution, but it is considered
when determining the amount of unmet demand that this optimal solution can face in its worst
case. On the other hand, the non-uniform distribution of excess capacity uses the same settings
butallowsassigningthepre-determinedreserveofonevehicletoanother,thatistheuniformbuffer
amount of one vehicle can be doubled at the expense of not having a pre-determined reserve in
another vehicle. Since there are 4 vehicles in total, at most 2 of them can have doubled pre-
determined reserve. This assignmentas well as the distribution of the remaining amount of buffer
are again automatically determined by the solution procedure.
To compare the quality of the solutions, we use the percent unmet demand ratio r
ud
but
generalize the percent cost ratio to r
0
rd
=
z
r
−z
bc
z
d
. Here we simply replaced z
d
in the numerator
with z
bc
, which is the optimal objective value for the deterministic solution with reduced vehicle
capacity (buffer capacity). Figures 4.11 and 4.12 display respectively for the uniform and non-
uniformdistributions theaverageresultsover30instancesfordifferentvaluesofpercentclustering
ratio r
c
for set 1 with 15% deviation in the uncertainty scenarios. Similar trends are observed in
69
the other randomly generated sets and percent deviation values for our three types of uncertainty
sets.
0 100 200 300 400
−1
0
1
2
3
4
5
6
7
Buffer Amount
Average Unmet Demand Ratio
Clustering
Ratio
0
2
4
6
8
0 100 200 300 400
−1
0
1
2
3
4
5
6
7
Buffer Amount
Average Cost Ratio
Clustering
Ratio
0
2
4
6
8
Figure 4.11: Comparison of uniform buffer capacity and robust solutions for random set 1
In Figure 4.11, fora given value ofpercent clusteringratio, it is clear thatincreasingthe buffer
amount makes a uniform distribution of slack have less unmet demand but with an increased cost
which may exceed the cost of the robust solution in some cases, giving negative values for the
percent cost ratio r
0
rd
. When we compare the quality of the two solutions, we see that when the
buffer amount is smaller than 200, this reserve capacity is insufficient to handle the uncertain
demand. When the buffer amount is equal to 200, the uniform distribution of slack leads to a less
costly solution than the robust with the same zero unmet demand. When the buffer amount is
equal to 300, the two methods have the same cost with the same zero unmet demand. After this
transition point (when the buffer amount is greater than 300), if we increase the buffer capacity
unnecessarily, the resulting solution is more costly than the robust and is not preferable. These
trends become less pronounced as the percent clustering ratio increases. That is, clustering is
good for a uniform distribution of slack, which makes sense since such an even distribution of
slack benefits by having each vehicle assigned to distinct, far awayclusters with the same demand
70
and uncertainty, as in the case of r
c
= 8. However, if the instance is random (r
c
= 0), then an
even distribution of slack is not necessarily advantageous and the robust solution allocates the
slack in a better way by taking into account the particular structure of the network.
0 100 200 300 400
−1
0
1
2
3
4
5
6
7
Buffer Amount
Average Unmet Demand Ratio
Clustering
Ratio
0
2
4
6
8
0 100 200 300 400
−1
0
1
2
3
4
5
6
7
Buffer Amount
Average Cost Ratio
Clustering
Ratio
0
2
4
6
8
Figure 4.12: Comparison of non-uniform buffer capacity and robust solutions for random set 1
In Figure 4.12, the first observation is that the cost ratios are higher than the uniform case,
meaning that the optimal deterministic solution of non-uniform distribution is less costly than
the one of uniform distribution. This is expected since non-uniform distribution explores a wider
set of feasible solutions by not necessarily forcing a uniform pattern in the solution. The second
observationisthatthedemandratiosarealsohigherthantheuniformcase. Thisisinfacttheprice
ofhavingalesscostlydeterministicsolutionfornon-uniformcase. Duetoparticularassignmentsof
pre-determined buffer amount among vehicles, non-uniform distribution creates less costly routes
where some of the vehicles havedoubled pre-determined buffer capacitywhereassome others have
no pre-determined buffer capacity which in return may result in unmet demand. This effect of
having an uneven distribution of pre-determined buffer capacity is much more pronounced in the
caseof 400units ofbuffer capacitysince all the available1600units ofextracapacity areallocated
to vehicles by this non-uniform pattern. That is, when the solution procedure chooses to double
71
thebuffercapacityofavehicle,anothervehiclewillhaveexactlyzeroslackinthesolution. Overall,
this particular method of non-uniform distribution does not outperform the robustor the uniform
distribution although it is able to improve the cost ratio over uniform distribution.
4.6 Summary
In this chapter, we propose to use robust optimization to obtain efficient routing solutions for
problems under uncertainty. Our work has shown that robust optimization is an attractive al-
ternative for formulating routing problems under uncertainty as it does not require distribution
assumptions on the uncertainty or a cumbersome representation through scenarios.
We derived three robust counterparts for the VRP with uncertainty in demand and showed
how to derive a robustcounterpartwith uncertainty in both demand and traveltimes fora special
type of travel time uncertainty. Our VRP formulation and definition of uncertainty sets resulted
in computationally amenable RVRPs. In the case of demand uncertainty, we need to solve only a
single CVRP with modified demand data to obtain the solution of the robust counterpart and in
the caseofdemand and traveltime uncertaintyonlyapolynomialnumberofCVRP with modified
demand and travel time data.
For the case of demand uncertainty, we used an open source VRP solver to experimentally
investigate the trade-offs between a robust and deterministic solutions in terms of the increased
cost of the robust solution and the possible unmet demand in the worst case of the deterministic
solution. We firstsolved instances from the literature and postulated some insights about the net-
work structures affecting the quality of the robust solution. We then generated random instances,
with particular network structures and different degree of clustering, to validate our findings.
Our results showed that if the network structure allows a strategic distribution of the slack in
the vehicles throughout the network in such a way that the vehicles can collaborate with ease by
sharing their slacks in case of uncertainty, then the robust solution is favorable on average. Such
72
a network structure appears in a problem with clustered zones far from the depot with a dense
random zone near the depot.
We also verified that the robust solution is superior than uniformly distributing among all
vehicles the excess vehicle capacity under the expected demand. We showed that such a solution
only competes with the robust solution if the network structure is highly clustered where the
expected demand of each cluster is about the same. Lastly, we showed that a particular method
of non-uniformly distributing this excess vehicle capacity is not a viable solution.
73
Chapter 5
Robust Optimization with Recourse for VRP
5.1 Methodology
We have addressed the challenges of developing RVRP in the previous chapters and in particular
we showed how to formulate and solve a RVRP with demand and travel time uncertainty. In this
chapter, we address the challenges of developing a recourse model for RVRP that is tailored to a
specificapplication. WefirstdevelopanaprioriRVRPfollowingrobustoptimizationmethodology
inthesamespiritofthepreviouschapterandthenweaugmentthemodelbyincorporatingrecourse
actions following scenario based stochastic optimization. As a result, we benefit from both the
simplicity of robust optimization and the flexibility of recourse models. For the last challenge of
solving the problem efficiently, we develop an insertion based heuristic to address a large scale
courier delivery problem.
Our problem domain is the courier delivery problem which is a variant of the VRP with
time windows. This problem and its variants are faced by several daily courier companies in
the industry. The particular characteristics of the problem we address in this study follow. The
delivery requests arrive daily from potential customers with uncertain time windows and service
times at the beginning of each day. The locations of the customers are known but it is not known
a priori whether a particular customer requests a delivery at a given day.
74
Thefirstgoalistoconstructadailyscheduleofcourierstoserveasmanycustomersaspossible
by meeting the deadlines. The vehicles are allowed to wait at a customer site if they arrive earlier
thantheearlieststarttimewhichmightcauseanearlinesspenalty; andiftheyarrivelaterthanthe
latest start time then a lateness penalty might be incurred. Therefore, the additional challenges
in constructing daily schedules are to minimize these penalties and too minimize the arrival time
of the vehicles to the depot at the end of the day, which also accounts for the tour length of each
vehicle in time units including waiting and travel times.
The second goal is to construct an a priori master plan to be used in constructing daily
schedules by making modifications every day according to daily customer requests. The reason is
twofold. First, it is desired to have similarity in daily schedules to minimize the cost of re-planing
every day. In particular, similarity in daily schedules allows for better familiarity in the routes
by the drivers, thereby potentially improves the driver familiarity. Second, for a given planning
horizon, the company will need to plan its resources and budget beforehand. The master plan
constitutes a basis for the estimate of such resource and cost allocations.
To formulate the problem, we adapt the VRPTW formulation and use a robust optimization
approach to develop a RVRP with uncertainty in service times and time windows. This also
provides the master plan which has robust routes to service customers with high probability of
occurrences. Inaddition,weusescenariobasedstochasticoptimizationtoaddresstheprobabilistic
nature of the customers by considering each day in the planning horizon as a scenario. Thus, we
incorporate recourse actions in the problem formulation which provides the daily schedules. The
objective of the problem is, for each scenario, to maximize the coverageof customers, to minimize
the total length of the routes and the total penalty, and to maximize the similarity of routes with
the master plan.
To solve the problem, we develop a two-phase approximate solution procedure using insertion
as an efficient heuristic subroutine to address the large scale real life problem. The first phase is
the construction phase where a master plan and initial daily schedules are generated. The second
phase is iterative. At each iteration, first the master plan is modified based on the feedback given
75
by the daily schedules and second the daily schedules are updated based on the new master plan.
The second phase and thus the two-phase algorithm end when there is no further improvement in
the objective function at the end of an iteration reaching a local optimum solution.
In the next three sections, we gradually develop the problem formulation. Section 5.5 presents
the problem data and the nature of uncertainty. In Section 5.6 we propose the two-phase solution
procedure and we present the experimental results in Section 5.7. The last section concludes the
chapter.
5.2 Deterministic VRP Formulation
In this section, we formulate a deterministic VRP for the courier delivery problem which is a
variantofthe VRPTW. Given a networkof customers and the location ofthe depot, the objective
is to routethe fleetofvehiclesto servicethe maximum number ofcustomersbasedon their service
times and time windows. We consider the soft time windows in this study. That is, if a vehicle
arrives to a customer earlier than its earliest start time, there will be a waiting penalty; and
similarly if it arrives later than its latest start time, there will be a lateness penalty. The vehicles
start and end their routes at the depot. The length of a route is composed of travel and waiting
times. There is a common due date for vehicles to return to the depot at the end of the day which
is a hard constraint. A customer can be serviced by only a single vehicle and there is no capacity
considerations in the problem. The primary objective is to maximize the number of customers
serviced. The secondary objectives are to minimize the route length, the earliness penalty, and
the lateness penalty.
The mathematical formulation of this problem is given below in Problem 4. K is the number
of vehicles. Let C be the set of customer nodes in the network with |C| = c, corresponding to
customers indexed as: 1...c. Then V is the set of all the nodes in the network with |V| = n =
c+1+K, including the depot at index 0 and K artificial nodes (identical copies of the depot).
Also let set A be the indices of these artificial nodes: c+1...c+1+K. The time to traverse
76
the arc from node i to node j is given by t
ij
and s
i
is the service time of node i. In particular,
s
0
= 0; ∀i ∈ A, t
i0
= t
0i
= s
i
= 0; ∀i ∈ A,∀j ∈ C, t
ij
= t
0j
and t
ji
= t
j0
; and ∀i,j ∈ A, i 6= j,
t
ij
=t
ji
=0. a
i
is the earliest start time and b
i
is the latest start time to service customer i. The
common due date for all vehicles is L. If customer i is selected to be serviced in the solution, then
z
i
is one, otherwise z
i
is zero. If the arc from node i to node j is selected in the solution then x
ij
is one, otherwise x
ij
is zero. The continuous variable y
i
is the arrival time to node i except the
depot in which case it is the departure time, i.e. y
0
= 0. In particular, y
i
for i∈A corresponds to
the arrival time to the depot of a vehicle. Note that y
i
= 0 for a customer i which is not serviced
in the solution, i.e. when z
i
= 0. The continuous variable e
i
is the earliness penalty and l
i
is
the lateness penalty of customer i; similarly e
i
= l
i
= 0 when z
i
= 0. M
ts
is a sufficiently large
number which can be set to M
ts
= L+max
(i,j)
t
ij
+max
i
s
i
; also M
ts
≤ 3L. Similarly, M
a
can be
set to M
a
= max
i
a
i
; also M
a
≤ 3L. The objective is to cover the maximum number of customers
and to minimize the total length of the routes and the total earliness and lateness penalty. The
positive constants α
1
, α
2
, α
3
and α
4
should be set according to the priority of these four goals.
In particular, α
1
should be adjusted with respect to α
3
and α
4
so that not visiting a customer
and hence avoiding earliness or lateness penalty would not be desirable. Note that in many cases
α
4
= 0 since there is no explicit penalty for waiting time but it indirectly increases the route
lengths.
The objective function (4.1) maximizes the coverage of customers while minimizing the total
route lengths and the total penalties (earliness and lateness). Note that minimizing the total
lengthsincludesminimizingthetraveltimesandwaitingtimes,whichisdifferentthanthestandard
VRP formulations which consider only minimization of travel times. Constraints (4.2-4.5)are the
routingconstraints. Constraint(4.6)forcesan artificialnode tobe servicedatthe end ofthe route
to keep track of the route length of vehicles. Constraint (4.7) is in the family of the MTZ subtour
elimination constraints. Additionally, it enforces the arrival time constraints when a customerj is
serviced right after customer i. Otherwise, it is redundant. Constraint (4.8) imposes the earliness
penalty, which is redundant if customer i is not serviced in the solution. Constraint (4.9) imposes
77
the lateness penalty and constraint (4.10) the common due date of vehicles. Lastly, constraints
(4.11-4.15) are bounds on the variables.
max α
1
X
i∈C
z
i
−α
2
X
i∈A
y
i
−α
3
X
i∈C
l
i
−α
4
X
i∈C
e
i
(4.1)
s.t.
X
j∈V, i6=j
x
ji
=z
i
i∈C (4.2)
X
j∈V, j6=i
x
ij
=z
i
i∈C (4.3)
X
i∈V\{0}
x
0i
=K (4.4)
X
i∈V\{0}
x
i0
=K (4.5)
x
i0
= 1 i∈A (4.6)
y
i
+t
ij
+s
i
≤y
j
+M
ts
(1−x
ij
)
i∈V
j ∈V \{0}, i6=j
(4.7)
a
i
≤y
i
+e
i
+M
a
(1−z
i
) i∈C (4.8)
y
i
≤b
i
+l
i
i∈C (4.9)
y
i
≤Lz
i
i∈C (4.10)
e
i
≥ 0 i∈C (4.11)
l
i
≥ 0 i∈C (4.12)
x
ij
∈{0,1} i,j ∈V, i6=j (4.13)
y
i
≥ 0 i∈V (4.14)
z
i
∈{0,1} i∈C (4.15)
(4)
5.3 A Priori Robust VRP Formulation
In this section, we consider service times, s
i
, and time windows, a
i
and b
i
, in Problem (4) to be
uncertain togetherwith probabilisticcustomers in orderto formulate an a prioriVRP formulation
whichistobeused lateringeneratingamasterplanwhichinreturnwillbe thebasisofgenerating
daily schedules. A variety of approaches can be used in addressing these uncertainties. However,
78
since one of the primary goals in formulating the master plan is to improve similarities in daily
schedules which are different for each outcome of the uncertainty, an intuitive approachis to use a
method which would resultin a masterplan which would stay “stable”for all possible realizations
of the uncertainty and thus which will require the least modifications in the resulting daily sched-
ules. This suggests using a robust optimization approach since using a master plan generated by
an expected value approach for instance may require significant modifications in generating daily
schedules for certain realizations of the uncertainty resulting in a poor performance of similarity.
Here, we use a robust optimization methodology to address the uncertainties in service times
and time windows as in the previous chapter. We treat each uncertain parameter independently
from each other and for each one we consider as before deviation vectors around a nominal value.
Therefore, for the number of scenarioswe use S
ts
for service times, S
a
for earlieststart times, and
S
b
forlateststarttimes. Wealsoindicatethenominalcasescenariowithindex0. Lastly, weadapt
the same three uncertaintysets as before: convexhull, box, and ellipsoidal. The methodologyand
the proofs follow directly from Chapter 3. In addition, we address the probabilistic nature of the
customers in our robust formulation by using a simple method of introducing a probability cut.
5.3.1 Uncertainty in Service Times
To introduce uncertainty in service times, in Problem (4) we replace equation (4.7) with (4.16)
for convex hull; with (4.17) for box; and with (4.18) for ellipsoidal.
y
j
−y
i
−t
ij
+M
0
ts
(1−x
ij
)≥s
0
j
+max{max
k
s
k
j
,0} i,j ∈V \{0}, i6=j (4.16)
y
j
−y
i
−t
ij
+M
0
ts
(1−x
ij
)≥s
0
j
+
Sts
X
k=1
|s
k
j
| i,j ∈ V \{0}, i6=j (4.17)
79
y
j
−y
i
−t
ij
+M
0
ts
(1−x
ij
)≥s
0
j
+
v
u
u
t
Sts
X
k=1
(s
k
j
)
2
i,j ∈ V \{0}, i6=j (4.18)
In the above equations, the value of big M should be updated accordingly as in the following:
M
0
ts
=L+max
(i,j)
t
ij
+max
i
s
0
i
+max
j
δ
j
where δ
j
is equal to max{max
k
s
k
j
,0} for convex hull,
Sts
X
k=1
|s
k
j
|
for box, and
v
u
u
t
Sts
X
k=1
(s
k
j
)
2
for ellipsoidal.
5.3.2 Uncertainty in Time Windows
Similar to the uncertainty in service times, to introduce uncertainty in time windows to Problem
(4), we replace equations (4.8) and (4.9) with (4.19) and (4.20) for convex hull; with (4.21) and
(4.22) for box; and with (4.23) and (4.24) for ellipsoidal.
y
i
+e
i
+M
0
a
(1−z
i
)≥a
0
i
+max{max
k
a
k
i
,0} i∈ C (4.19)
y
i
−l
i
≤b
0
i
+min{min
k
b
k
i
,0} i∈C (4.20)
y
i
+e
i
+M
0
a
(1−z
i
)≥a
0
i
+
Sa
X
k=1
|a
k
i
| i∈ C (4.21)
y
i
−l
i
≤b
0
i
−
S
b
X
k=1
|b
k
i
| i∈ C (4.22)
y
i
+e
i
+M
0
a
(1−z
i
)≥a
0
i
+
v
u
u
t
Sa
X
k=1
(a
k
i
)
2
i∈ C (4.23)
80
y
i
−l
i
≤b
0
i
−
v
u
u
t
S
b
X
k=1
(b
k
i
)
2
i∈C (4.24)
In the above equations, the value of big M should be updated accordingly as in the following:
M
0
a
= max
i
(a
i
+δ
i
) where δ
i
is equal to max{max
k
a
k
i
,0} for convex hull,
Sa
X
k=1
|a
k
i
| for box, and
v
u
u
t
Sa
X
k=1
(a
k
i
)
2
for ellipsoidal.
5.3.3 Probabilistic Customers
In constructing robust routes, it may not be realistic nor feasible to consider all the potential
customers for largescale problem instances. To model the probabilistic nature ofthe customers in
a robust formulation, one idea could be considering only a subset of all potential customers. Let
O
d
bethendimensionaldatavectorindicating theoccurrenceofeachnodeonagivenscenario(or
day)d outofD totaldays. In this data, ifnodei appearsin dayd then O
d
i
= 1, otherwiseO
d
i
=0.
Then the general probability of occurrence of customer i is given by: p
i
=
P
d
O
d
i
D
. Now, using a
general probabilistic model, we can consider only customer with high probability of occurrences
defined by a cut ς. That is, customer i∈C if p
i
>ς.
5.4 Recourse Model for Robust Counterpart of VRP
In this section, we consider customers in Problem (4) to be probabilistic and we use scenario
based stochastic optimization to address this uncertainty. For a given planning horizon of D
days, we consider each day as a scenario and we assume that the data indicating customers that
appear in each day is available (i.e. O
d
i
). Accordingly we formulate stochastic recourse models
in which some of the variables are fixed at the end of the first stage (a priori variables) and the
others are determined in the second stage based on the fixed variables (adjusted variables). We
also construct daily schedules based on recourse actions such that the similarity of routes among
81
days is maximized. We eventually combine this stochastic optimization approach with the robust
optimization approach used in the previous section to formulate the courier delivery problem.
Althoughdifferenttypesofrecourseactionshavebeeninvestigatedintheliterature,thereisnot
much research done in incorporating the recourse action into the mathematical formulation with
additionalvariablesandconstraintsratherthanincorporatingitintotheprobabilisticcomputation
as an expected value or than using real-time re-optimization techniques. We propose in the
next two subsections VRP models with stochastic customers capturing the recourse action in the
formulation through scenarios. We also put everything together to construct the recourse model
for the robust counterpart of the courier deliver problem considered in this chapter.
5.4.1 Recourse Action of Omitting Customers
In this subsection, we consider a classical recourse action in the literature, omitting customers.
We incorporate this recourse action by augmenting Problem (4) with additional constraints and
variables. We treat the routes of Problem (4) as the master plan (first stage a priori routes) and
we create additional recourse routes for each day (second stage adjusted routes). We force the
definition of omitting customers in the following way. For a given scenario, we only service cus-
tomers which arealreadyserviced in the master plan. Thatis, the customerswhich do not appear
in the master plan but which are part of the particular day data are omitted in constructing that
daily schedule. The other occurring customers are visited in the daily schedule in the same order
as they are visited in the master plan. Therefore, this recourse action is very simple and avoids
rescheduling of daily routes. However, the resulting daily schedules are completely dependent on
the master plan and there is no flexibility to deviate from the master plan to schedule customers
which are part of that daily data but which are not scheduled in the master plan.
Let set N
D
be the set of integers from 1 to D. The binary r
d
ij
are recourse variables associated
with each day d and each arc (i,j) in the network to indicate the adjusted routes traversedby the
vehicles based on the a priori routes given by x
ij
variables. Then we need to relate the recourse
82
variables to a priori variables in the stochastic VRP based on D scenarios. We add the following
setofconstraintstoProblem(4)tospecifyexactlyallthearcstobe selectedin theadjusted route:
O
d
i
O
d
j
x
ij
≤r
d
ij
i,j ∈ V, i6=j, d∈N
D
where on a particular day d, all the arcs between the occurring nodes on the a priori route are
selected for the adjusted route;
O
d
i
(1−O
d
j
)x
ij
+(1−O
d
j
)O
d
k
x
jk
≤1+r
d
ik
i,k ∈V, j ∈V \{0}, i6=j 6=k, d∈N
D
where on a particular day d, if nodes i, j, and k are visited in the a priori route in this order but
node j does not occur, then the arc (i,k) is selected for the adjusted route.
In general, on a particular day d, if a sequence of nodes i
1
, i
2
, ..., i
q−1
, i
q
are visited in the a
priori route in this order but the inner sequence of nodes i
2
, ..., i
q−1
do not occur, then to force
selecting arc (i
1
,i
q
) in the adjusted route the below constraint should be added to Problem (4)
for i
1
,i
q
∈ V, i
2
,...,i
q−1
∈V \{0}, i
1
6=i
2
6=...6=i
q−1
6=i
q
, d∈ N
D
.
O
d
i1
(1−O
d
i2
)x
i1i2
+
q−2
X
j=2
(1−O
d
ij
)(1−O
d
ij+1
)x
ijij+1
+(1−O
d
iq−1
)O
d
iq
x
iq−1iq
≤q−2+r
d
i1iq
For a given day d, there are total of n
2
additional variables and
n−1
X
i=2
n!
(n−i)!
additional con-
straints. The total additional complexity due to the recourse is O(Dn
2
) variables and O(Dn!)
constraints. These additional constraints enforce lower bounds on the recourse variables depend-
ing on the occurrence data. When a constraint is active, then the related r
d
ij
is set to 1, meaning
that this arc (i,j) is part of the adjusted route at day d; otherwise the constraint is redundant
and there is no lower bound enforced (i.e. r
d
ij
∈{0,1}). The constraints which become active are
due to selecting the same arcs in the adjusted route as in the a priori route (between consecutive
occurring customers) and due to enforcing the definition of the recourse action by selecting arcs
83
in the adjusted route which are not part of the a priori route (between non-consecutive occurring
customers which have non-occurring customers in between that are to be skipped). For the latter
case, the lower bound is correctly enforced by considering all possible subsets of the nodes in the
network and by checking for each subset all possible combinations of skipping nodes to identify
which nodes to skip. Therefore these active constraints form valid adjusted routes and moreover
they exactly specify the totality of the r
d
ij
variables to be set to 1. Now all we need to do is to set
all the remaining r
d
ij
variables to 0 to prevent from selecting unnecessary arcs in the solution. For
this we need to add the following constraint to Problem (4) to specify the total number of arcs
that should be selected in the adjusted route based on the a priori route and the occurrences so
that for all the remaining arcs which should not be selected in the adjusted route, r
d
ij
is correctly
set to zero:
K +
X
j∈V\{0}
O
d
j
X
i∈V, i6=j
x
ij
=
X
i∈V
X
j∈V, j6=i
r
d
ij
d∈N
D
It is easy to see that the LHS of this equality is the sum of the number of vehicles and the number
of occurring customers on the a priori route, and the RHS is the number of arcs to be selected by
the adjusted route. If this number is m, then the above O(Dn!) constraints set exactly m of the
r
d
ij
variables to 1 (by forming K routes over m−K occurring customers) and the last constraint
imposes that a total of m r
d
ij
variables must be 1. These two conditions combined imply that all
the other r
d
ij
variables are set to 0. Here, the only implicit assumption coming from the initial
formulation to form the a priori routes is that all the K vehicles must be used to form exactly K
routes on any given day (i.e. on a given day at least one customer along every a priori route must
occur).
84
To maximize the similarity of the adjusted routes with the a priori routes, let’s consider the
total travel time deviation from the a priori route due to the recourse action. We then need to
further add the following set of constraints to Problem (4):
X
i∈V
X
j∈V, j6=i
t
ij
(r
d
ij
−x
ij
)≥γ d∈N
D
Here,γ isafreevariabletoaccountforthepenaltyfordeviation. Thisvalueisalwaysnon-positive
since the nodes in the networkobey the triangularinequality. Lastly, let α
5
be a positive constant
to prioritorize the importance of similarity of routes in the objective function and replace the
objective function in Problem (4) with the following:
max α
1
X
i∈C
p
i
z
i
−α
2
X
i∈A
y
i
−α
3
X
i∈C
l
i
−α
4
X
i∈C
e
i
+α
5
γ
Depending on the probability of occurrences p
i
the stochastic model picks as many customers
as possible by minimizing the route lengths and the total penalty (earliness and lateness) as long
as they can be feasibly serviced with respect to the common due date L and it constructs a priori
routes based on D scenarios in such a way that the maximum deviation required by a recourse
action (which is to skip non occurring customers on the route) is minimized, with the expense
of polynomially many additional variables and exponentially many additional constraints. Note
that constraints (4.7-4.9)are independent from the above modifications and therefore they can be
replaced with their robust counterparts. Also, to identify the customers to be addressed in these
robust a priori routes the probability cut should be set to ς = 0 allowing robustroutes to schedule
all the potential customers since the recourse action does not provide scheduling of non-occurring
customers in a given scenario.
Lastly, although the number of constraints in this stochastic formulation is exponential in n,
the actual number of active constraints for a certain scenario d is much smaller with a given
85
solution x, as we discussed before. This suggests that a natural way of solving this formulation
could be a constraint generation procedure.
5.4.2 Recourse Action of Rescheduling Routes
In this subsection, we consider a flexible recourse action which allows complete rescheduling of
daily routes to increase the similarity among the days. This recourse action is more sophisticated
and realistic than omitting customers since all the customers occurring on a particular day are
considered to be part of that daily schedule. However this requires more efforts in formulating
the problem. As a general methodology, based on Problem (4) we define a nominal VRP to
determine the a priori routes for customers with high probability of occurrence. We indicate the
variables, uncertain parameters, and related sets with superscript 0. Accordingly, the value of the
probabilitycutς forcustomersinC
0
becomesaparameterofthe model. Then foreachscenariod,
we formulate a similar VRP with variables, uncertain parameters, and related sets superscripted
withd. NotethatwedonotneedtodetermineaprobabilitycutforcustomersinC
d
sincetheyare
determined by the occurrence data O
d
i
. Lastly, we combine all the models and relate the recourse
variables to a priori variables. Thus, the size of the stochastic VRP is D times the size of the
nominal VRP. The nominal problem based on Problem (4) follows below in Problem (5).
86
max α
1
X
i∈C
0
z
0
i
−α
2
X
i∈A
0
y
0
i
−α
3
X
i∈C
0
l
0
i
−α
4
X
i∈C
0
e
0
i
(5.1)
s.t.
X
j∈V
0
, i6=j
x
0
ji
=z
0
i
i∈ C
0
(5.2)
X
j∈V
0
, j6=i
x
0
ij
=z
0
i
i∈ C
0
(5.3)
X
i∈V
0
\{0}
x
0
0i
=K (5.4)
X
i∈V
0
\{0}
x
0
i0
=K (5.5)
x
0
i0
= 1 i∈ A
0
(5.6)
y
0
i
+t
ij
+s
0
i
≤y
0
j
+M
ts
(1−x
0
ij
)
i∈V
0
j ∈ V
0
\{0}
i6=j
(5.7)
a
0
i
≤y
0
i
+e
0
i
+M
a
(1−z
0
i
) i∈ C
0
(5.8)
y
0
i
≤b
0
i
+l
0
i
i∈ C
0
(5.9)
y
0
i
≤Lz
0
i
i∈ C
0
(5.10)
e
0
i
≥0 i∈ C
0
(5.11)
l
0
i
≥ 0 i∈ C
0
(5.12)
x
0
ij
∈{0,1} i,j ∈ V
0
, i6=j (5.13)
y
0
i
≥0 i∈ V
0
(5.14)
z
0
i
∈{0,1} i∈ C
0
(5.15)
(5)
Let C
d
be the set of customers that occur in a given scenario d, that is customer i ∈ C
d
if
O
d
i
= 1. As before, let A
d
be the set of integers from|C
d
|+1 to|C
d
|+K. Similarly, let V
d
be the
set of nodes that occur in a given scenario d, that is V
d
= C
d
∪A
d
∪{0}. Let N
D
be the set of
87
integers from 1 to D, as before. Then the corresponding additional set of constraints to be added
in Problem (5) is given below:
X
j∈V
d
, i6=j
x
d
ji
=z
d
i
i∈C
d
, d∈ N
D
(5.16)
X
j∈V
d
, j6=i
x
d
ij
=z
d
i
i∈C
d
, d∈ N
D
(5.17)
X
i∈V
d
\{0}
x
d
0i
=K d∈ N
D
(5.18)
X
i∈V
d
\{0}
x
d
i0
=K d∈ N
D
(5.19)
x
d
i0
= 1 i∈A
d
, d∈ N
D
(5.20)
y
d
i
+t
ij
+s
d
i
≤y
d
j
+M
ts
(1−x
d
ij
)
i∈ V
d
, j ∈V
d
\{0}
i6=j, d∈ N
D
(5.21)
a
d
i
≤y
d
i
+e
d
i
+M
a
(1−z
d
i
) i∈C
d
, d∈ N
D
(5.22)
y
d
i
≤b
d
i
+l
d
i
i∈C
d
, d∈ N
D
(5.23)
y
d
i
≤Lz
d
i
i∈C
d
, d∈ N
D
(5.24)
e
d
i
≥ 0 i∈C
d
, d∈ N
D
(5.25)
l
d
i
≥0 i∈C
d
, d∈ N
D
(5.26)
x
d
ij
∈{0,1} i,j ∈V
d
, i6=j, d∈N
D
(5.27)
y
d
i
≥ 0 i∈V
d
, d∈N
D
(5.28)
z
d
i
∈{0,1} i∈C
d
, d∈ N
D
(5.29)
In order to increase the similarity of routes, we consider the method of rewarding common
arcs in scenarios. For this, let the binary c
d
ij
be additional recourse variables associated with each
day d and each arc (i,j) in the network to indicate the common arc (i,j) which is traversed both
in the scenario d and in the a priori route. Then to relate recourse variables to a priori variables
88
in order to identify the common arcs in scenarios with the a priori routes, we further add the
following set of constraints:
x
0
ij
+x
d
ij
≥ 2c
d
ij
i,j ∈ V
d
, j 6=i, d∈N
D
(5.30)
c
d
ij
∈{0,1} i,j ∈ V
d
, i6=j, d∈ N
D
(5.31)
Lastly, we replace the objective function in equation (5.1) with the following, where α
5
is a
positive constant to prioritorize the importance of common arcs in the objective function:
max
X
d∈ND∪{0}
(α
1
X
i∈C
d
z
d
i
−α
2
X
i∈A
d
y
d
i
−α
3
X
i∈C
d
l
d
i
−α
4
X
i∈C
d
e
d
i
)+α
5
X
d∈ND
X
i∈V
d
X
j∈V
d
, j6=i
c
d
ij
t
ij
(5.32)
For the master plan and for each scenario the above model maximizes the total number of
customers which can be feasibly serviced by minimizing the route lengths and the total penalty
(earliness and lateness)and it constructs the a priori routes and the adjusted routes in such a way
that for each scenario the total length of the parts of the adjusted routes which are common with
the a priori routes is maximized.
If we consider the problem given by (5.2-5.32), modify the constraints (5.7-5.9) for a given
uncertainty set for the robust optimization methodology as described in the previous section and
determine the probability cut ς for C
0
, then we will obtain the recourse model for the robust
counterpart of the courier delivery problem. In this formulation, we construct a master plan
based on uncertain service times and time windows due to a robust optimization approach and
we construct daily schedules based on probabilistic customers represented through scenarios due
to stochastic optimization approach and on a recourse action of complete rescheduling of daily
routes to increase the similarity among scenarios. In other words, we generate a robust master
plan which is used as a basis in generating flexible daily schedules using a recourse action of
89
rescheduling routes. This formulation combines the simplicity of robust optimization and the
flexibility of stochastic optimization with recourse.
5.5 Real Life Problem Data and Uncertainty
In this section, we present the real life problem data of the courier delivery problem in industry
which is faced by a worldwide delivery company. We discuss both the characteristics of the data
and investigate the nature of the uncertainty.
The problem is a large scale real life problem with 3715 potential customers in an urban area
of northwest California. The locations of the customers are known and most of them can be
enclosed in a box of 25 miles
2
as shown in Figure 5.1. At the beginning of each day, any of these
potential customers can put a delivery request with a service time and time windows (earliest
start time and latest start time). Therefore on a given day, customers occur probabilistically with
uncertain service times and time windows. The challenge is to construct daily schedules to satisfy
as many requests as possible for each day at a minimum cost while increasing the similarity in
daily schedules during the planning horizon and thus decreasing the efforts of daily re-planing.
It is desired to have an a priori master plan for the planning horizon to be used as a basis for
daily schedules both to increase the similarity in routes and to estimate the resource allocations
for that horizon. Since for a particular customer, the service time and time windows are known
but changesfor eachday, they are uncertainin the masterplan. Also, the costfor the masterplan
and the daily schedules are the lateness penalty and the total time spent by the couriers which is
the sum of travel and waiting times.
The data of delivery requests is collected for a horizon of 29 days which is also the planning
horizon. A typical day appears in Figure 5.2 depicting the occurring customers on that particular
day among all the potential customers. The average number of occurring customers per day is
471. There are a total of 4 couriers to service the customers each day. The couriers are ready to
leave the depot at 9am and they must return to the depot by 8pm. For each of the 29 days, the
90
40 40.2 40.4 40.6 40.8 41 41.2 41.4
−124.25
−124.2
−124.15
−124.1
−124.05
−124
−123.95
−123.9
−123.85
−123.8
−123.75
Longitude
Latitude
Figure 5.1: Location of potential customers
40.7 40.75 40.8 40.85 40.9 40.95 41
−124.25
−124.2
−124.15
−124.1
−124.05
Longitude
Latitude
Figure 5.2: A typical daily occurrence data of customers
91
[0−7.5) [7.5−15) [15−22.5) [22.5−30) [30−37.5) [37.5−45)
0
1000
2000
3000
4000
5000
6000
7000
8000
Interval
Number of Data Points
Data
Fit
Figure 5.3: Actual and fitted service time distributions
0 500 1000 1500 2000 2500 3000 3500
0
0.1
0.2
0.3
0.4
0.5
0.6
0.7
0.8
0.9
1
Customer ID
Probability of Occurrence
P1
P2
P3
P4
Figure 5.4: Actual and modified probability distributions of the occurrences of customers
92
realizations of service times and time windows of occurring customers are available. The travel
time is considered deterministic and to convert the distance measures to time units, it is assumed
that the couriers travel with an average speed of 35 mph in the city.
We adapt the formulation given by (5.2-5.32) to address the courier delivery problem. We
consider each day as a scenario in formulating the uncertainty. For the robust counterparts
of constraints (5.7-5.9), we use the averages of service times and time windows values of each
customer as its nominal values and calculate the deviations in each scenario from these averages.
This representation of uncertainty is then used in constructing robust routes in the master plan.
In addition, we have also analyzed the service time characteristics of all the potential customers
in general. Figure 5.3 shows the actual distribution of the data with bins together with fitted
distribution to this observed data. As it is usually the case in the routing literature and in
the industry, service times follow a lognormal distribution, see Dessouky et al. (1999). In our
particular case, we use a shifted and scaled lognormal distribution with mean μ
s
= 0.0953 and
standard deviation σ
s
= 0.25. The mean and the standard deviation of the actual sample data
are respectively about 4 and 7 minutes.
To generate the daily schedules, we use the robust routes in the master plan as a basis and
we adapt a hybrid recourse action of partial rescheduling of the routes for each day, which is a
combination of omitting non-occurring customers and complete rescheduling of routes. Here, a
given day of data (occurrences of customers, service times and time windows) constitutes again
a scenario. In particular, we have analyzed the probability distribution of the occurrences of all
the potential customers in general. In Figure 5.4, P2 is the actual discrete distribution and P4
is its continuous approximation which is a power function. Without loss of generality, customers
are represented by non-increasing IDs with respect to their probability of occurrence. Note that
there are customers with probability 1 (i.e. occurring in all of the scenarios). In addition, we
present two more discrete distributions, P1 and P3, which are generated by modifying P4. In
P1, the probabilities of occurrences are decreased with respect to original P2; and in P3, they are
93
increased. These three distributions, P1, P2 and P3, will beused later in our experimentalsection
in sampling scenarios.
5.6 Two-Phase Heuristic Solution Procedure
To address the challenge of solving the large scale real life problem, we develop a heuristic solu-
tion procedure based on insertion. Insertion has been proved to be a very efficient, simple, and
therefore widely used heuristic for large scale VRP instances with time windows. In this section,
we propose a problem specific solution procedure for our courier delivery problem in the industry.
In our two-phase heuristic, we develop concurrently a robust master plan and daily schedules for
scenarioswhich are constructed based on the master plan accordingto a hybrid recourseaction of
partial rescheduling of routes. This recourse action combines the two recourse actions of omitting
customers and rescheduling routes.
The first phase is the construction of initial solutions. First of all, a robust master plan is
constructed by making insertions of customers to initially empty routes by minimizing the cost
of insertion greedily. That is, among all the feasible insertions with respect to the common due
date the cheapest one is done at each step. The cost of an insertion of a customer to a location in
a vehicle is determined by the increase in the route length and the increase in the total lateness
penalty of all the customers serviced by that vehicle. Note that the cost of insertion is dependent
on the location of the inserted customer on the route. Second of all, the routes of the master plan
are updated by each scenario following the two recourse actions to construct the daily schedules.
First, a scenario modifies the initial robust routes of the master plan by following the recourse
actionofomittingcustomersandthenthereschedulingofthesemodifiedroutesisdonebyinserting
brand new customers in that particular day data by using again the greedy insertion.
The second phase is iterative. At each iteration, first, scenarios give feedback to the master
plan about the customers that could not be scheduled in their daily routes; second, based on this
feedbackthemasterplanprioritorizestheseunscheduledcustomersbythenumberofscenariosthat
94
they appear but could not be scheduled. Then the master plan updates its routes by performing
feasible maximum priority insertions in the cheapest way. Note that the selection of customers
is based on the priority not on the cost of insertion. However, once a customer is selected then
the cheapest possible insertion is done for this particular customer. Then the new master plan
is re-dispatched to the scenarios which construct their daily schedules with respect to the hybrid
recourse action as before. At the end of each iteration, the objective function is evaluated over all
the scenarios based on the four measures in the following order of priority: number of customer
serviced, total length of the routes, total lateness penalty, and similarity ofroutes with the master
plan. The iterationsofthe secondphaseandthus the overalltwo-phasealgorithmstop whenthere
is no improvement in the objective.
Note that for the master plan, the first phase is cost driven whereas the second phase is
priority driven. For the scenarios, both phases are cost driven based on the current master plan.
The maximization of the number of customers serviced is mainly due to the recourse action of
rescheduling of the routes; the minimization of the route lengths and lateness penalty is mainly
due to the greedy insertions; and maximization of the similarity of routes is due to the recourse
action of omitting customers and to the feedback procedure between the master plan and daily
schedules to prioritorize the customers in the iterative second phase.
Betweenthe twophasesofthe algorithm, wealsoimplementthe ideaofbuffer capacitystudied
in the previous chapter. For the first phase, we reserve a “buffer capacity” by decreasing the
common due date of vehicles. This slack time is later used in the second stage to schedule
additional customers. This parameter of the algorithm is actually a tool to balance the cost
driven stage and the priority driven stage in constructing the master plan which indirectly effects
the daily schedules as well. We later explore experimentally the effect of this parameter on the
aspects of cost and similarity of routes.
Another parameter of the algorithm is the set of customers to be considered in the master
plan during the first phase. Recall that the average number of customers per day is 471 whereas
total number of potential customers is 3715. Attempting to schedule all the customers is neither
95
realistic nor feasible. Instead, a subset of customers should be selected. One intuitive idea is
to select some of the customers which appear the most in the scenarios, customer with a high
probability of occurrence. This should increase the similarity of routes since the master plan will
be composed of mainly the customers which will not be omitted in most of the daily schedules.
A reasonable and realistic number to pick for the cardinality of the set of customers with high
probability of occurrences could be a value around the average number of customers per day. We
also explore later the effect of this parameter on the quality of the solution.
The last parameters of the algorithm are the weights used to balance the cost of increase in
route length and the cost of increase in lateness penalty during the greedy insertion. Changing
these weights in the selection of cheapest insertion will alter the outcome by emphasizing more
on minimizing either the route length or the lateness penalty. However, to keep the number of
control parameters small we will weight these two measures equally.
Thealgorithmicdescriptionofthetwo-phasesolutionprocedureisgiveninAlgorithm1-3. Note
that as a result of the second phase the master plan is generated based on historical scenarios.
Then the master plan can be used during the planning horizon for each future scenario in the
following way: first non-occurring customers in the data of a particular scenario are dropped and
then Algorithm 1 is executed with these initial routes in order to generate the daily schedule.
This is repeated in each day during the planning horizon to generate daily schedules with respect
to the master plan. The implicit assumption here is that the new realizations of scenarios will
be similar to the past data, i.e. there is no seasonal trends or irregular outcomes. Also, at the
end of each planning horizon, based on new observations of that particular planning horizon the
two-phase algorithm can be used in generating a new master plan for the next planning horizon.
96
Algorithm 1 Insertion Routine
Require: Initial routes
Calculateinsertion costofpossible insertionlocationsin the routes fornon-scheduledcustomers
repeat
Pick the cheapest insertion (a feasible insertion of a non-scheduled customer which will result
in the leastincreaseofthe total length ofthe routes andthe total latenesspenalty)and insert
it into the vehicle at the proper location
Update insertion cost of possible insertion locations of that vehicle
until All the occurring customers are inserted or no feasible insertion is possible
return The resulting routes
Algorithm 2 Phase One
Require: Distance matrix of the network, scenario data for each day (customers appearing in
thatdaytogetherwithearlieststart,lateststart,andservicetimedata),masterdata(customers
to be scheduled in the master plan during Phase One), percent buffer capacity
Reserve a buffer capacity in the length of the routes
Setup the masterplandata(bycalculatingworstcasevaluesofservicetimesand time windows)
Call Insertion Routine for the master plan with initially empty routes
for Each scenario i do
Take the master routes as initial solution and drop the non-occurring customers in the data
of scenario i
Call Insertion Routine for scenario i with these initial routes
end for
Calculate the objective value and save the current solution
return The current solution
97
Algorithm 3 Phase Two
Require: Solution of Phase One as the initial solution
Remove reserved buffer capacity in the length of the routes from consideration
repeat
Calculate priorities of customers which could not be inserted in the scenarios and create a
sorted list with non-increasing priorities
if The priority list is empty then
return The current solution
end if
for Each customer i in the priority list do
if Customer i can be feasibly inserted then
Make the least cost insertion of customer i in the master schedule
end if
end for
if No feasible insertion was possible then
return The current solution
end if
for Each scenario i do
Takethe masterroutesas initial solutionand dropthe non-occurringcustomersin the data
of scenario i
Call Insertion Routine for scenario i with these initial routes
end for
Calculate the objective value and save the current solution
until The objective value is not improved
return The current solution
98
5.7 Experimental Results
In this experimental section, we consider the courier delivery problem developed in this chapter
adapted from (5.2-5.32) and we use the two-phase solution procedure described above to obtain
efficientapproximatesolutions. Inparticular,wereplacetheconstraints(5.7-5.9)withtheirrobust
counterparts for the convex hull uncertainty as in (4.16) for uncertainty in service times and as in
(4.19-4.20) for uncertainty in time windows. We try different values for the probability cut ς to
determine C
0
. Also, in the objective function we set the parameterα
4
= 0 since our problem does
not consider earliness penalty. For the other 4 goals, we prioritorize α
1
as the highest, we weight
α
2
and α
3
equally but put less priority than α
1
, and we give α
5
the smallest priority. Recall that
to increase the similarity of routes in scenarios during our solution procedure, we use the hybrid
recourse action of partial rescheduling of routes, i.e. first omitting non-occurring customers and
then rescheduling routes for the brand new customers in the scenarios.
We analyze the sensitivity of the two-phase algorithm to the problem data, effect of changing
the service time distribution and effect of changing the probability distribution of the occurrences
of customers, as well as to its parameters, effect of changing the set of customers to be scheduled
in the master plan during the first phase and effect of changing the buffer capacity. Recall that
we set the two weights to be used in calculating the cost of insertion equal to each other and we
fix them throughout the experiments to keep the number of control parameters small.
Wealsocomparethequalityofthesolutionofthetwo-phasealgorithmwithasolutionobtained
by independently scheduling scenarios without a master plan which we refer as the independent
daily insertion algorithm. This solution is obtained by treating each scenario as an independent
courierdeliveryproblemwithnoscenariosandzerobuffercapacity. Thismeansthateachscenario
in the original problem is considered as a “master plan” and only the first stage of phase one of
the two-phase algorithm is executed without reserving any buffer capacity to construct the initial
“master plan” with greedy insertions, which results in the daily schedule of the scenario in the
original problem. Therefore this algorithm does not make any considerations in increasing the
99
similarity of the resulting independent daily schedules but maximizes the number of customers
serviced while minimizing total route length and total lateness penalty. We also evaluate the
solution of the two-phase algorithm on the real problem data and show that our solution can be
better than the current practice.
Throughout the experimental analysis, we assume that the available past data from the pre-
vious planning horizon stays the same for the current planning horizon. That is, this data is used
both in generating and evaluating the performance of the solution of the two-phase algorithm.
In particular for the real-life instance, this is inevitable due to lack of replication. Lastly, all the
experiments are performed on a Dell Precision 670computer with a 3.2 GHz Intel Xeon Processor
and 2 GB RAM running Red Hat Linux 9.0.
5.7.1 Two-Phase versus Independent Daily Insertion
In this first set of experiments, we explore the effect of changing problem data on the two-phase
algorithm and we compare the quality of the solutions with the independent daily insertion algo-
rithm, leavingthe effectofchangingtheparametersofourheuristictothe nextsetofexperiments.
Thereforeforthecurrentexperimentsweneedtofixthevalueoftheseparametersofthetwo-phase
algorithm, the buffer capacity and the set of customers to be scheduled in the master plan during
phase one. For the latter, we refer to the distribution P2 in Figure 5.4 and we define the value of
thecutς toselectthecustomerswithhighprobabilityofoccurrence. Inordertokeeptheresulting
number of total customers around the average number of customer per day, 471, we fix the value
of this cut to 0.25. That is, only the customers which have a higher probability of occurrence
than 0.25 are considered in scheduling the master plan during phase one. This total number of
customers turns out to be 480 for the real problem data. For the former problem parameter, we
setthe buffer capacityto 50%, meaning that only halfof the total allowedtime for the tour length
is considered during the first phase. This way, we weight the cost driven first phase and priority
driven second phase equally.
100
Toexploretheeffectofuncertaintyoftheproblemparameters,wefocusonlyontheuncertainty
in customers and service times keeping the values of time windows of customers at their nominal
values. The reason is twofold. First, some of the customers in the problem data have fix time
windows which are not subject to uncertainty and some of them have very little deviations,
which suggest that the effect of uncertainty in time windows is less significant compared to other
uncertainties in the problem. Second, to keep the number of control parameters small for the
experiments, we prefer removing uncertainty in time windows from consideration in our first set
of experiments.
For the uncertainty in service time and probabilistic customers, we consider a base case with
respect to the fitted lognormal distribution in Figure 5.3 and the probability distribution P2 in
Figure 5.4. For each problem instance, we sample data of scenarios with respect to these two
distributions. First, the occurrence data of each scenario is generated by randomly selecting
customers according to P2 until 471 customers are selected in each scenario. Then each customer
is assigned a random service time following the lognormal distribution. Recall that the time
windows have been already fixed for each customer. Thus, all the required data for a problem
instance is generated by this process.
We generate additional cases by deviating from the base case in two ways. First, we change
the standard deviationσ
s
ofthe lognormaldistribution to see the effect of increasedservice times,
with σ
s
= 0.500, and decreased service times, with σ
s
= 0.125. Second, instead of P2 we sample
customers from P1 and P3 in Figure 5.4. The effect of sampling from P1 with respect to P2 is
the following: in each scenario most of the customers are selected among the customers with high
probability of occurrences since the others have a significantly less chance to occur. Note that
sincethecustomersconsideredinthemasterplanduringthefirstphaseofthetwo-phasealgorithm
also have high probability of occurrences, the similarity in routes of this heuristic solution should
improve on such instances. The effect of sampling from P3 is simply the reverse: this time the
customers which have relatively smaller probability of occurrences have in fact a higher chance to
101
occur in each scenario with respect to P2. Thus, in the solution of the two-phase heuristic, the
similarity measure is expected to be low for such instances.
For each case, we generate 10 random problem instances and report the averages of the solu-
tions. When moving from one case to another we modify only one parameter at a time keeping
the rest of the problem instance the same, which allows observing the sole effect of changing this
particular parameter. In Table 5.1, we provide these experimental results: left part is the input
parameters and the right part is the output measures. The headings are as follows: “Alg” is the
algorithm used in generating solutions, “TP” is for two-phase and “II” is for independent daily
insertion; “Prob” is the probability distribution in Figure 5.4 used to sample customers, 1 for P1,
2 for P2, and 3 for P3; “Std” is the value of the standard deviation of lognormal distribution of
service times, σ
s
; “NS” is the total number of customers which could not be serviced in the daily
schedules; “Length” is the total route lengths of daily schedules; “Penalty” is the total lateness
penalty in daily schedules; “Arcs” is the total length of common arcs in daily schedules; and
“Time” is the solution time in seconds. Since we compare the quality of the solutions with the
independent daily insertion algorithm which does not generate a master plan, we slightly mod-
ify the similarity measure “Arcs” to make the solutions comparable. Instead of calculating the
similarity of a scenario with respect to the master plan, we calculate this number based on the
similarity of scenarios with each other.
The general observations based on Table 5.1 suggest that the independent daily insertion
covers customers with smaller route lengths but with an increased lateness penalty and decreased
similarity of routes compared to the two-phase. Therefore, generally speaking the two-phase
improves the similarity of routes and lateness penalty at the expense of increased route length.
The improvement on similarity is expected due to the feedback process in the two-phase between
the master plan and the scenarios. However, the improvement on the lateness penalty is mainly
due to the buffer capacity. During the greedy insertion routine, most of the lateness penalty
is incurred at the latest insertions in the routes since the cheapest ones are performed in the
beginning delaying the costly ones whose penalty increases even more due to previous cheaper
102
Table 5.1: Comparison of two-phase with independent daily insertion
Alg Prob Std NS Length Penalty Arcs Time
TP 1 0.125 0 200526.8 84.5 1029.2 796.3
TP 1 0.250 0 206824.6 190.3 919.7 909.0
TP 1 0.500 42.1 226855.3 1327.9 593.5 1024.1
TP 2 0.125 0 200361.8 141.7 594.0 840.2
TP 2 0.250 0 207336.7 274.2 613.7 928.1
TP 2 0.500 43.0 226928.0 1194.9 401.8 1059.4
TP 3 0.125 0 201682.4 213.9 114.2 923.3
TP 3 0.250 0 207239.0 333.6 98.9 996.7
TP 3 0.500 39.0 226538.9 1088.7 77.1 1115.3
II 1 0.125 0 193343.4 10941.6 37.5 2402.8
II 1 0.250 0 203853.2 6088.3 40.7 3271.0
II 1 0.500 7.4 222818.0 278.0 35.3 2117.1
II 2 0.125 0 193401.0 9479.1 29.5 2393.6
II 2 0.250 0 203198.2 5349.8 27.2 3279.5
II 2 0.500 8.1 222843.8 336.3 23.5 2134.3
II 3 0.125 0 193568.9 4761.1 16.4 2393.3
II 3 0.250 0 203008.8 2994.0 17.0 3278.1
II 3 0.500 7.8 222577.2 344.0 9.5 1980.1
103
insertions. Reserving a buffer capacity prevents the insertion routine from incurring these high
penalties during phase one; and during the iterative phase two since some of the insertions are
done based on priorities, again the greedy nature of the insertion routine is avoided up to some
extend. Another general comment is that high values of σ
s
result in customers which could not
be serviced in the solution. In such cases, independent daily insertion performs better in covering
customers than two-phase because there is no effort done in creating common arcs in scenarios,
which provides flexible schedules to service more customers. Lastly, independent daily insertion
takes longer computational time to find a solution since two-phase benefits from the master plan
in constructing daily schedules, which decreases efforts of re-planing in each scenario.
WhenweanalyzeTable5.1inparticularforthetwo-phase,weseethatingeneralincreasingthe
probability of occurrence for a given standard deviation increases route length, lateness penalty,
solution time, and decreases similarity, making all these measures worse. The result is intuitive
and expected since sampling from a wider likely range of customers results in diverse scenarios,
which increases the cost aspects and decreases the similarity. On the other hand, increasing
the standard deviation for a given probability of occurrences has also the same effect on these
measures as before with the addition that high values result in unserviced customers. This is also
expectedbecauseincreasedanddispersedservicetimesmaketheproblemworsewithrespecttoall
measures. As a generalconclusion, for the two-phasealgorithmit is desired to havea combination
of minimum values of σ
s
and “Prob” in the problem data.
Lastly, when we analyze Table 5.1 in particular for the independent daily insertion algorithm,
the results are more obscure and there is no overall general trend. However, increasing the prob-
ability of occurrence for a given standard deviation increases route length, but decreases lateness
penalty; and increasing the standard deviation for a given probability of occurrence has the same
effect with the addition that high values result in unserviced customers. This decrease in lateness
penalty with respect to increased standard deviation is mainly due to the load of customers in
the vehicles. High standard deviation results in high service times yielding more balanced routes
because of the tie-breaking and the nature of the cheapest insertion during the independent daily
104
insertion algorithm. However, the effect of low standard deviation is just the opposite. This is
intuitive since under high standard deviation, inserting several customers with high service times
in one vehicle is not feasible and therefore these customers are evenly distributed to all the vehi-
cles. On the other hand, in case of low standard deviation, several customers with small service
times can be fit in a single vehicle leaving only few other customers with small service times to
be inserted into the other vehicles. As a result, when the vehicles are balanced in case of high
standard deviation, it is easier to meet the deadlines (and hence to incur smaller lateness penalty)
at the expense of increased route length compared to the case of low standard deviation having
an uneven balance of customers in the vehicles.
5.7.2 Two-Phase versus Real Life Solution
In this second set of experiments, we both explore the effect of changing parameters of the two-
phase algorithm and comparethe quality ofits solution with the currentpractice. For the percent
buffercapacity,weexplorethefollowingrange: 0, 20,40,60,80,100. Notethat0%buffercapacity
bypasses the first phase of the heuristic by making insertions in the master plan only based on
priorities; and 100% is the reverse making all the insertions based on costs. For the cut ς of the
probability distribution P2 in Figure 5.4, we use 0.25 as the base case resulting in 480 customers
and consider the following deviations: 0.20 resulting in 563 customers and 0.30 resulting in 384
customers.
Table 5.2 shows the solutions on the real life problem instance for different parameter settings
of our heuristic. The new headings are “Cut” for the cut of the probability distribution and “BF”
for the percent buffer capacity. Also, “NA” stands for not applicable. We omit the column for
unserviced customers since in this particular real life problem all the customers could be feasibly
serviced in each solution. In addition, we provide the solution obtained by the independent daily
insertionalgorithmandthereallifesolution, aswellasaverageresultsofeachcutoverthe6buffer
capacity values for the two-phase algorithm indicated by “Avg” under the “BF” column. Also we
105
Table 5.2: Comparison of two-phase with the real life solution
Alg Cut BF Length Penalty Arcs Time
Real NA NA 212453.3 13964.6 692.0 NA
II NA NA 206364.0 1697.0 102.4 3920.6
TP 0.20 0 212193.1 8081.4 599.6 535.8
TP 0.20 20 212533.1 15818.6 713.9 534.9
TP 0.20 40 213238.5 11572.0 757.4 1010.8
TP 0.20 60 212803.5 12338.7 1614.1 1022.9
TP 0.20 80 213397.1 11688.6 1826.6 862.6
TP 0.20 100 213836.7 17357.6 2329.7 624.2
TP 0.20 Avg 213000.3 12809.5 1306.9 765.2
TP 0.25 0 212626.4 10643.8 796.1 530.6
TP 0.25 20 212610.4 22857.2 964.8 1030.4
TP 0.25 40 212281.1 8854.2 1196.9 1076.4
TP 0.25 60 212376.8 11747.6 1662.9 989.7
TP 0.25 80 212649.8 10072.9 1530.5 732.9
TP 0.25 100 213836.7 17357.6 2329.7 521.8
TP 0.25 Avg 212730.2 13588.9 1413.5 813.6
TP 0.30 0 212482.5 8095.5 1288.6 510.9
TP 0.30 20 213099.4 16979.1 1004.2 980.8
TP 0.30 40 212870.1 14662.4 1507.1 1100.7
TP 0.30 60 212446.5 12495.7 1513.8 972.1
TP 0.30 80 212715.6 10287.7 1793.8 753.6
TP 0.30 100 213836.7 17357.6 2329.7 523.2
TP 0.30 Avg 212908.5 13313 1572.9 806.9
use as before the slightly modified version of the similarity measure “Arcs” to make the solutions
comparable.
When we look at the individual solutions, we see that both the real life solution and the
two-phase are able to provide a solution with improved similarity at the expense of increased
route length and lateness penalty with respect to independently scheduling scenarios as in the
independent daily insertion algorithm. When we compare the quality of the two-phase solution
with the current practice, we see that in general the two-phase provides a better improvement
in similarity with a higher route length and/or lateness penalty. However, we observe that in
106
the following three cases the two-phase outperforms the real life solution in all of the measures:
Cut=0.25, BF=40; Cut=0.25, BF=60; and Cut=0.3, BF=60. This suggest that our heuristic can
be tuned to provide improvements over the current practice. Moreover, the first two improved
solutions are obtained around the base case values of these two parameters.
When we look at the average results for the two-phase algorithm, we see that as the cut
increasesthesimilaritymeasureisslightlyimprovedwhereastheeffectonroutelengthandlateness
penalty is not clear. The reason in the increase of similarity is due to the fact that the customers
with relatively less probability of occurrences are eliminated by increasing the cut, resulting in
more common arcs in the master plan with the scenarios. When it comes to the buffer capacity,
there is no general trend observed. However for a given value of cut, increasing buffer capacity
increasessimilarityofroutesin generalsincethis would putmoreemphasize onthe prioritydriven
second phase of the heuristic improving the similarity measure. Lastly note that for a given value
of cut, 0% buffer capacity behaves just like the indepenent daily insertion algorithm resulting in
the smallest route length and lateness penalty with the worst similarity of routes in general; on
the other hand, 100% is just the opposite, resulting in the highest similarity of routes with the
worst route length and lateness penalty in general.
5.8 Summary
In this chapter, we considered a real life courierdelivery problem, a variantof VRPTW, for which
we proposed a recourse model for the robust counterpart of the problem formulation and we
developed an efficient two-phase heuristic based on insertion. We addressed the uncertainty in
service times and time windows by using robust optimization and the probabilistic nature of the
customers by using scenariobased stochastic optimization with recourse. Thus, we benefited from
the simplicity of the robust model and the flexibility of recourse actions.
We first adapted an a priori robust model for the courier delivery problem. We then defined a
problem specific recourse action of partial rescheduling of routes which is a hybrid of two known
107
recourse actions in the literature: omitting non-occurring customers and complete rescheduling
of routes. We showed how to develop recourse models by incorporating each of these recourse
actions into the problem formulation to provide a master plan for the planning horizon as well as
daily schedules for scenarios which are modified versions of the master plan based on the recourse
action. For each scenario, our overall courier delivery problem formulation maximizes number of
customers serviced while minimizing route length and lateness penalty by increasing similarity of
routes with the robust routes of the master plan. For this model with hybrid recourse action of
partial rescheduling of routes, we eventually developed a two-phase heuristic using insertion as a
subroutine to solve efficiently the large scale real life problem.
We exploredexperimentallythe sensitivity ofourheuristic touncertain problemparametersas
well as to its own parameters. In particular, we showed that it is desired to have a small standard
deviation for the service time distribution of customers and small probability of occurrences for
the less frequent customers in scenarios. We also compared the quality of the solution with an
independent daily insertion algorithm which does not make efforts in providing a master plan
and thus in increasing similarity of routes. We observed that the two-phase heuristic improves in
general the similarity measure and the lateness penalty at the expense of increased route length.
We also compared the solution of our heuristic with the real life problem solution. We showed
that by tuning the parameters of our heuristic, it is possible to outperform the current practice in
all of the measures.
108
Chapter 6
Conclusions and Future Research Directions
In this study, we considered the robust counterpart of the VRP. We followed the robust optimiza-
tion methodology in formulating uncertainty in several different problem parameters. We had to
overcome the main challenges of solving a RVRP by identifying a suitable VRP formulation, by
formulatingan amenableRVRP with respectto meaningful uncertaintysets, and by adapting effi-
cientsolution procedures. This processresulted in an a prioriRVRP formulationwith uncertainty
in demand and traveltimes which isnotmoredifficult to solvethan the deterministic counterpart.
We extended our a priori RVRP formulation by incorporating scenario based stochastic opti-
mization with recourse for a real life courier delivery problem with uncertainty in service times
and time windowstogetherwith probabilisticcustomers. We had to overcomethe mainchallenges
of solving the recourse model of this RVRP by identifying a suitable a priori RVRP formulation,
by defining a realistic and practical recourse action of partial rescheduling of routes for the prob-
abilistic customers which is a combination of omitting non-occurring customers and complete
rescheduling of routes, by formulating an efficient recourse model which provides a robust mas-
ter plan for the planning horizon and daily schedules which are modified versions of the robust
master plan in each scenario based on the hybrid recourse action, and by developing a two-phase
insertion based heuristic to solve the large scale real life problem efficiently. This process resulted
in a recourse model which benefits from the simplicity of robust optimization and the flexibility
109
of scenario based stochastic optimization, as well as in solutions which are better than the cur-
rent practice with respect to the following four measures in each scenario: number of customers
serviced, length of the routes, lateness penalty, and similarity of routes.
It is always possible to extend this study by considering several other VRP formulations and
uncertainty sets in formulating robust counterparts, in particular by considering dependency in
uncertaintysets. Forinstance,anincreaseintraveltimeuncertaintymayresultinaccumulationof
demand atthe customerlocationwhich in returnmay increasethe demand uncertainty. As shown
in this study, developing and solving different RVRP will require addressing the main challenges
of robust optimization.
AninterestingfutureresearchdirectioncouldbereconsideringRVRPwithdemanduncertainty
for which robust optimization implicitly assumes the feasibility of each scenario. In case scenarios
result in infeasibility for the robust counterpart, a different model and a different approach will
be necessary to address this problem. One idea could be constructing robust routes minimizing
the unmet demand in these problematic scenarios.
Weshowedthatrobustoptimizationcouldbeacompetitiveapproachinaddressinguncertainty
for the VRP when compared to the deterministic solution. Our experimental analysis was based
on comparing the quality of the robust and deterministic solutions under the worst case and the
expected case scenarios by using performance measures. Another experimental approach could
be evaluating the quality of these solutions under all (or some) of the possible realizations of the
uncertainty by simulating the uncertain outcomes for fixed solutions.
In addition, throughoutour experiments we mainly focused on investigating the effect of prob-
lem parameters on the robust solution. Another approach could be investigating the effect of the
uncertainty set. It is clear that as the size of the uncertainty set increases, the robust solution
becomes more conservative. However, this increase gradually changes the quality of the robust
solution. One interesting idea could be identifying the minimum increase in the uncertainty which
would make the robust solution preferable compared to deterministic solution, or similarly iden-
tifying the maximum increase which would still keep the robust solution preferable.
110
We first contrasted robust optimization and scenario based stochastic optimization and then
wecombinedboth methods forourreallife courierdeliveryproblem. Thus, weshowedthatrobust
optimizationcanbeusedforbothaprioriandrecoursemodels. Forthisparticularcourierdelivery
problem, we showed how to generate a robust master plan which would increase the similarity
of daily schedules. One idea could be using different approaches in generating the master plan
and to compare the quality of the solutions. For instance, using an expected value approach
for the uncertainty in master plan may not be a viable option in terms of similarity but other
probabilistic approaches can be used such as deviating from the expected value by a fraction of
the standard deviation of the uncertain parameter instead of considering the worst case value of
robust optimization. Depending on the particular application, such a method might compete well
with robust optimization.
Lastly, robust optimization and stochastic optimization can be adapted to address the same
uncertainty in the model as alternative approaches to each other. However, this will require
different setof assumptions to follow, which would result in differentdefinitions of the uncertainty
set. An interesting idea could be exploring the compatibility of these approaches to formulate
comparableuncertainty sets to measurethe quality oftheir solutions. As shown in this study, this
will require addressing the main challenges of solving a RVRP as well as the main challenges of
solving a recourse model, which in addition requires the compatibility of these two formulations
with respect to the set of assumptions and to the definitions of the uncertainty set.
111
References
Achuthan, N. R., L. Caccetta, and H. S. P. (1996). A new subtour elimination constraint for the
vehicle routing problem. European Journal of Operational Research 91, 573–586.
Alfa, A. S., S. S. Heragu, and M. Chen (1991). A 3-opt based simulated annealing algorithm for
vehicle routing problems. Computers & Industrial Engineering 21, 635–639.
Applegate, D., W. Cook, S. Dash, and A. Rohe (2002). Solution of a min-max vehicle routing
problem. INFORMS Journal on Computing 14(2).
Arora,S. (1996). Polynomialtime approximationscheme forEuclideanTSPand othergeometric
problems. In 37th Ann. IEEE Symp. on Foundations of Comput. Sci., pp. 2–11.IEEEComputer
Society.
Augerat, P., J. Belenguer, E. Benavent, A. Corbern, D. Naddef, and G. Rinaldi (1995). Compu-
tationalresultswith abranchandcutcodeforthe capacitatedvehicleroutingproblem. Research
report 949-m, Universite Joseph Fourier, Grenoble, France.
Averbakh, I. (2001). On the complexity of a class of combinatorial optimization problems with
uncertainty. Mathematical Programming 90, 263–272.
Averbakh, I. and O. Berman (2000). Algorithms for the robust 1-center problem on a tree.
European Journal of Operational Research 123(2), 292–302.
Baldacci, R., E. Hadjiconstantinou, and A. Mingozzi (2004). An exact algorithm for the capaci-
tated vehicle routing problem based on a two-commodity network flow formulation. Operations
Research 52(5), 723–738.
Balinski, M. and R. Quandt (1964). On an integer program for a delivery problem. Operations
Research 12, 300–304.
Beasley,J.E.andN.Christofides(1997).Vehicleroutingwithasparsefeasibilitygraph.European
Journal of Operational Research (98), 499–511.
Ben-Tal, A., L. El-Ghaoui, andA. Nemirovski(2000, March). Robustsemidefinite programming.
InR.Saigal,L.Vandenberghe,andH.Wolkowicz(Eds.),Handbook of Semidefinite Programming.
Waterloo: Kluwer Academic Publishers.
Ben-Tal, A., A. Goryashko, E. Guslitzer, and A. Nemirovski (2003). Adjustable robust solu-
tions of uncertain linear programs. Technical Report, Minerva Optimization Center, Technion.
http://iew3.technion.ac.il/Labs/Opt/index.php?4.
112
Ben-Tal, A., T. Margalit, and A. Nemirovski (2000). Robust modeling of multi-stage portfolio
problems. In K. Roos, T. Terlaky, and S. Zhang (Eds.), High Performance Optimization, pp.
303–328. Dordrecht: Kluwer Academic Publisher.
Ben-Tal, A. and A. Nemirovski (1998). Robust convex optimization. Mathematics of Operations
Research 23(4), 769–805.
Ben-Tal, A. and A. Nemirovski (1999). Robust solutions to uncertain programs. Operations
Research Letters 25, 1–13.
Ben-Tal, A. and A. Nemirovski (2000a). Robust solutions of linear programming problems con-
taminated with uncertain data. Mathematical Programming 88(3), 411–424.
Ben-Tal, A. and A. Nemirovski(2000b, March). Structural design via semidefinite programming.
InR.Saigal,L.Vandenberghe,andH.Wolkowicz(Eds.),Handbook of Semidefinite Programming.
Waterloo, Canada: Kluwer Academic Publishers.
Bent, R. W. and P. V. Hentenryck (2004). Scenario-basedplanning for partially dynamic vehicle
routing with stochastic customers. Operations Research 52(6), 977–987.
Bertsimas,D.(1988). Probabilistic combinational optimization problems. Ph.D.thesis,Operation
Research Center, Massachusetts Institute of Technology, Cambridge, MA.
Bertsimas, D. (1992). A vehicle routing problem with stochastic demand. Operations Re-
search 40(3), 574–585.
Bertsimas, D., P. Jaillet, and A. R. Odoni (1990). A priori optimization. Operations Re-
search 38(6), 1019–1033.
Bertsimas,D. andM. Sim (2003). Robustdiscreteoptimizationandnetworkflows. Mathematical
Programming 98(1-3), 49–71.
Bertsimas, D. and M. Sim (2004). The price of robustness. Operations Research 52(1), 35–53.
Bertsimas, D. and D. Simchi-Levi (1996). A new generation of vehicle routing research: robust
algorithms, addressing uncertainty. Operations Research 44(2), 286–304.
Bodin, L., B. Golden, A. Assad, and M. Ball (1983). Routing and scheduling of vehicles and
crews: the state of the art. Computers & Operations Research 10(2), 63–211.
Branke,J.,M.Middendorf,G.Noeth,andM.M.Dessouky(2005).Waitingstrategiesfordynamic
vehicle control. Transportation Science 39, 298–312.
Carraway, R., T. Morin, and H. Moskowitz (1989). Generalized dynamic programming for
stochastic combinatorial optimization. Operations Research 37(5), 819–829.
Caseau, Y. and F. Laburthe (1999). Heuristics for large constrained vehicle routing problems.
Heuristics 5, 281–303.
Charikar,M.,S.Khuller,andB.Raghavachari(2001). Algortihmsforcapacitatedvehiclerouting.
SIAM J. Comput. 31(3), 665–682.
Cheeseman, F., B. Kanefsky, and W. M. Taylor (1991). Where the really hard problems are. In
Proceedings NCAI-91, Sydney, pp. 331–337.
113
Christofides, Mingozzi, and Toth (1979). The vehicle routing problem. In Christofides, Mingozzi,
Toth, and Sandi (Eds.), Combinatorial Optimization, pp. 315–338. Chichester: Wiley.
Christofides, Mingozzi, and Toth. (1981a). Exact algorithms for the vehicle routing problem,
based on spanning tree and shortest path relaxations. Mathematical Programming 20(3), 255–
282.
Christofides,Mingozzi,andToth.(1981b). Statespacerelaxationproceduresforthecomputation
of bounds to routing problems. Networks 11, 145–164.
Christofides,N.(1976).Worst-caseanalysisofanewheuristicforthetravellingsalesmanproblem.
TechnicalReport388,GraduateSchoolofIndustrialAdministration,Carnegie-MellonUniversity,
Pittsburgh.
Clarke, G. and J. Wright (1964). Scheduling of vehicles from a central depot to a number of
delivery points. Operations Research 12(4), 568–581.
Cordeau, J. F., M. Gendreau, G. Laporte, J. Y. Potvin, and F. Semet (2002). A guide to vehicle
routing heuristics. Journal of the Operational Research Society 53(5), 512–522.
Daganzo, C. F. (1978). An approximate analytic model of many-to-many demand responsive
transportation systems. Transportation Research 12, 325–333.
Daganzo, C. F. (1984). The distance traveled to visit n points with a maximum of c stops per
vehicle: an analytical model and an application. Transportation Science 18, 331–350.
Daganzo, C. F. (1988). Shipment composition enhancement at a consolidation center. Trans-
portation Research 22B(2), 103–124.
Dantzig, G. B. (1955). Linear programming under uncertainty. Management Science, 197–206.
Dantzig,G. B.andJ.H. Ramser(1959). The truckdispatchingproblem. Management Science 6,
80–91.
De Backer, B., V. Furnon, P. Kilby, P. Prosser, and P. Shaw (2000). Solving vehicle routing
problems using constraint programming and metaheuristics. Heuristicsi 6, 501–523.
Desrochers, M. and G. Laporte (1991). Improvements and extensions to the millertucker-zemlin
sub-tour elimination constraints. Operations Research Letters 10, 27–36.
Desrosiers,J.,F.Soumis,andM.Desrochers(1984).Routingtimewindowsbycolumngeneration.
Networks 14, 545–565.
Dessouky, M., R. Hall, A. Nowroozi, and K. Mourikas (1999). Bus dispatching at timed transfer
transit stations using bus tracking technology. Transportation Research Part C 7, 187–208.
Dror, M. (1993). Modeling vehicle routing with uncertain demands as a stochastic program:
propertiesof the correspondingsolution. European Journal of Operational Research 64, 432–441.
Dror,M.,G.Laporte,andP.Trudeau(1989).Vehicleroutingwithstochasticdemands: properties
and solution framework. Transportation Science 23(3).
Dumas, Y., J. Desrosiers, and S. F. (1991). The pickup and deliveryproblem with time windows.
European Journal of Operational Research 54, 7–22.
114
Eilon, S., D. D. T. Watson-Gandy, and N. Christofides (1971). Distribution Management: Math-
ematical Modeling and Practical Analysis. London: Griffin.
El-Ghaoui, L. and H. Lebret (1997). Robust solutions to least-squareproblems to uncertain data
matrices. SIAM J. Matrix Anal. Appl. 18, 1035–1064.
El-Ghaoui, L., M. Oks, and F. Oustry (2003). Worst-case value at risk and robust portfolio
optimization: A conic programming approach. Operations Research 51(4), 543–556.
El-Ghaoui, L., F. Oustry, and H. Lebret (1998). Robust solutions to uncertain semidefinite
programs. SIAM J. Optimization 9(1).
Erera, A. (2000). Design of Large-scale Logistics Systems for Uncertain Environments. Ph. D.
thesis, University of California, Berkeley.
Ermoliev, Y. and R.-B. Wets (1988). Numerical Techniques for Stochastic Optimization. Berlin:
Springer-Verlag.
Federgruen, A. and D. Simchi-Levi (1995). Analytical analysis of vehicle routing and inventory-
routing problems. In M. Ball, T. Magnanti, C. Monma, and G. Nemhauser (Eds.), Network
Routing, Volume 8 of Handbooks in Operations Research and Management Science, pp. 297–374.
Amsterdam: Elsevier Science.
Fisher, M. (1995). Vehicle routing. In M. Ball, T. Magnanti, C. Monma, and G. Nemhauser
(Eds.), Network Routing, Volume 8 of Handbooks in Operations Research and Management Sci-
ence, pp. 1–33. Amsterdam: Elsevier Science.
Fisher, M. L. and R. Jaikumar (1981). A generalized assignment heuristic for vehicle routing.
Networks 11, 109–124.
Foster, B. and D. Ryan (1976). An integer programming apprach to the vehicle scheduling
problem. Operational Research Quarterly 27, 367–384.
Fukasawa, R., H. Longo, J. Lysgaard, M. P. de Arag˜ ao, M. Reis, E. Uchoa, and F. R. Werneck
(2006). Robust branch-and-cut-and-price for the capacitated vehicle routing problem. Mathe-
matical Programming 106, 491–51.
Gendreau, M., A. Hertz, and G. Laporte (1991). A tabu search heuristic for the vehicle routing
problem. Publication 777, Centre de Recherche sur les Transports, Montreal.
Gendreau, M., G. Laporte, and R. Seguin (1995). An exact algorithm for the vehicle routing
problem with stochastic customers and demands. Transportation Science 29(2), 143–155.
Gendreau, M., G. Laporte, and R. Seguin (1996). Stochastic vehicle routing. European Journal
of Operational Research 88(1), 3–12.
Gent, I. P. and T. Walsh (1996). The TSP phase transition. Artificial Intelligence 88, 349–358.
Gillett, G. and L. Miller (1974). A heuristic algorithm for the vehicle dispatch problem. Opera-
tions Research 22, 340–349.
Golden, B. and A. Assad (Eds.) (1988). Vehicle Routing: Methods and Studies. Amsterdam:
North Holland Publication.
115
Goldfarb, D. and G. Iyengar (2003). Robust portfolio selection problems. Mathematics of Oper-
ations Research 28(1), 1–38.
Guslitser, E. (2002). Uncertainty-immunized solutions in linear programming. Master’s thesis,
Minerva Optimization Center, Technion. iew3.technion.ac.il/Labs/Opt/index.php?4.
Hall, R. W. (1996). Pickup and delivery systems for overnight carriers. Transportation Research
A (30), 173–187.
Hall, R. W. and J. G. Partyka (1997). Vehicle routing, on the road to efficiency. OR/MS Today.
June.
Haughton,M.(1998). Theperformanceofroutemodificationanddemandstabilizationstrategies
in stochastic vehicle routing. Transportation Research 32B(8), 551–566.
Haughton,M.(2000). Quantifyingthebenefits ofroutereoptimizationunderstochasticcustomer
demands. Journal of Operational Research Society (51), 320–322.
Haughton, M. and A. Stenger (1998). Modeling the customer service performance of fixed-rotes
delivery systems under stochastic demands. Journal of Business Logistics (9), 155–172.
Homberger, J. and H. Gehring (1999). Two evolutionary metaheuristics for the vehicle routing
problem with time windows. INFOR 37(3), 297–318.
Jaillet, P. (1987). Stochastic routing problem. In G. Andreatta, F. Mason, and P. Serafini
(Eds.), Stochastics in Combinatorial Optimization, pp. 178–186. New Jersey: World Scientific
Publishing.
Jaillet, P. (1988). A priori solution of a traveling salesman problem in which a random subset of
the customers are visited. Operations Research 36(6), 929–936.
Jaillet, P. and A. Odoni (1988). The probabilistic vehicle routing problem. In B. L. Golden and
A.A.Assad(Eds.),Vehicle Routing: Methods and Studies.North-Holland,Amsterdam: Elsevier.
J´ ez´ equel,A. (1985). Probabilistic vehicle routing problems. Master’s thesis, Department of Civil
Engineering, Massachusetts Institute of Technology.
Jones, K. L., I. J. Lustig, J. M. Farwolden, and W. B. Powell (1993). Multicommodity network
flows: The impact of formulation on decomposition. Mathematical Programming 62, 95–117.
Jula, H., M. M. Dessouky, and P. Ioannou (2006). Truck route planning in non-stationary
stochastic networks with time-windows at customer locations. IEEE Transactions on Intelligent
Transportation Ssystems 37, 51–63.
Kamoun, M. and R. W. Hall (1996). Design of express mail services for metropolitan regions.
Journal of Business Logistics 17(2), 265–320.
Koskosidis, Y., W. Powell, and M. Solomon (1992). An optimization-based heuristic for vehicle
routing and scheduling with soft time window constraints. Transportation Science 26(2), 69–85.
Kouvelis, P. and G. Yu (1997). Robust discrete optimization and its applications. Norwell, MA:
Kluwer Academic Publishers.
116
Lambert, V., G. Laporte, and F. Louveaux (1993). Designing collection routes through bank
branches. Computers & Operations Research 20(7), 783–791.
Laporte, G. (1992). The vehicle routing problem: An overview of exact and approximate algo-
rithms. European Journal of Operational Research 59, 345–358.
Laporte, G., M. Gendreau, J. Potvin, and F. Semet (2000). Classical and modern heuristics for
the vehicle routing problem. International Transactions in Operations Research 7, 285–300.
Laporte, G., F. Louveaux, and H. Mercure (1992). The vehicle routing problem with stochastic
travel times. Transportation Science 26(3), 161–170.
Laporte, G., F. V. Louveaux, and L. V. Hamme (2002). An integer l-shaped algorithm for
the capacitated vehicle routing problem with stochastic demands. Operations Research 50(3),
415–423.
Laporte, G., H. Mercure, and Y. Nobert (1992). A branch and bound algorithm for a class
of asymmetrical vehicle routeing problems. Journal of the Operational Research Society 43(5),
469–481.
Laporte, G., Y. Norbert, and D. M. (1985). Optimal routing under capacity and distance re-
strictions. Operations Research 33, 1050–1073.
Laporte, G. and I. Osman (1995). Routing problems: a bibliography. Annals of Operations
Research 61, 227–262.
Lenstra, J. K. and A. H. G. Rinnooy Kan (1981). Complexity of vehicle routing and scheduling
problems. Networks 11, 221–227.
Lu, Q. and M. Dessouky (2004). An exact algorithm for the multiple vehicle pickup and delivery
problem. Transportation Science 38(4), 503–514.
Lu, Q. and M. M. Dessouky (2006). A new insertion-based construction heuristic for solving
the pickup and delivery problem with hard time windows. European Journal of Operational
Research 175, 672–687.
Lutgens, F. and J. Sturm (2002, November). Robust option modelling. Technical report, Uni-
versity of Maastricht.
Lysgaard, J., N. A. Letchford, and W. R. Eglese (2004). A new branch-and-cut algorithm for
the capacitated vehicle routing problem. Mathematical Programming 100, 423–445.
Meyer(1996). Mohaddesassociates,inc,gatewaycitiestruckingstudy. Technicalreport,Gateway
Cities Council of Governments Southeast Los Angeles County.
Morales,J.C.(2006). Planning Robust Freight Transportation Operations. Ph.D.thesis, Gerogia
Institute of Technology.
Mulvey, J. M., R. J. Vanderbei, and S. A. Zenios (1995). Robust optimization of large-scale
systems. Operations Research 43, 264–281.
O˜ guz, O. (2000). Bounds on the opportunity cost of neglecting reoptimization in mathematical
programming. Management Science 46(7), 1009–1012.
117
Ord´ o˜ nez, F. and J. Zhao (2004). Robust capacity expansionofnetwork flows. USC-ISE Working
paper 2004-01. To appear in Networks.
Orloff (1976). Route-constrained fleet scheduling. Transportation Science 10, 149–168.
Orponen, P. and H. Mannila (1987). On approximation preserving reductions: Complete prob-
lems and robust measures. Technical Report C-1987-28, Department of Computer Science, Uni-
versity of Helsinki.
Osman,I.H.(1993). Metastrategysimulatedannealingandtabusearchalgorithmsforthevehicle
routing problem. Annals of Operations Research 41, 421–451.
Palmer,K., M. M. Dessouky, and T. Abdelmaguid (2004). Impacts ofmanagementpracticesand
advanced technologies on demand responsive transit systems. Transportation Research, Part A:
Policy and Practice 38, 495–509.
Pinar, M. (2002). Linear huber m-estimator under ellipsoidal data uncertainty. BIT Numerical
Mathematics 42(4), 856–866.
Potvin, J.-Y. and S. Bengio. (1996). The vehicle routing problem with time windows. part ii:
Genetic search. INFORMS Journal on Computing 8(2), 165–172.
Powell, W., P. Jaillet, and A. Odoni (1995). Stochastic and dynamic networks and routing. In
M. Ball, T. Magnanti, C. Monma, and G. Nemhauser (Eds.), Network Routing, Volume 8 of
Handbooks in Operations Research and Management Science. Amsterdam: Elsevier Science.
Psaraftis, H. N. (1995). Dynamic vehicle routing: status and prospects. Annals of Operations
Research 61, 143–164.
Quadrifoglio, L. and M. M. Dessouky (2007). An insertion heuristic for scheduling mobility
allowance shuttle transit (mast) services. Journal of Scheduling 10, 25–40.
Ralphs, T. (2003). Parallel branch and cut for vehicle routing. Parallel Computing 29(5), 607–
629.
Ralphs, T., L. Kopman, W. Pulleyblank, and L. Trotter (2003). On the capacitated vehicle
routing problem. Mathematical Programming 94(2-3), 343–359.
Rao,M.R.andS.Zionts(1968). Allocationoftransportationunitstoalternativetrips-acolumn
generation scheme with out-of-kilter subproblems. Operations Research 16, 52–63.
Savelsbergh,M.W. P.andM.Goetschalckx(1995). Acomparisonoftheefficiencyoffixedversus
variable vehicle routes. Journal of Business Logistics (19), 163–187.
Savelsbergh, M. W. P. and N. Sol (1998). Drive: Dynamic routing of independent vehicles.
Operations Research 46(4), 474–490.
Schultz, R., L. Sougie, and M. van der Vlerk (1998). Solving stochastic programs with complete
integer recourse by enumeration: A framework using Grobner basis reductions. Mathematical
Programming 83, 229–252.
Secomandi, N. (2001). A rollout policy for the vehicle routing problem with stochastic demands.
Operations Research 49(5), 796–802.
118
Sexton, T. and Y. Choi (1986). Pickup and delivery of partial loads with soft time windows.
American Journal of Mathematical and Management Sciences 6, 369–398.
Solomon, M. M. (1987). Algorithms for the vehicle routing and scheduling problems with time
window constraints. Operations Research 35, 254–265.
Soyster, A. L. (1973). Convex programming with set-inclusive constraints and applications to
inexact linear programming. Operations Research 21(5), 1154–1157.
Stein, D. (1978). An asymptotic probabilistic analysis of a routing problem. Mathematics of
Operations Research 3, 89–101.
Swihart, M. R. and J. D. Papastavrou (1999). A stochastic and dynamic model for the single-
vehicle pick-up and delivery problem. European Journal of Operational Research 114, 447–464.
Taillard, E. (1993). Parallel iterative search methods for vehicle routing problems. Networks 23,
661–673.
Taillard, E., P. Badeau, M. Gendreau, F. Guertin, and J. Potvin (1997). A tabu search heuristic
for the vehicle routing problem with soft time windows. Transportation Science 31(2), 170–186.
Takriti, S. and J. Birge (2000). Lagrangian solution techniques and bounds for loosely coupled
mixed-integer stochastic programs. Operations Research 48(1), 91–98.
Toth, P. and D. Vigo (2002a). Models, relaxations and exact approaches for the capacitated
vehicle routing problem. Discrete Applied Mathematics 123, 487–512.
Toth, P. and D. Vigo (Eds.) (2002b). The Vehicle Routing Problem. Monographs on Discrete
Mathematics and Applications. SIAM.
Wassan, N. and I. Osman (2002). Tabu search variants for the mix fleet vehicle routing problem.
Journal of the Operational Research Society 53(7).
Waters, C. D. J. (1989). Vehicle-scheduling problems with uncertainty and omitted customers.
The Journal of the Operational Society 40(12), 1099–1108.
Willard, J. A. G. (1989). Vehicle routing using r-optimal tabu search. Master’s thesis, Manage-
ment School, Imperial College, London.
Wren, A. (1971). Computers in Transport Planning and Operation. London: Ian Allan.
Wren, A. and A. Holliday (1972). Computer scheduling of vehicles from one or more depots to
a number of delivery points. Operations Research Quarterly 23, 333–344.
Yaman, H., O. Kara¸ san, and M. Pinar (2001). The robust spanning tree problem with interval
data. Operations Research Letters 29(1), 31–40.
Zhang, W. (2004). Phase transitions and backbones of the asymmetric traveling salesman prob-
lem. Journal of Artificial Intelligence Research 21, 471–497.
Zhang, W. and R. E. Korf(1996). A study ofcomplexity transitions on the asymmetric traveling
salesman problem. Artif. Intell. 81(1-2), 223–239.
119
Zhong, H., R. W. Hall, and M. M. Dessouky (2007). Territory planning and vehicle dispatching.
Transportation Science (41), 74–89.
Zhou, K., J. Doyle, and K. Glover (1996). Robust and optimal control. Upper Saddle River, NJ:
Prentice Hall.
120
Abstract (if available)
Abstract
We consider a Vehicle Routing Problem (VRP) under uncertainty and seek to find efficient solutions in practice. Routing problems contain a fairly high degree of uncertainty in real life
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Vehicle routing and resource allocation for health care under uncertainty
PDF
New approaches for routing courier delivery services
PDF
Routing and inventory model for emergency response to minimize unmet demand
PDF
Routing problems for fuel efficient vehicles
PDF
Routing for ridesharing
PDF
Dynamic programming-based algorithms and heuristics for routing problems
PDF
Models and algorithms for energy efficient wireless sensor networks
PDF
Models and solution approaches for facility location of medical supplies for large-scale emergencies
PDF
Train routing and timetabling algorithms for general networks
PDF
Models and algorithms for pricing and routing in ride-sharing
PDF
Discounted robust stochastic games with applications to homeland security and flow control
PDF
Train scheduling and routing under dynamic headway control
PDF
A stochastic employment problem
PDF
Continuous approximation formulas for cumulative routing optimization problems
PDF
An online cost allocation model for horizontal supply chains
PDF
Robust adaptive control for unmanned aerial vehicles
PDF
Computational geometric partitioning for vehicle routing
PDF
Solving the empty container problem using double-container trucks to reduce vehicle miles
PDF
Continuous approximation formulas for location and hybrid routing/location problems
PDF
Strategies for effective rail track capacity use
Asset Metadata
Creator
Sungur, Ilgaz
(author)
Core Title
The robust vehicle routing problem
School
Viterbi School of Engineering
Degree
Doctor of Philosophy
Degree Program
Industrial and Systems Engineering
Publication Date
07/23/2007
Defense Date
06/14/2007
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
insertion heuristic,large scale stochastic optimization,OAI-PMH Harvest,robust optimization,vehicle routing
Language
English
Advisor
Ordonez, Fernando (
committee chair
), Dessouky, Maged M. (
committee member
), Koenig, Sven (
committee member
)
Creator Email
sungur@usc.edu
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-m642
Unique identifier
UC1147163
Identifier
etd-Sungur-20070723 (filename),usctheses-m40 (legacy collection record id),usctheses-c127-527694 (legacy record id),usctheses-m642 (legacy record id)
Legacy Identifier
etd-Sungur-20070723.pdf
Dmrecord
527694
Document Type
Dissertation
Rights
Sungur, Ilgaz
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Repository Name
Libraries, University of Southern California
Repository Location
Los Angeles, California
Repository Email
cisadmin@lib.usc.edu
Tags
insertion heuristic
large scale stochastic optimization
robust optimization
vehicle routing