University of Southern California Dissertations and Theses

Applications of explicit enumeration schemes in combinatorial optimization

(USC Thesis Other)

Applications of explicit enumeration schemes in combinatorial optimization

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content Applications of Explicit Enumeration Schemes in Combinatorial Optimization
by
Julien Yu "Yue"
A Dissertation Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(INDUSTRIAL AND SYSTEMS ENGINEERING)
May 2024
Copyright 2024 Julien Yu "Yue"

Dedication
I dedicate this dissertation to the enduring enigma and profound wisdom of Beardwood, Halton, and Hammersley, whose work continues to unravel mysteries and inspire curiosity after decades.
ii

Acknowledgements
I wish to extend my deepest gratitude to those who have been instrumental in my journey through the
PhD program at the University of Southern California. Their support, guidance, and inspiration have been
the bedrock of my academic and personal growth.
Foremost, I express my profound appreciation to my faculty mentor, Professor John Gunnar Carlsson.
His exceptional research acumen, admirable personality, and fervent zest for life have been a constant
source of motivation. Working under Professor Carlsson’s guidance has been an extraordinary privilege;
his wisdom and mentorship have been pivotal in shaping my academic journey.
I am immensely thankful to the esteemed members of my committee: Professor Meisam Razaviyayn,
Professor Ketan Savla, Professor Renyuan Xu, and Professor Maged Dessouky. Their insightful feedback
and invaluable expertise have been crucial in steering my dissertation towards excellence. Their contributions have significantly enriched my research and academic experience.
I acknowledge the Center for Advanced Research Computing (CARC) at the University of Southern
California for providing computing resources that have contributed to the research results reported within
this dissertation.
A special note of thanks goes to my core collaborators, Dr. Bo Jones and Ke Xu, whose brilliance and
collaborative spirit have greatly enriched my research. Additionally, I am grateful to my peer mentor, Dr.
Ziyu He, and all my labmates – Dr. Xiangfei Meng, Dr. Jiachuan Chen, Dr. Mohammadjavad Azizi, Dr.
iii

Haochen Jia, Dr. Shannon Sweitzer, Dr. Ying Peng, and Dr. Han Yu – for their camaraderie and shared
insights.
On a personal note, my heartfelt thanks to my family – my mom, dad, grandma, cousins, extended
family members, and my dearly missed grandparents. Their unwavering support and belief in me have
been the foundation of my strength and perseverance throughout these five years. Their encouragement
have been a guiding light in my life.
Lastly, I reserve a special mention for my girlfriend, Yixin Yi. Her unconditional care, love, and companionship have been my greatest comfort and joy. Meeting Yixin has been the most beautiful chapter of
my life, providing a source of happiness and balance amidst the rigors of academic pursuit.
iv

Table of Contents
Dedication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . ii
Acknowledgements . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . iii
List of Tables . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii
List of Figures . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . viii
Abstract . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . xi
Chapter 1: Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 1
Chapter 2: Literature Review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5
2.1 Foundation of the Beardwood-Halton-Hammersley theorem . . . . . . . . . . . . . . . . . 5
2.2 Exact algorithms and relaxation for the travelling salesman problem . . . . . . . . . . . . . 6
2.3 Approximations for the travelling salesman and vehicle routing problems . . . . . . . . . . 8
2.4 Heuristic approaches . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 9
2.4.1 Heuristic algorithms for the travelling salesman problem . . . . . . . . . . . . . . . 9
2.4.2 Heuristic algorithms for the capacitated vehicle routing problem . . . . . . . . . . 10
2.4.3 Heuristic algorithms for the split-delivery vehicle routing problem and its variant
with time windows . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.5 Estimates of the Euclidean travelling salesman constant . . . . . . . . . . . . . . . . . . . . 13
2.6 More proportionality constants and their upper bounds . . . . . . . . . . . . . . . . . . . . 15
2.7 Notational conventions . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17
Chapter 3: Preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
3.1 Exact bounds of the Euclidean travelling salesman constant β2 . . . . . . . . . . . . . . . . 18
3.1.1 Initial upper bound by Beardwood, Halton, and Hammersley . . . . . . . . . . . . 18
3.1.2 Improvement of the exact upper bound . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.1.3 Improvement of the exact lower bound . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.2 A priori splitting rule for the split-delivery vehicle routing problem . . . . . . . . . . . . . 24
Chapter 4: Minimum Size Coalescing Partitions . . . . . . . . . . . . . . . . . . . . . . . . . . 28
Chapter 5: Advancements in TSP through Enumeration . . . . . . . . . . . . . . . . . . . . . 38
5.1 Enumeration of path permutations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
5.2 Construction of high-dimensional integral . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
5.3 Separation of the integral . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
v

5.4 Algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 53
Chapter 6: Computational Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
6.1 A new upper bound of the travelling salesman constant . . . . . . . . . . . . . . . . . . . . 58
6.1.1 Selection of hyperparameters . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
6.1.2 Construction of hyperrectangles with a decision tree . . . . . . . . . . . . . . . . . 61
6.1.3 Computational results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69
6.2 A lossless splitting rule for split-delivery routing problems . . . . . . . . . . . . . . . . . . 72
6.2.1 Limiting the non-split problem size . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
6.2.2 Adding time windows . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 80
Chapter 7: Conclusion and Future Work . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
7.1 Lower bound of the travelling salesman constant . . . . . . . . . . . . . . . . . . . . . . . . 85
Bibliography . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
Chapter 8: Appendix . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
vi

List of Tables
4.1 The size of the minimum coalescing partition µ across various choices of n and k.
Fundamentally, µ embodies the partition with the least number of elements that can
coalesce to all partitions λ ⊢ n where |λ| ≤ k. . . . . . . . . . . . . . . . . . . . . . . . . . 35
6.1 For zigzag configurations involving 5-tuples and 6-tuples, a prioritized task is to
determine the improvement limit for the upper bound of β2 across different settings of
hyperparameter h. The right-hand side of (5.1) is simulated for each h and k pair, with the
95% confidence interval obtained through bootstrapping. . . . . . . . . . . . . . . . . . . . 59
6.2 The best upper bounds of β2 obtained by our method, for increasing values of N. The
overall best upper bound β2 < 0.9038 is highlighted. . . . . . . . . . . . . . . . . . . . . . 71
6.3 LKH-3 results for SDVRP with rounded costs and limited fleet on the instances of
Belenguer et al. [20] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
6.4 LKH-3 results for SDVRP with limited fleet on the instances of Archetti et al. [13] . . . . . 82
6.5 LKH-3 results for SDVRPTW with limited fleet on the Group C and RC instances of
Solomon [90] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 83
vii

List of Figures
3.1 Partitioning of the unit square into horizontal strips. In 3.1a, the unit square is partitioned
into a number of strips, traversed horizontally from left to right. The additional length
needed to link these paths into one continuous tour is shown in 3.1b using dashed lines.
It can be routinely demonstrated that their contribution to the total travelling salesman
tour length is almost surely o(
√
n). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
5.1 A depiction of a 5-tuple (k = 4) where the permutation (0,2,3,1,4) among the six possible
permutations achieves the lowest-cost path. . . . . . . . . . . . . . . . . . . . . . . . . . . 39
5.2 A representation of all six permutations for a 5-tuple (k = 4), highlighting that while the
identity permutation often results in the lowest traversal cost, alternative permutations
may also offer shorter paths. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 39
5.3 The number of traversal options within each horizontal strip increases superexponentially
with k. Additionally, a larger k allows for greater flexibility in selecting the optimal route
for traversing the horizontal strip, as illustrated in 5.3b and 5.3c. . . . . . . . . . . . . . . . 40
5.4 5.4a shows simulations for 5-tuples, with the x-axis representing the ratio of the shortest
path out of 6 options to the identity permutation (0,1,2,3,4). 5.4b illustrates simulations
for 6-tuples, where the x-axis indicates the ratio of the shortest path out of 24 options
to the identity permutation (0,1,2,3,4,5). The juxtaposition suggests that an increase in
k diminishes the likelihood of the identity permutation being the most cost-effective for
achieving the shortest path. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41
5.5 The piecewise quadratic bounding method is utilized to determine reliable lower and
upper bounds for the final summand in (5.4), applicable to any chosen valid hyperrectangle
B. The trapezoidal rule is employed for determining the upper bounds of G(xk) and
e
−xk , and the midpoint rule is adopted for calculating their lower bounds. Importantly,
during computational experiments, the lower bounding technique is exclusively applied
to the identity permutation, while the upper bounding technique is reserved for assessing
all alternative permutations. The capability to secure positive contributions under such
stringent criteria unequivocally substantiates the improvement of β2’s upper bound. . . . 51
viii

6.1 The computational experiments for h
2 values of 3, 3.5, 4, 4.5, 5, and 5.5, across a range
of k from 3 to 9, showcase the potential improvement limits for the upper bound of β2.
This investigation reveals that the optimal choice of h tends to increase with larger k,
suggesting the inadequacy of the default h =
√
3 when computational capabilities permit
the exploration of higher k values. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
6.2 Illustration of a decision tree in an experiment with a small sample of N = 1000
hyperpoints. Each node within the tree represents a subset of rows in the original
dataset, with non-leaf nodes carrying decisions based on specific features and threshold
values. Leaf nodes, which can emerge at any depth beyond the root, are defined by a
homogeneity in their labels (permutation indices). Additionally, every node is associated
with an ancestor attribute, a list that traces its lineage back to the root, with its length
corresponding to the node’s depth in the tree. The root node’s list of ancestors is notably
empty, highlighting its origin status in the tree structure. . . . . . . . . . . . . . . . . . . . 63
6.3 Given k = 4, an example hyperrectangle B ⊂ D ⊂ R
9
is depicted as a combination
of one line segment and four rectangular regions. Additionally, the representation of
a 9-dimensional hyperpoint as five different points on the 2-dimensional coordinate
plane bridges the gap between abstract high-dimensional data and its tangible graphical
interpretation. For this specific hyperpoint, the visual contrast between the identity
permutation (0,1,2,3,4) and the optimal alternative permutation (0,1,3,2,4) is highlighted.
For a specified hyperrectangle B, the goal is to discern an alternative permutation that,
on average, outperforms the identity permutation across the infinite set of hyperpoints
within this hyperrectangle. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
6.4 Comparative analysis of the upper bound of β2 against sample size N, for hyperparameter
values h
2 = 3, 3.25, 3.5. The solid lines depict the new guaranteed upper bounds of
β2, which decrease with increasing N and significantly outperform Steinerberger’s
benchmark. The bars indicate computation time, revealing an expected increase with
larger sample sizes. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
6.5 The a priori splitting rule, when applied to Instance Set 1 from [20], results in varying sizes
of non-split capacitated vehicle routing problem instances for different (¯k, q) combinations.
The generated problem size increases with larger ¯k or smaller q. It is important to note
that when q ≥ dmax/Q, the resulting problems effectively disallow splitting, hence such
instances are excluded from Figure 6.5. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 76
6.6 The gap between the solution obtained using the a priori splitting rule and the best known
solution, for varying selections of (
¯k, q), analyzed across all instances of Set 1 [20]. . . . . 77
7.1 Figure 7.1a demonstrates that a unit square containing n points, each independently and
uniformly sampled, can be divided into n/k square cells. These cells are designed to
permit cost-free traversal along their boundaries, and each has a side length calculated
as p
k/n. The quantity of points per cell adheres to a Poisson distribution with mean k.
Likewise, the unit square can be divided into n/k hexagonal regions, facilitating more
cost-free movements along their boundaries, as depicted in Figure 7.1b. In these hexagons,
the distribution of points remains Poisson with mean k. The side length for each hexagon
is determined by 2
0.5
· 3
−0.75p
k/n. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
ix

7.2 The simulated lower bound of β2 was explored for both square and hexagonal configurations. For each value of k, spanning from 5 to 39, a total of 10, 000 instances were
simulated. The sample averages and standard deviations derived from these simulations
were then utilized to approximate the ℓi values for each scenario. The yellow region
highlights the 95% confidence interval for the average lower bound of β2, revealing that
achieving a value of k = 29 is necessary to exceed the best known lower bound reported
in [46]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
x

Abstract
Enumeration schemes play a pivotal role in addressing the challenge of complex optimization problems.
This dissertation embarks on a detailed exploration of two vital topics in combinatorial optimization, refining the upper bound of the Euclidean Travelling Salesman Problem (TSP), and providing novel insights
into efficient resource allocation in the Split-Delivery Vehicle Routing Problem (SDVRP) and the SplitDelivery Vehicle Routing Problem with Time Windows (SDVRPTW). The integration of these two themes
offers a novel perspective in the field of combinatorial optimization, combining theoretical innovation with
practical problem-solving approaches.
At the outset, this dissertation focuses on a fundamental challenge in combinatorial optimization:
improving the upper bound of the Euclidean TSP constant. Let X1, X2, . . . , Xn be n independent and
uniformly distributed random points in a compact region R ⊂ R
2 of area 1. Let TSP(X1, . . . , Xn) denote the length of the optimal Euclidean travelling salesman tour that traverses all these points. The
classical Beardwood-Halton-Hammersley theorem proves the existence of a universal constant β2 such
that TSP(X1, . . . , Xn)/
√
n → β2 almost surely, which satisfies 0.625 ≤ β2 ≤ 0.92116027. This research presents a computer-aided proof using numerical quadrature and tree-based machine learning that
β2 ≤ 0.9038, marking a significant stride in optimizing the Euclidean TSP and underscoring the critical
role of enumeration schemes in unraveling combinatorial optimization problems. This advancement also
opens the door for further improvements with the evolution of computational hardware over time. In
addition to the explicit upper bound improvement, this research also suggests a novel enumeration-based
xi

direction for enhancing the lower bound of the Euclidean TSP constant, although it appears computationally intractable at present.
Simultaneously, the dissertation ventures into the intricate territory of splittable demand in resource
allocation problems, a challenging aspect of combinatorial optimization. It presents an innovative a priori
splitting rule, minimizing the number of pieces while maintaining the feasibility of all possible demand
splitting patterns. This strategy effectively transforms splittable problems into unsplittable ones, leveraging the advanced heuristics developed for the latter. This part of the research demonstrates the versatility
and efficacy of enumeration schemes in devising practical solutions, as validated by computational experiments on large-scale instances, including those involving vehicle routing problems with time-windows.
On the whole, this dissertation underscores the indispensable role of enumeration schemes in tackling
combinatorial optimization problems, paving the way for innovative solutions.
xii

Chapter 1
Introduction
The Beardwood-Halton-Hammersley (BHH) theorem stands as a cornerstone in the probabilistic analysis of combinatorial optimization, offering a fundamental characterization of the optimal length for the
Euclidean travelling salesman tour [19]. As one of the pioneering limit theorems in combinatorial optimization, the BHH theorem is comprehensively articulated as follows:
Theorem 1 (Beardwood-Halton-Hammersley theorem). Let {Xi} be a sequence of independent random
variables uniformly distributed in the unit square, and let TSP({X1, . . . , Xn}) denote the length of the
optimal Euclidean travelling salesman tour of X1, . . . , Xn. There exists a universal constant β2 satisfying
0.625 ≤ β2 ≤ 0.92116027, such that
limn→∞
TSP({X1, . . . , Xn})
√
n
→ β2
almost surely. Moreover, when the elements of {Xi} are independently sampled from a distribution µ with
compact support R ⊂ R
2
, we have
limn→∞
TSP({X1, . . . , Xn})
√
n
→ β2
Z
R
p
f(x) dx ,
1

almost surely, where f(x) is the density of the absolutely continuous part of the distribution µ.
More generally, for every integer d ≥ 2, when the elements of {Xi} are independently sampled from a
distribution µd with compact support R ⊂ R
d
, there exists a universal dimension-dependent constant βd such
that
limn→∞
TSP({X1, . . . , Xn})
n
d−1
d
→ βd
Z
R
f(x)
d−1
d dx ,
almost surely, where f(x) is the density of the absolutely continuous part of the distribution µd.
The value of βd is presently unknown, although previous studies have provided numerical Monte Carlo
estimates for β2 [41, 63, 67, 75, 77, 93, 96], as well as for β3 and β4 [63, 77]. These estimations primarily rely
on the Held-Karp linear programming relaxation [56, 57]. The most contemporary estimates obtained in
[8, 34] are derived from simulating exceedingly large instances, and show with a high degree of confidence
that β2 ≈ 0.71. However, it wasn’t until over five decades after the proposal of the classical BHH theorem
that the initial bounds of 0.625 ≤ β2 ≤ 0.92116027 from [19] were refined, as demonstrated in the studies
[46, 94].
This dissertation presents a meticulous computer-aided proof of the theorem stated below:
Theorem 2. β2 < 0.9038.
Employing a sophisticated combination of computer algebra tools [44, 62, 83], coupled with numerical
quadrature and tree-based machine learning techniques, this segment of the dissertation delves into the
critical role played by enumeration in determining a better guaranteed upper bound for the Euclidean
travelling salesman constant. While the enhancement achieved thus far may still be considered modest,
the methodology’s capacity to leverage the advancing power of computer hardware holds great promise
for significant future developments.
2

The significant role of enumeration schemes in advancing the theoretical aspects of combinatorial
optimization is emphatically highlighted, and equally, their importance in the field of resource allocation,
especially regarding demand splitting, cannot be overstated.
Problems with a splittable formulation, exemplified by cases such as the maximum flow problem [66],
typically offer greater computational feasibility compared to their unsplittable counterparts, primarily due
to their nature as linear programming relaxation. Within the context of vehicle routing, the option to split
demand not only has the potential to markedly enhance the solution quality but also introduces greater
complexity to the problem [14]. Such complexity is the major factor why many advanced VRP solvers did
not include the SDVRP functionality [61, 78, 85, 89].
Drawing inspiration from the research in [29], this dissertation introduces a lossless a priori splitting
rule for the SDVRP and its variant with time windows (SDVRPTW), ingeniously leveraging enumeration
schemes to reduce computational overhead while ensuring high-quality solutions. Splitting each customer
demand into the fewest possible components, this method guarantees the inclusion of the optimal solution and the accommodation of all viable demand allocations to the fleet. This demand splitting strategy
is philosophically similar to the principles underlying sparse signal recovery that motivates the ℓ1 approximation of ℓ0 minimization [26, 100], aiming for high fidelity reconstruction with a minimal set of
components.
Beyond vehicle routing, this method is applicable to a wide range of problems where splitting offers
improvements in solution quality but increases computational complexity compared to the non-split formulation of the problem. Such problems include arc routing [73], facility location [16], and bin packing [21,
27]. However, this research mainly focuses on the split-delivery vehicle routing problem and its variant
with time windows, demonstrating how enumeration schemes could enhance the problem-solving efficiency. In essence, this section of the dissertation underscores the versatility and practical relevance of
enumeration schemes in addressing real-world combinatorial optimization challenges.
3

The dissertation unfolds through a structured narrative beginning with a comprehensive literature
review in Chapter 2, setting the stage with a detailed examination of prior works. It advances into Chapter
3, shedding light on the foundational aspects of the travelling salesman and split-delivery vehicle routing
problems that are essential for this dissertation’s contributions. The core contributions are presented in
Chapters 4 and 5, where novel enumeration-based methodologies for demand splitting in SDVRP and
refining the TSP upper bound are introduced, respectively. Computational experiments validating these
approaches against existing benchmarks are detailed in Chapter 6. The conclusion in Chapter 7 wraps
up the discussion by summarizing the key findings and offering insights, effectively encapsulating this
dissertation’s contributions to the field of combinatorial optimization.
4

Chapter 2
Literature Review
2.1 Foundation of the Beardwood-Halton-Hammersley theorem
The study of the shortest path through n points in a Euclidean space precedes the proposal of the BeardwoodHalton-Hammersley theorem [19]. Consider {xi}, a sequence of points in the unit square [0, 1]2
. Let
TSP(x1, . . . , xn)represent the length of the optimal Euclidean travelling salesman tour through the points
x1, . . . , xn. Prior research [40, 95, 98] observed that for sufficiently large n, TSP(x1, . . . , xn) scales as
O(
√
n). Specifically, it is bounded as follows:
TSP(x1, . . . , xn) ≤ C1
√
n + C2,
for some explicit constants C1, C2. Detailed analysis by [98] demonstrated the existence of a path through
all n points in the unit square with a length not exceeding (2.8n)
1/2 + 2. This bound was later refined
by [40] to (2n)
1/2 + 1.75. Furthermore, [40] proved that the shortest path through all n points in a
unit d-dimensional hypercube would scale as O(n
d−1
d ), forming the foundational basis for the classical
Beardwood-Halton-Hammersley theorem [19].
5

2.2 Exact algorithms and relaxation for the travelling salesman
problem
The Held-Karp dynamic programming algorithm, as elucidated in [55], marks a groundbreaking advancement in addressing the travelling salesman problem. This exact algorithm represents a significant improvement over the superexponentiality of the then-existing exact algorithms, achieving the time complexity of
Θ(n
22
n
). It can optimally solve small travelling salesman instances of 40 cities.
In their subsequent endeavors [56, 57], Held and Karp introduced the concept of minimum weight
1-tree, the composite of a spanning tree augmented with extra edges, to develop a linear programming relaxation for the symmetric travelling salesman problem. This methodology significantly improves computational efficiency and has become a well-established approach in combinatorial optimization. Specifically,
the Held-Karp relaxation for the symmetric travelling salesman problem is characterized by a notable modification: replacement of the last integer programming constraint with a linear constraint. The formulation
can be succinctly represented as follows:
min
xe
X
e∈E
wexe s.t.
X
e∈N(v)
xe = 2 ∀v ∈ V,
X
e∈N(S)
xe ≥ 2 ∀∅ ̸= S ⊂ V,
0 ≤ xe ≤ 1 ∀e ∈ E,
where G = (V, E) is a weighted, undirected network. Here, N(v) denotes the set of all neighboring edges
of a given node v, and N(S) represents the set of all edges that have at least one endpoint in the nonempty
node set S.
6

The Held-Karp relaxation is recognized for its effective lower bound construction for the optimal travelling salesman tour in numerous instances, offering a strong theoretical foundation that has greatly influenced subsequent studies. This relaxation is a crucial component of the branch-and-bound algorithm,
effectively pruning the search space. It is well-established that the Held-Karp relaxation accounts for at
least 2
3
of the optimal Euclidean travelling salesman tour length [87], with various conjectures proposing
this figure to be closer to 3
4
(implying an integrality gap of 4
3
). This statement has been substantiated under
certain conditions [24, 71], but has not been proved in general. The practical utility and theoretical importance of the Held-Karp lower bound persist as a driving force in combinatorial optimization research,
especially in the bottom-up asymptotic analysis [63] of the Euclidean travelling salesman constant, a critical aspect of this dissertation’s major contributions.
Moreover, recent advancements in computational capabilities and algorithmic techniques have led to
a noticeable enhancement in the solvability of larger TSP instances. The current fastest exact algorithm
for the travelling salesman problem, inspired by quantum speedup for dynamic programming across the
Boolean hypercube, achieves a reduced time complexity of O(1.728n
) [5]. Despite this improvement, the
algorithm’s practicality is still limited when addressing large-scale problems. In contrast, the Concorde
TSP solver [7, 9], which integrates the branch-and-cut algorithms, is widely acknowledged as the leading
exact solver for the travelling salesman problem. It has demonstrated its capability by solving a problem
with 85,900 cities [8]. However, this feat necessitates an extensive computational effort of over 136 CPU
years, highlighting the ongoing challenge in finding computationally tractable exact solutions for largescale travelling salesman problems.
7

2.3 Approximations for the travelling salesman and vehicle routing
problems
The challenge of computing exact solutions for the travelling salesman problem has catalyzed extensive
research into approximation and heuristic algorithms. This intensified research interest stems from the
need for practical solutions that are computationally attainable. A significant development in approximation algorithms was the Christofides-Serdyukov algorithm in the late 1970s [31, 86], noted for its ability
to guarantee a solution within 1.5 times the optimal TSP value for any metric space, with a computational
complexity of O(n
3
). This algorithm has been foundational in advancing the field.
Recent algorithmic improvements [47, 71, 72] have enriched the understanding of special instances
of the travelling salesman problem, particularly the graphic travelling salesman problem and the graphic
travelling salesman path problem (TSPP). Notably, [72] refined the approximation factor for the graphic
TSP to 13
9
and introduced an ϵ-approximation for the graphic TSPP. Despite these advancements, the
approximation factor for the metric TSP remained at 1.5 for over forty years. It was only in the 2020s
that a minor improvement was made for the metric TSP, introducing a randomized 3
2 − ϵ approximation
algorithm for some ϵ > 10−36 [65]. The enduring nature of the Christofides-Serdyukov approximation
factor indirectly implies the challenge in establishing better bounds for the Euclidean TSP constant, a
primary focus of this dissertation.
Regarding the approximation algorithms of the capacitated vehicle routing problem (CVRP), the foundational approach has also remained largely unchanged since the 1980s [3, 53]. This method involves
computing a travelling salesman tour and optimally partitioning it into segments, each served by a tour
with the capacity of 1. If the approximation factor of the metric TSP is α, then the CVRP’s approximation
factor is α+ 2, whereas for the split-delivery (SDVRP) version, it is α+ 1. Recently, [22] has demonstrated
8

a slightly improved ϵ-approximation for both the capacitated vehicle routing problem and its split-delivery
relaxation, for any given α > 1.
2.4 Heuristic approaches
2.4.1 Heuristic algorithms for the travelling salesman problem
Heuristic methods, in contrast to approximation algorithms, have shown a consistent capability to generate
high-quality solutions for the travelling salesman problem, the vehicle routing problem, and their variants.
These heuristic approaches are recognized for their practical efficiency in producing near-optimal solutions, and are particularly valuable when computational resources are limited. This efficiency becomes
critical in handling large-scale instances of the travelling salesman and vehicle routing problems, where
exact solution derivation is often impractical.
Central to heuristic approaches for the travelling salesman problem are strategies such as simulated
annealing [1], tabu search [41], and the more recently developed evolutionary algorithm [49, 74]. These
methods are especially valuable when there are strict time constraints, insufficient computational capabilities, and the necessity for immediate high-quality solutions.
In the late 1980s and early 1990s, there were noteworthy contributions to this field. [1] provided an
extensive analysis of simulated annealing, focusing on its application to a 100-city travelling salesman
tour. [41] extended the tabu search algorithm, originally designed for various combinatorial optimization
problems [48], to large-scale travelling salesman problem instances. This adaptation highlighted the algorithm’s long-term memory and suitability for parallel computing, offering an asymptotic estimation of
the Euclidean TSP constant, the quantity of interest in this dissertation. Moreover, [70] investigated the
integration of simulated annealing with tabu search in the travelling salesman problem, enhancing the
solution process through parallel computing and achieving significant computational speed-ups.
9

The application of genetic algorithms to the travelling salesman problem was another significant development, as demonstrated by [49]. This research emphasized the role of representation and recombination
operators in developing near-optimal tours. [28] further affirmed that in the context of the travelling
salesman problem, the genetic algorithm’s global convergence property was preserved using a generalized mutation operator. In the 21st century, Nagata’s work [74] introduced a highly effective evolutionary
algorithm for the TSP, employing the edge assembly crossover within an evolutionary framework. This
method proved particularly effective for large-scale instances surpassing 10,000 cities, and set new records
for several previously unsolved national TSP instances, marking a significant leap in TSP heuristic development.
Among the various heuristic algorithms for the travelling salesman problem, the dynamic local search
algorithm developed by Lin and Kernighan [68] is particularly noteworthy for its solution quality and
applicability to the vehicle routing problem and its variants. The Lin-Kernighan algorithm begins with an
existing Hamiltonian cycle and iteratively examines neighboring tours to find shorter routes until a local
minimum is achieved. This algorithm is distinct from fixed-step methods such as 2-opt or 3-opt in that it
does not limit itself to a preset number of edge rearrangements, although it typically involves fewer edge
replacements, and uniquely allows for the reversal of the travel direction. The algorithm’s time complexity
is O(n
⌊
p1
2
⌋
), where p1 represents the backtracking depth.
2.4.2 Heuristic algorithms for the capacitated vehicle routing problem
The Clark-Wright algorithm was introduced [32] during the early 1960s. This algorithm represents a significant advancement in logistics and distribution networks, addressing the complexities of the capacitated
vehicle routing problem (CVRP), which are more intricate compared to the travelling salesman problem.
The Clark-Wright algorithm was pioneering in its introduction of an iterative process for rapidly identifying optimal or near-optimal routes. Its capacity to manage the exponentially growing number of potential
10

routes as delivery points increase was particularly notable, achieved by merging shorter trips into longer
ones.
The versatility and efficiency of the Clark-Wright algorithm have established it as a foundation for
subsequent researches, offering pragmatic solutions to complex routing challenges that are essential in
logistics. Further enhancements [4, 80] have been made to the original Clark-Wright savings algorithm,
strengthening its application to the capacitated vehicle routing problem. Moreover, adaptations of the
Clark-Wright algorithm to diverse VRP variants, such as the vehicle routing problem with backhauls [36],
the stochastic vehicle routing problem [39], and the open vehicle routing problem [79], have also demonstrated competitive results.
Regarding the classic Lin-Kernighan heuristic for the travelling salesman problem [68], Helsgaun’s
subsequent enhancements have significantly bolstered its importance [58]. These improvements led to
the creation of the Lin-Kernighan-Helsgaun (LKH) solver [59, 60], which has become a standard in the
industry. The latest update, LKH-3 [61], facilitates heuristic solutions to a broad spectrum of complex
variants of the travelling salesman and vehicle routing problems, accommodating constraints related to
capacity, time windows, pickup-and-delivery, and distance. However, it is notable that the split-delivery
vehicle routing problem (SDVRP), a topic of particular interest in this dissertation, is not directly solvable
with the current LKH-3 solver.
2.4.3 Heuristic algorithms for the split-delivery vehicle routing problem and its variant
with time windows
The research on the split-delivery vehicle routing problem (SDVRP) has become a significant focus in
logistics, distinguishing itself from its unsplittable counterparts. The key distinction between splittable
and unsplittable problems lies in the ability to divide customer demands across multiple deliveries, leading
to a more complex problem space that demands advanced heuristics and algorithms. Notable developments
11

in this field include remarkable solutions, exemplified by the influential works of [2, 10, 11, 20, 23, 25, 29,
54, 88]. Additionally, for the SDVRP variant that incorporates time window constraints, the split-delivery
vehicle routing problem with time windows (SDVRPTW), effective algorithms have also been proposed,
as illustrated in [12, 38, 69, 84].
The iterated local search heuristic by [88] marked a significant advancement in the research of splitdelivery vehicle routing problem, enhancing the majority of best-known solutions at that time. This multistart iterated local search (ILS) heuristic adapts well to scenarios with both limited and unlimited fleet sizes.
The algorithm’s effectiveness is attributed to its specialized local search operators and its innovative perturbation mechanism, leading to superior performance. Through computational experiments, it was shown
that this algorithm consistently finds better solutions than previous heuristics, with less requirement for
parameter tuning.
The recent study by [54] introduced a memetic algorithm for the split-delivery vehicle routing problem, combining the general edge assembly crossover (gEAX) with local search techniques. This algorithm
excels in creating promising offspring solutions and refining them with extensive local search. Features
like a feasibility-restoring procedure, a diversification-oriented mutation, and a pool updating technique
based on quality and distance all contribute to enhancing its effectiveness. When tested against various
benchmarks for the split-delivery vehicle routing problem, this algorithm outperformed the best-known
upper bounds, primarily set by [88], in numerous instances, showcasing its high performance and practicality.
For the split-delivery vehicle routing problem with time windows, [38] introduced an exact branchand-price-and-cut method, where the column generation subproblem is essentially a resource-constrained
elementary shortest-path problem combined with the linear relaxation of a bounded knapsack problem.
The key contributions of this research include a novel decomposition method for the SDVRPTW and a
specialized label-setting algorithm for solving the subproblem. [12] built upon this with enhancements
12

like a tabu search algorithm for the column-generation subproblem, the extension of valid inequalities,
and a new separation algorithm for k-path inequalities, with computational results demonstrating their
effectiveness.
Last but not the least, the study by [69] addressed the split-collection vehicle routing problem with time
windows and linear weight-related cost, a generalized version of the split-delivery vehicle routing problem with time windows. This research stands out for integrating vehicle weight into cost calculations,
thereby increasing the model’s applicability and generalizability. This research progressed beyond earlier
branch-and-price-and-cut algorithms, effectively addressing complex variants of the split-collection routing problem with time windows, thereby providing competitive solution quality for the SDVRPTW and
surpassing many best-known solutions set by [12, 38].
The evolution of heuristics for the split-delivery vehicle routing problem (SDVRP) and its variant with
time windows (SDVRPTW), as depicted in these studies, demonstrates a movement towards solutions that
are more sophisticated, efficient, and applicable to real-world logistics and supply chain management.
2.5 Estimates of the Euclidean travelling salesman constant
The Beardwood-Halton-Hammersley (BHH) theorem is now recognized as a classic result with significant
implications in the fields of combinatorial optimization and probability theory. Its influence and applicability are well documented in various works [33, 34, 42, 52, 91, 97, 102]. Notably, the theorem’s relevance
extends beyond the unit square to any compact region R ⊂ R
2
, from the uniform distribution to the absolutely continuous part of any distribution µ, from a 2-dimensional square to a hypercube in any dimension
(substituting β2 with a more general but distinct constant βd), and from planar to toroidal topology [63,
77], specifically under periodic boundary conditions that align opposite sides of the unit square.
Numerous studies have provided numerical estimates for the Beardwood-Halton-Hammersley constant in two or higher dimensions [19, 41, 63, 67, 75, 77, 93]. However, these estimates vary based on
13

different methodologies. Initially, Beardwood et al. [19] approximated β2 as 0.749 in their seminal paper,
albeit without any confidence level analysis. In the late 1970s, [93] provided the first simulation-based estimate of β2 at 0.765, with the limitation that computational power at the time could only handle instances
of up to 100 points.
In the late 1980s and early 1990s, [75] analyzed the asymptotic performance of existing TSP heuristics,
using the 3-optimal heuristic to estimate β2 as 0.7425. [67] employed multicanonical annealing on the
travelling salesman problem and approximated β2 to be 0.7211 using simple linear regression analysis.
The most significant estimate until then, by [41], suggested that β2 was upper bounded by 0.731, with a
confidence level of 0.01. This was a secondary result of their study adapting tabu search algorithms to TSP
instances of up to 100,000 cities, which were considered large-scale at the time.
The approximation of the Beardwood-Halton-Hammersley constant βd for dimensions higher than two
(i.e., d > 2) became increasingly attainable with the significant advancements in computing power during
the late 1990s. One of the key studies in this area was conducted by [63], whose work was built upon the
Held-Karp lower bound [56, 57] and undertook a thorough asymptotic experimental analysis of the bounds
of the travelling salesman problem. Their research provided refined estimates of β2 ≈ 0.7124 ± 0.0002,
β3 ≈ 0.6980±0.0003, and β4 ≈ 0.7234±0.0003. These findings effectively invalidated previous estimates
of β2 such as those by [19, 93], and challenged the efficacy of many subsequently proposed heuristics for
the travelling salesman problem.
Concurrently, Percus and Martin [77] conducted research on the finite size and dimensional dependence in the Euclidean travelling salesman problem within a d-dimensional hypercube, incorporating periodic boundary conditions. Their approach, which utilized a mean-field methodology and exploited a notable universality in the k-th nearest neighbor distribution, resulted in estimates of β2 ≈ 0.7120 ± 0.0002
and β3 ≈ 0.6979 ± 0.0002. Other researchers also contributed estimates for the Euclidean travelling
salesman constant, although some were later deemed less accurate than the aforementioned studies. For
14

instance, [96] used iterative methods to estimate β2 ≈ 0.70787±0.0008. These studies collectively marked
a significant advancement in the understanding and approximation of β2 and βd in higher dimensions, contributing substantially to the field of combinatorial optimization and the study of the travelling salesman
constant.
More recently, [34] leveraged their Concorde TSP Solver [7, 9] to provide an estimate of β2 ≈ 0.712403±
0.000007. This latest estimation, supported by extensive simulations of large instances, rendered many
earlier estimates less relevant. Empirical evaluations also suggest a decreasing trend in the estimate of β2
over the number of traversed cities. This trend can be attributed to the advancements in computational capabilities, enabling simulations of larger and more complex instances, and thereby enhancing the accuracy
of these estimates.
2.6 More proportionality constants and their upper bounds
In continuous approximation analysis, particularly within the context of logistics systems, establishing
bounds for proportionality constants is a common and important practice [6, 35, 43]. In fewer problem
settings, the focus is often on obtaining exact expressions for these constants. Notably, the minimum spanning tree problem as discussed in [17] and specific versions of the location-routing problem as explored in
[99] allow for precise formulations of these constants.
Significant progress in this area was marked by [92], which established a limit theorem and provided
proof for a uniform convergence theorem for subadditive Euclidean functionals. Observable in a range
of related problems, such limiting behaviors are applicable to numerous problem settings, including the
Steiner tree, rectilinear Steiner tree, and minimum weight matching problems, as demonstrated in studies
like [81, 91, 92, 102]. This highlights the broad applicability of this analytical approach and inspires ensuing
researches on mathematical constants.
15

Specifically, the research by [45] focuses on a few limiting constants derived from asymptotic formulas,
examining their relationship with the d-dimensional travelling salesman constant denoted as β
d
T SP (or
βd in our notation). Key constants discussed include the minimum spanning tree constant β
d
MST , the
minimum matching constant β
d
MM, the 2-factor constant β
d
T F , and the constant β
d
HK associated with the
Held-Karp relaxation [56, 57]. This research underscores the interconnected nature of these constants and
their relevance to the broader field of combinatorial optimization.
[45] elucidates that, in the context of n random points distributed within a d-dimensional hypercube
[0, 1]d
, a set of crucial strict inequalities consistently hold for any dimension d ≥ 2. The inequalities are
expressed as follows:
β
d
MST < βd
T SP ,
2β
d
MM < βd
T SP ,
β
d
T F < βd
T SP ,
β
d
HK < βd
T SP .
These inequalities highlight the intricate connections between different constants in combinatorial optimization within multi-dimensional spaces, contributing significantly to the comprehension of asymptotic
behaviors.
Furthermore, improving the upper bound of the Euclidean travelling salesman constant β
d
T SP in any dimension has implications beyond the travelling salesman problem itself. Such an improvement inherently
refines the upper bounds of related problems, including the minimum spanning tree, minimum matching,
and 2-factor. Thus, the focus of this dissertation on refining the upper bound of the Euclidean TSP constant is not only pivotal for advancing our understanding of the travelling salesman problem but also has a
16

ripple effect, potentially enhancing the theoretical bounds of these associated combinatorial optimization
problems.
2.7 Notational conventions
In this dissertation, the notation {Xi} is employed with meanings that vary depending on the context.
Initially, within the scope of the Beardwood-Halton-Hammersley Theorem 1, {Xi} denotes a sequence of
independent random variables in a d-dimensional hypercube, where d is an integer, typically chosen to be
2.
Later, in the section addressing the exact upper bound of the 2-dimensional Euclidean travelling salesman constant, the notation {(Xi
, Yi)} is used intermittently to denote a sequence of 2-dimensional independent random variables within the unit square. In this context, Xi and Yi are 1-dimensional random
variables. As explicated in Section 3.1, the sequence {Xi+1 − Xi} conforms to a scaled exponential distribution as n approaches infinity.
Throughout the dissertation, the term TSP({X1, . . . , Xn}) is consistently used to represent the length
of the shortest travelling salesman tour that traverses all points X1, . . . , Xn. For a path that traverses a
fixed number of points, 1-based indexing is utilized in Section 3.1 to establish definitions and elucidate
concepts more clearly. Following Chapter 5, where the contributions of this dissertation are discussed, the
indexing convention transitions to 0-based.
17

Chapter 3
Preliminaries
3.1 Exact bounds of the Euclidean travelling salesman constant β2
Estimations of the 2-dimensional Euclidean travelling salesman constant, denoted as β2, along with its
higher-dimensional counterparts βd, have been a subject of intense research for many decades. Despite
this, precise bounds for these constants have not been much improved since the seminal work by [19]. This
section focuses on discussing the theoretical upper and lower bounds of the Euclidean travelling salesman
constant.
3.1.1 Initial upper bound by Beardwood, Halton, and Hammersley
Beardwood et al. [19] effectively utilized the notion that when n points are independent and uniformly
distributed in a unit square [0, 1]2
, it mirrors a Poisson process with an intensity of n. This connection
between Poisson process [18] and the Euclidean travelling salesman problem was detailed in their initial
work and captured in the following lemma:
Lemma 3. For a Poisson process Pn with intensity n within a unit square [0, 1]2
, the following holds:
limn→∞
E(TSP(Pn))
√
n
→ β2 ,
18

where T SP(Pn) denotes the length of the optimal travelling salesman tour that connects all points in the
Poisson process Pn, and β2 is a constant.
In the Poisson process, the mean and variance of the number of points are both n. This implies that
the deviation of point counts is typically on the order of √
n, which is minor compared to n. The scenario
where the deviation is of higher order of magnitude occurs with an exponentially decreasing probability,
supporting the validity of Lemma 3. The inverse statement that Lemma 3 implies Theorem 1 was also
essentially established in [19].
Utilizing this perspective that relates the travelling salesman problem to Poisson process, [19] further
demonstrated the key principle that for n independent and uniformly sampled points X1, . . . , Xn within
the unit square [0, 1]2
, it holds that
limn→∞
E(TSP({X1, . . . , Xn}))
√
n
→ β2 .
This finding is crucial as it shifts the emphasis to calculating expected tour lengths, rather than depending
on the almost sure context as outlined in Theorem 1, thus facilitating subsequent analyses.
Building on this insight, [19] established their initial upper bound for the Euclidean travelling salesman
constant by constructing a “reasonably good” tour. This was achieved by partitioning the unit square into
horizontal strips and connecting the independent and uniformly distributed points within each strip from
left to right. This process creates individual paths, as depicted in Figure 3.1a. A single tour is then formed
by merging these paths with an additional length that is almost surely o(
√
n), as shown in Figure 3.1b.
The primary goal of this methodology then becomes computing the expected total length of all horizontal
paths within the unit square.
As this dissertation focuses on the discussion and analysis of the Euclidean travelling salesman constant in the 2-dimensional framework, we will adopt a simplified notation for ease of reference. The n
19

{
(a) (b)
Figure 3.1: Partitioning of the unit square into horizontal strips. In 3.1a, the unit square is partitioned into
a number of strips, traversed horizontally from left to right. The additional length needed to link these
paths into one continuous tour is shown in 3.1b using dashed lines. It can be routinely demonstrated that
their contribution to the total travelling salesman tour length is almost surely o(
√
n).
independent and uniformly distributed points within the unit square [0, 1]2 will henceforth be denoted as
(X1, Y1), . . . ,(Xn, Yn). Without loss of generality, the set of points in a given horizontal strip within the
unit square can be denoted as Sh and defined as follows:
Sh =
n
(x, y) : 0 ≤ x ≤ 1, 0 ≤ y ≤
h
√
n
o
∩
n
(X1, Y1), . . . ,(Xn, Yn)
o
.
The set Sh encompasses all points within a horizontal strip of height h/√
n and width 1, considering n
as a variable that could potentially reach infinity. Notably, the parameter h is adjustable, and subsequent
analysis in this dissertation will explore how varying h could lead to favorable computational outcomes
in refining the upper bound of β2.
20

We define a subset of k points in Sh, represented as(X1, Y1), . . . ,(Xk, Yk), to be x-adjacent if and only
if their order statistics along the x-axis are consecutive. Assume n
′
points are in Sh. The conventional
notation for order statistics is employed as follows:
X(1) = min{X1, . . . , Xn
′ }, . . . , X(n
′
) = max{X1, . . . , Xn
′ } .
More precisely, the point set (X1, Y1), . . . ,(Xk, Yk) is defined as x-adjacent if and only if there exists an
integer i such that 0 ≤ i ≤ n
′
− k, and
{X1, X2, . . . , Xk} = {X(i+1), X(i+2), . . . , X(i+k)} .
Now consider a path of k points (X1, Y1), . . . ,(Xk, Yk) that are x-adjacent, where we assume without
loss of generality that X1 ≤ · · · ≤ Xk. It is straightforward that Yi = √
h
n
Ui
, where Ui
is uniform on
the interval [0, 1]. Additionally, the sequence {Xi} conforms to a Poisson process with an intensity of
h
√
n over the interval [0, 1]. As demonstrated in [19], the difference between consecutive x-coordinates,
Xi+1 − Xi
, adheres to a scaled exponential distribution as n → ∞. This can be represented as:
h
√
n(Xi+1 − Xi) ∼ Zi
,
where Zi follows an exponential distribution with a mean of 1. Consequently, the expected distance between each consecutive pair of points (Xi
, Yi) and (Xi+1, Yi+1) can be expressed as:
E



Xi+1 − Xi
Yi+1 − Yi



=
1
h
√
n
E



Zi
h
2
(Ui+1 − Ui)



.
21

Here, Ui+1 and Ui are uniformly distributed over the interval [0, 1]. By multiplying the right hand side of
the above equation by n, an upper bound for the travelling salesman tour length is derived. As a result,
for any given value of h, the following inequality holds:
β2 ≤
1
h
E



Z
h
2
(U0 − U1)



(3.1)
=
1
h
Z ∞
0
Z 1
0
Z 1
0
e
−zp
z
2 + h
4(u0 − u1)
2 du0 du1 dz
=
1
3h
5
Z ∞
0
e
−z

3h
2
z
2
log(h
2
/z +
p
h
4/z2 + 1) + 2z
3 + (h
4 − 2z
2
)
p
h
4 + z
2

dz .
Essentially proposed by [19], the above expression of the upper bound of β2 attains its minimum value at
h =
√
3, yielding an integral value approximately equal to 0.92116027. Despite their significant contributions in this field, Beardwood et al. [19] initially made a miscalculation in their numerical integration
for determining the upper bound of β2 using Simpson’s rule. This resulted in an incorrect estimate of
0.92037. This miscalculation went unnoticed and was repeatedly referenced in subsequent research and
publications, until it was corrected over fifty years later in [94].
3.1.2 Improvement of the exact upper bound
The research in [94] offers a modest refinement to the upper bound of β2. It highlights that the conventional left-to-right ordering described in [19] is not universally optimal, introducing an alternative configuration, termed a “zigzag structure”, that demonstrates improved results in scenarios where an alternative
ordering of points is preferable. This improvement in the upper bound of β2 is quantified by establishing
lower bounds for the probability of occurrences of the zigzag configuration and the resultant improvement. Although the numerical enhancement remains minor, the study underscores the significance of the
conceptual advancement more than the numerical enhancement itself.
22

More specifically, the improvement to the upper bound of β2 proposed by [94] can be expressed as:
β2 ≤ βUB − ϵ0 ,
for some explicit
ϵ0 >
9
16
10−6
.
It is asserted in [94] that this conceptual framework could lead to an estimated improvement of approximately 0.0148 in the upper bound of β2, based on simulations. This claim has been substantiated by
computational experiments conducted in this dissertation. However, it is similar to preceding research
[19, 34, 41, 63, 67, 75, 77, 93, 96] and mainly offers numerical estimations for β2 without practical implications. The primary challenge faced by [94] was in the explicit calculation of expectations, an area where
this dissertation makes a key contribution.
This dissertation introduces a theoretical framework and computer-aided proofs that significantly reduce the upper bound of the 2-dimensional Euclidean travelling salesman constant β2. The results of this
research suggest a notable tightening of the upper bound of β2 to 0.9038. Additionally, the developed
methodology has broader applications, extending to the refinement of the upper bounds for the generalized βd and other problems characterized by subadditive Euclidean functionals.
3.1.3 Improvement of the exact lower bound
In contrast to the complications involved in deriving the upper bound, the original lower bound for β2,
set at 5
8
, was established more directly, employing probabilistic methods and analysis of expected minimums. In their landmark study, Beardwood et al. [19] provided a lower bound for the expected length of
a travelling salesman tour encompassing all points within a unit square, with points sampled according to
a Poisson process of intensity n as specified in Lemma 3.
23

Building upon this, consider the scenario with n independent and uniformly distributed points within
a unit square. In a feasible travelling salesman tour, each point is linked to exactly two other points. By
calculating the expected distance from a point to its nearest neighbor and multiplying this by √
n, a value of
1/2 is obtained. Similarly, calculating the expected distance to the second nearest neighbor and multiplying
by √
n results in 3/4. Since both of the two nearest points contribute equally to the expected length of the
travelling salesman tour, a preliminary lower bound of β2 is deduced to be 5/8.
There has been limited research focused on enhancing the lower bounds of β2 and its higher dimensional counterparts, βd. Notable contributions in this domain have been made by [94] and later by [46].
The study by [94] improved the initial lower bound of β2 by 19/10368, which was then further advanced
by the subsequent research of [46].
The advancements made by [46] revolve around identifying scenarios where short cycles are prevalent
and proposing methods to adjust for these cycles by tightening the triangular inequalities. This approach
leads to a more refined lower bound for β2, which can be represented as:
0.6277 ≤
5
8
+
19
20736
+
3072√
2 − 4325
10752
≤ β2 .
This enhancement to the lower bound of β2 is indeed more substantial compared to the improvements
made to the explicit upper bound. However, [46, 94] acknowledge the need for novel approaches to achieve
further advancements in the lower bounds of β2 and higher dimensional βd.
3.2 A priori splitting rule for the split-delivery vehicle routing problem
This segment of the dissertation presents an overview of the utilization of enumeration schemes in addressing complex forms of combinatorial optimization challenges, specifically focusing on the split-delivery
24

vehicle routing problem (SDVRP) and its time-constrained variation, the SDVRP with time windows (SDVRPTW).
Prior to delving into the concept of a priori splitting scheme, which plays a critical role in the analysis
herein, it is essential to formally define the split-delivery vehicle routing problem. The SDVRP, an intricate
relaxation of the capacitated vehicle routing problem, can be delineated as follows:
Definition 4. [Split-Delivery Vehicle Routing Problem] Consider a directed graph G = (V, E), where
V = {V0, . . . , Vn} denotes a set comprising a depot V0 and n customer nodes. The set of arcs E =
(Vi
, Vj ) : Vi
, Vj ∈ V, i ̸= j signifies all potential routes between these nodes. Each customer node Vi has
a fixed demand Di
, and each arc (Vi
, Vj ) incurs a specific travel cost cij . Given k trucks with capacity Q
such that kQ suffices to meet all demands, the objective is to assign routes that start and end at the depot,
satisfying all customer demands while minimizing the total travel cost. The division of customer demands
among multiple trucks is permitted, allowing various trucks to serve a single customer.
Next, we consider the split-delivery vehicle routing problem with time windows, which incorporates
additional temporal constraints to the SDVRP:
Definition 5. [Split-Delivery Vehicle Routing Problem with Time Windows] Consider a directed graph
G = (V, E), where V = {V0, . . . , Vn} denotes a set comprising a depot V0 and n customer nodes. The set
of arcs E = (Vi
, Vj ) : Vi
, Vj ∈ V, i ̸= j signifies all potential routes between these nodes. Each customer
node Vi
is associated with a demand Di
, a specific time window [ai
, bi], and a service time si
. For each arc
(Vi
, Vj ), there is an associated travel cost cij and travel time tij . With k trucks of capacity Q, the objective
is to assign routes that start and end at the depot, fulfilling all customer demands and minimizing the total
travel cost. Customer demands can be split among trucks, with the additional requirement that all trucks
visiting a customer must arrive within the customer’s time window and stay during the service time.
25

To address the split-delivery vehicle routing problem, the research conducted by [29] proposes two a
priori rules for the segmentation of customer demands. The first rule, referred to as the 20/10/5/1 rule,
partitions the demand Di
in relation to the predetermined truck capacity Q. This rule is elaborated as
follows:
m20 = max{m ∈ N ∪ {0} : 0.20Qm ≤ Di},
m10 = max{m ∈ N ∪ {0} : 0.10Qm ≤ Di − 0.20Qm20},
m5 = max{m ∈ N ∪ {0} : 0.05Qm ≤ Di − 0.20Qm20 − 0.10Qm10},
m1 = max{m ∈ N ∪ {0} : 0.01Qm ≤ Di − 0.20Qm20 − 0.10Qm10 − 0.05Qm5}.
The alternative 25/10/5/1 rule differs primarily in the initial fraction of vehicle capacity designated for
demand segmentation. It is articulated as:
m25 = max{m ∈ N ∪ {0} : 0.25Qm ≤ Di},
m10 = max{m ∈ N ∪ {0} : 0.10Qm ≤ Di − 0.25Qm25},
m5 = max{m ∈ N ∪ {0} : 0.05Qm ≤ Di − 0.25Qm25 − 0.10Qm10},
m1 = max{m ∈ N ∪ {0} : 0.01Qm ≤ Di − 0.25Qm25 − 0.10Qm10 − 0.05Qm5}.
In their computational experiments, [29] established that for numerous benchmark instances of the splitdelivery vehicle routing problem, their splitting heuristic not only achieved rapid convergence but also
yielded solutions that were comparable or superior to those of preceding algorithms [2, 11, 14, 15, 30, 37].
However, it is critical to acknowledge that most of today’s best-known solutions for SDVRP are attributed
to [54, 88].
26

This part of the dissertation delves into the exploration of a priori splitting rules, recognizing their
potential as a robust approach for addressing the SDVRP and its variant with time windows. Their appeal
lies in their straightforwardness relative to mainstream SDVRP heuristics, their adept conversion to the
CVRP framework, and their innovative employment of enumeration techniques. Chapter 4 will introduce a
novel, lossless a priori splitting rule, offering a contrast to the potentially suboptimal rule proposed by [29].
Drawing inspiration from [29, 64] and building upon them, the approach formulated in this dissertation
achieves superior computational outcomes for most benchmark instances compared to [29], meanwhile
circumventing the necessity for advanced metaheuristics [54, 88].
27

Chapter 4
Minimum Size Coalescing Partitions
This chapter delves into the methodology of decomposing a positive integer n into a sum of smaller positive
integers. This foundational concept is essential for developing the lossless a priori splitting rule, which is
utilized to efficiently address the split-delivery vehicle routing problem later on.
For illustration, one way to partition n = 11 is {4, 3, 2, 2}. Such a partition, when comprising k
elements, is referred to as a k-partition. The notation µ ⊢ n is conventionally employed in combinatorics
to indicate that µ is a partition of n, ensuring that P
i µi = n. A key objective of this chapter is to introduce
and prove a method for solving the minimum size coalescing partition problem (MSCP) for any specified
n and k, provided k ≤ n. The discussion herein elaborates on and extends the research introduced in [64],
defining the minimum size coalescing partition problem as follows:
Definition 6 (Minimum Size Coalescing Partition Problem). Given integers n and k, with k < n, a partition µ of n coalesces to a partition λ of n if it is possible to combine elements of µ in a specific manner
to achieve λ. This implies the existence of a decomposition of µ into disjoint subsets µS1
, . . . µS|λ|
, where
for each i ∈ {1, . . . , |λ|}, µSi ⊢ λi
. If µ can coalesce to every λ ⊢ n with |λ| ≤ k, then it is said that µ
coalesces to all k-partitions of n. The objective of the minimum size coalescing partition problem (MSCP)
is to identify the partition µ with the smallest size capable of coalescing to all k-partitions of n.
28

Taking a straightforward example, consider n = 7 and k = 3. The solution to the minimum size
coalescing partition problem is the partition µ = {3, 2, 1, 1}. This partition demonstrates the ability to
coalesce to all eight k-partitions of the positive integer 7 for k ≤ 3, as illustrated:
{7} = {3 + 2 + 1 + 1}, {6, 1} = {3 + 2 + 1, 1}, {5, 2} = {3 + 2, 1 + 1}, {4, 3} = {3 + 1, 2 + 1},
{5, 1, 1} = {3 + 2, 1, 1}, {4, 2, 1} = {3 + 1, 2, 1}, {3, 3, 1} = {3, 2 + 1, 1}, {3, 2, 2} = {3, 2, 1 + 1}
The minimal nature of µ is evident since it only has 4 elements, indicating that it is not possible to coalesce
to all 3-partitions of n = 7 using a partition with fewer elements.
Now, envision n representing the total demand of a customer for goods delivery in practical demand
fulfillment contexts, with k indicating the maximum fleet size. This represents the maximum number of
divisions into which a customer’s demand can be segmented. The investigation into the minimum size
coalescing partition problem thus becomes crucial for formulating a lossless a priori splitting rule tailored
to the split-delivery vehicle routing problem.
The methodology for constructing the partition µ involves systematically breaking down the demand
n into segments in the following manner: Initiate with µ1, determined by the ceiling of n divided by k.
For each additional element, calculate it as the ceiling of the leftover demand divided by k. The greedy
algorithm proceeds until the aggregate of all {µi} is obtained, effectively coalescing to all k-partitions of
n.
This chapter aims to validate the proposed solution’s optimality, commencing with the affirmation
of its feasibility through pivotal lemmas. Subsequent discussions will verify that the proposed method
could achieve a partition of minimal size while preserving the lossless property. This strategy will be
further adapted and utilized to effectively address the split-delivery vehicle routing problem and its variant
incorporating time windows.
29

First, it is evident from the pigeonhole principle that for every partition λ ⊢ n with |λ| ≤ k, there
exists an index i for which λi ≥ ⌈n/k⌉. Consequently, if the partition {µ1, . . . , µl} coalesces to all kpartitions of (n − ⌈n/k⌉), then the completed set {⌈n/k⌉, µ1, . . . , µl} must coalesce to all k-partitions of
n. Expanding upon this notion, the item ⌈n/k⌉ can invariably be a member of a set that coalesces to all
k-partitions of n, for any specified value of n.
It is important to note that ⌈n/k⌉ represents the uppermost value for which this claim holds true,
as the selection of any larger value would inhibit coalescing to partitions where ⌈n/k⌉ is the maximum
value. This validates the feasibility of the suggested approach. The next phase involves demonstrating
that this solution is indeed optimal, which will be based on the specific properties of feasible solutions to
the minimum size coalescing partition problem.
Lemma 7. Consider the scenario where the partition λ coalesces to all k-partitions of n, and the partition γ
is defined such that γ ⊢ n and |γ| ≤ k. Assume without loss of generality that λ = {λ1 ≥ λ2 ≥ · · · ≥ λl},
and γ = {γ1 ≥ γ2 ≥ · · · ≥ γr}. The construction of γ from λ can be described through the following greedy
algorithm: First, include λ1 in the sum yielding γ1. Then, for each subsequent j ∈ {2, . . . , l}, iteratively
include λj in the sum that forms the largest remaining γi which can still accommodate λj .
Proof. Assuming, for the purpose of contradiction, that in the course of constructing γ, the greedy algorithm encounters a situation where it is unable to incorporate the value λj
∗ into any ri
. Utilizing the
disjoint set decomposition ∪
|γ|
i=1Si = {1, . . . , j∗ − 1} to depict the stage immediately prior to the addition
of λj
∗ , the scenario can be described as
γi =


X
j∈Si
λj

 + αi ∀i ∈ {1, . . . , |γ|},
30

with αi representing the remaining capacity in γi to be filled. For a contradiction to manifest, the following
condition must be satisfied:
αi < λj
∗ ≤ λj
∗−1 ≤ . . . ≤ λ1 ∀i ∈ {1, . . . , |γ|}.
Meanwhile, consider that λ can coalesce to any k-partition of n, including the partition
α = {α|γ|
, α|γ|−1
, · · · , α2, n −
X
|γ|
i=2
αi}.
Given that every αi
is strictly smaller than λj
∗ , when coalescing λ to form the partition α, elements greater
than or equal to λj
∗ cannot be used in sums yielding any αi
. Thus, all λj for j ∈ {1, 2, . . . , j∗} must be
dedicated to forming the last component of α, necessitating
n −
X
|γ|
i=2
αi ≥


X
|γ|
i=1
X
j∈Si
λj

 + λj
∗ .
Also, by incorporating the property of γi
,
n =


X
|γ|
i=1
X
j∈Si
λj

 +
X
|γ|
i=1
αi
,
n −
X
|γ|
i=2
αi =


X
|γ|
i=1
X
j∈Si
λj

 + α1.
Through the integration of the above two arguments, it can be inferred that
λj
∗ ≤ α1,
leading to an undeniable contradiction.
31

Hence, in the context of constructing γ according to Lemma 7, for each j ∈ {1, 2, . . . , |λ|}, it is always
possible to fit λj into a particular γi without ever exceeding it. This establishes the soundness of Lemma
7.
For example, by the proposed greedy algorithm, the partition {5, 4, 2, 2, 1, 1} coalesces to all 3-partitions
of 15. Specifically, to form the partition {6, 4, 3, 2}, the decomposition guaranteed by Lemma 7 can be expressed as
{6, 4, 3, 2} = {5 + 1, 4, 2 + 1, 2}.
Consider an arbitrary partition λ = {λ1 ≥ λ2 ≥ · · · ≥ λl}, which coalesces to all k-partitions of n. Let
ν be a partition of (n − λ1), with |ν| ≤ k. It is possible to construct another partition γ = {γ1 ≥ γ2 ≥
· · · ≥ γ|ν|} of n by setting
{γ1, γ2, . . . , γ|ν|} = {λ1 + ν1, ν2, . . . , ν|ν|},
where γ1 = λ1 + ν1 is the largest element of γ. Drawing from the conclusion of Lemma 7, γ can be
constructed by coalescing λ, starting with the inclusion of λ1 to form γ1, and continuing with the orderly
inclusion of λ2, . . . , λl
. Advancing this notion, the set {λ2, . . . , λl} coalesces to any specified partition
ν, thereby encompassing all k-partitions of (n − λ1). Meanwhile, for any m ∈ {2, . . . , l}, the subset
{λm, . . . , λl} coalesces to all k-partitions of (n−
Pm−1
j=1 λj ). The proof of optimality for the main theorem
can then be outlined as follows:
Theorem 8. Let n and k be integers, with k ≤ n. Construct the partition µ beginning with µ1 = ⌈n/k⌉.
For each i > 1, define µi = ⌈(n −
Pi−1
j=1 µj )/k⌉. Continue the process until no component remains. Then, µ
constitutes a solution to the minimum size coalescing partition problem (MSCP). That is, µ is a minimum size
partition that coalesces to all k-partitions of n.
32

Proof. Let λ be a partition that coalesces to all k-partitions of n. Suppose, for contradiction, there exists
such a λ with |λ| < |µ|. This implies:
X
|λ|
i=1
λi = n > X
|λ|
i=1
µi
.
However, it is argued that
Xm
i=1
λi ≤
Xm
i=1
µi ∀m ∈ {1, . . . , min(|µ|, |λ|)}, (4.1)
resulting in a contradiction with the previous statement. To substantiate Equation (4.1) via induction, first
consider m = 1. If λ1 > µ1, then λ1 > ⌈n/k⌉, suggesting λ1 cannot be an element of any partition γ
where γ1 = ⌈n/k⌉, rendering λ infeasible.
Assuming Equation (4.1) applies for m − 1, and considering the fact that λm is the largest element in
a partition coalescing to all k-partitions of (n −
Pm−1
i=1 λi), it follows that
λm ≤
&
(n −
mX−1
i=1
λi)/k'
.
33

Accordingly,
Xm
i=1
λi =
mX−1
i=1
λi
!
+ λm ≤
mX−1
i=1
λi
!
+
& n −
mX−1
i=1
λi
!
k
'
=
& mX−1
i=1
λi
!
+

n −
mX−1
i=1
λi
!
k
'
=
&
n
k
+
k − 1
k
mX−1
i=1
λi
'
≤
&
n
k
+
k − 1
k
mX−1
i=1
µi
'
=
mX−1
i=1
µi
!
+
& n −
mX−1
i=1
µi
!
k
'
=
mX−1
i=1
µi
!
+ µm =
Xm
i=1
µi
.
The induction-based proof of Equation (4.1) is thereby concluded. This completes the demonstration of
optimality for the partition µ, which is constructed through the greedy algorithm described in Theorem
8.
Now let |µ| = L(n, k) denote the size of the partition µ generated by the greedy algorithm defined in
Theorem 8, for any given values of n and k. While a closed-form expression for L(n, k) is not provided, it
is evident that the following recurrence relation is valid:
L(n, k) = 1 + L
k − 1
k
n

, k
, L(1, k) = 1.
Additionally, let L
′
(n, k) be defined to meet the recurrence relation:
L
′
(n, k) = 1 + L

k − 1
k
n, k
, L′
(1, k) = 1.
34

For any specified values of n and k, L
′
(n, k) ≥ L(n, k). Moreover, L
′
(n, k) is increasing in n, letting
|µ| = L(n, k) ≤ L
′
(n, k) ≤ L
′

k
k − 1
⌈logk/(k−1)(n)⌉
, k!
= ⌈logk/(k−1)(n)⌉ + 1.
Therefore, an upper bound on the size of the minimum coalescing partition can be formulated as
|µ| ≤ ⌈logk/(k−1)(n)⌉ + 1.
Note that this bound weakens as k increases relative to n. Nonetheless, the exact size of |µ| for specific
values of n and k can be directly determined by executing the algorithm. Table 4.1 presents the calculated
values of |µ| for n ranging from 1 to 22 and k from 1 to 11.
n
k
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22
1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
2 1 2 2 3 3 3 3 4 4 4 4 4 4 4 4 5 5 5 5 5 5 5
3 1 2 3 3 4 4 4 5 5 5 5 6 6 6 6 6 6 7 7 7 7 7
4 1 2 3 4 4 5 5 6 6 6 7 7 7 7 8 8 8 8 8 9 9 9
5 1 2 3 4 5 5 6 6 7 7 7 8 8 8 9 9 9 9 10 10 10 10
6 1 2 3 4 5 6 6 7 7 8 8 9 9 9 10 10 10 11 11 11 11 12
7 1 2 3 4 5 6 7 7 8 8 9 9 10 10 10 11 11 11 12 12 12 12
8 1 2 3 4 5 6 7 8 8 9 9 10 10 11 11 12 12 12 13 13 13 14
9 1 2 3 4 5 6 7 8 9 9 10 10 11 11 12 12 13 13 13 14 14 14
10 1 2 3 4 5 6 7 8 9 10 10 11 11 12 12 13 13 14 14 15 15 15
11 1 2 3 4 5 6 7 8 9 10 11 11 12 12 13 13 14 14 15 15 16 16
Table 4.1: The size of the minimum coalescing partition µ across various choices of n and k. Fundamentally, µ embodies the partition with the least number of elements that can coalesce to all partitions λ ⊢ n
where |λ| ≤ k.
In addressing the split-delivery vehicle routing problem (SDVRP), the subsequent analysis introduces
a splitting rule that balances between partition size and algorithmic efficiency within a theoretically robust
framework. This splitting strategy simply involves breaking down each customer demand di
into segments
35

according to the minimum size partition that coalesces to all k-partitions of di
, as thoroughly described in
Theorem 8.
This approach enables the segmentation of customer demands across multiple trucks, effectively transforming the SDVRP into a CVRP with an increased number of nodes. Trucks can be allocated to service
distinct combinations of segments from the same customer demand, thereby maintaining all feasible ways
to split demands within the solution space of the SDVRP. By ensuring no viable solutions are disregarded,
the SDVRP can be solved optimally, provided that the CVRP is optimally tackled. An upper bound for the
size of the reduction is illustrated as:
Number of CVRP nodes ≤
Xn
i=1
(⌈logk/(k−1)(di)⌉ + 1).
Yet, observations indicated in Table 4.1 suggest that the actual size of partitions tends to be significantly
smaller.
In real-world scenarios, it is reasonable to anticipate that demand splitting among trucks is relatively
limited. In other words, it is probable that for any given customer, the optimal SDVRP solution would
involve no more than some ¯k < k trucks. For example, even with access to a large fleet of k = 100
trucks, it is unlikely that any customer’s demand necessitates servicing from more than, say, ¯k = 3 trucks.
Thus, rather than requiring the splitting rule to coalesce to all k-partitions, focusing on the minimum size
partition that coalesces to all ¯k-partitions is more pragmatic. The resulting splitting mechanism still allows
for the discovery of viable solutions that adhere to the simplified assumption, while markedly reducing
the complexity of the resulting CVRP problem to be solved.
Moreover, to implement an a priori splitting rule, it is essential for all the demands to be integral, allowing for the division of commodities into discrete units of size 1. For continuous commodities, discretization
36

can be achieved by defining a unit of size 1 as the minimal amount of the commodity, with splitting permitted only in multiples of this unit. Selecting an appropriate level of granularity for this discretization
is crucial. Indeed, opting for a sufficiently small unit size for the commodities does not necessarily lead
to an excessive increase in the size of the generated CVRP problem, following an a priori split. Reducing
the unit size by half effectively doubles each demand di
, which, in turn, increases the size of the minimum
coalescing partition to at most
Xn
i=1
(1 + ⌈logk/(k−1)(2di)⌉) ≤
Xn
i=1
(1 + ⌈logk/(k−1)(di)⌉ + ⌈logk/(k−1)(2)⌉).
Thus, increasing the granularity by a set multiple only contributes a constant to the overall size of the generated CVRP problem, a factor that remains invariant regardless of the initial sizes of customer demands.
37

Chapter 5
Advancements in TSP through Enumeration
5.1 Enumeration of path permutations
To derive a strictly proven better upper bound of β2, the main idea is to substitute the horizontal traversal
strategy outlined in Section 3.1 with an approach that offers greater versatility. This necessitates a systematic methodology for constructing alternative routes that traverse clusters of x-adjacent points (also
defined in Section 3.1) in a more cost-effective manner.
The preferred approach entails using horizontal strips of height h/√
n, selecting x-adjacent (k + 1)-
tuples (Xi
, Yi), . . . ,(Xi+k, Yi+k) for a specified (small) k, and choosing the permutation that minimizes
the path length for those (k + 1)-tuples. It is critical to only consider permutations that start with the
leftmost point (Xi
, Yi) and end with the rightmost point (Xi+k, Yi+k) in each (k + 1)-tuple, ensuring
these paths can be integrated into a cohesive path that traverses all points in the horizontal strip.
To secure an improved upper bound of β2, the primary step is to pinpoint multiple scenarios where
x-adjacent (k + 1)-tuples (Xi
, Yi), . . . ,(Xi+k, Yi+k) exhibit a strictly shorter alternative traversal path
compared to the conventional left-to-right path, which corresponds to the identity permutation. Assuming without loss of generality that Xi ≤ . . . ≤ Xi+k, Figure 5.1 provides an example of a 5-tuple where
the alternative permutation (0,2,3,1,4) delivers the shortest path, outperforming that of the identity permutation (0,1,2,3,4).
38

{
{
{
{
Figure 5.1: A depiction of a 5-tuple (k = 4) where the permutation (0,2,3,1,4) among the six possible
permutations achieves the lowest-cost path.
Figure 5.2: A representation of all six permutations for a 5-tuple (k = 4), highlighting that while the
identity permutation often results in the lowest traversal cost, alternative permutations may also offer
shorter paths.
39

(a) k = 1, 2 (b) k = 3 (c) k = 4
Figure 5.3: The number of traversal options within each horizontal strip increases superexponentially
with k. Additionally, a larger k allows for greater flexibility in selecting the optimal route for traversing
the horizontal strip, as illustrated in 5.3b and 5.3c.
It should be noted that for k = 3, the analysis includes x-adjacent 4-tuples with only two permissible
permutations: (0,1,2,3) or (0,2,1,3). For k = 4, six permutation options are available for 5-tuples: (0,1,2,3,4),
(0,1,3,2,4), (0,2,1,3,4), (0,2,3,1,4), (0,3,1,2,4), and (0,3,2,1,4), as demonstrated in Figure 5.2. Generally, for
(k+ 1)-tuples, there are (k−1)! potential permutations. As k increases, the number of traversal strategies
within each horizontal strip grows superexponentially, as depicted in Figure 5.3.
Sampling a significant number of x-adjacent five-tuples facilitates a quantitative comparison between
the identity permutation (left-to-right traversal) and alternative permutations via simulation. The simulation results, as detailed in Section 2.5 of [94], reveal that with the strip height hyperparameter h = 3,
the theoretical limit of improvement on β2’s upper bound using five-tuples is approximately 0.8902. This
claim is corroborated and expanded upon in broader contexts through our simulations.
In Figure 5.4a, simulations of 5-tuples indicate that approximately 27% of our 30 million samples features a more cost-effective alternative permutation compared to the identity permutation, synonymous
with the conventional left-to-right path. The permutations (0,1,3,2,4) and (0,2,1,3,4) are deemed the most
economical among the six permutations in about 10.5% of the samples each. The permutations (0,2,3,1,4),
(0,3,1,2,4), and (0,3,2,1,4) are identified as the most cost-effective among the six permutations for approximately 2.1%, 2.1%, and 1.7% of the samples, respectively.
40

(a) (b)
Figure 5.4: 5.4a shows simulations for 5-tuples, with the x-axis representing the ratio of the shortest path
out of 6 options to the identity permutation (0,1,2,3,4). 5.4b illustrates simulations for 6-tuples, where the
x-axis indicates the ratio of the shortest path out of 24 options to the identity permutation (0,1,2,3,4,5). The
juxtaposition suggests that an increase in k diminishes the likelihood of the identity permutation being
the most cost-effective for achieving the shortest path.
In Figure 5.4b, simulations involving 6-tuples demonstrate that in over 34% of our 30 million samples, an
alternative path proves more beneficial than the left-to-right path. This result is expected, as the addition of
more points to the tuple increases the number of permutation options, indicating that a rise in k leads to a
higher proportion of samples where a non-identity permutation emerges as the optimal traversal strategy.
Following the establishment that alternative permutations are beneficial in numerous instances, the
next critical step is to provide a rigorous and mathematical justification that these instances of reduced
path lengths contribute positively towards the improvement of the upper bound of β2.
The classic Beardwood-Halton-Hammersley upper bound, as introduced by [19], is applicable to the
scenario of k = 1 (or k = 2, considering predetermined endpoints). The original bound described in (3.1)
is fundamentally
β2 ≤
1
h
E



Z
h
2
(U0 − U1)



,
41

where the adoption of k = 3 leads to a potentially tighter upper bound
β2 ≤
1
3h
E min




Z1
h
2
(U0 − U1)

+

Z2
h
2
(U1 − U2)

+

Z3
h
2
(U2 − U3)

Z1 + Z2
h
2
(U0 − U2)

+

Z2
h
2
(U2 − U1)

+

Z2 + Z3
h
2
(U1 − U3)




.
Consistent with the framework outlined in Section 3.1, all Zi variables are independent and adhere to an
exponential distribution with a mean of 1, while all Ui variables are independent and follow a uniform
distribution. The original upper bound presented by [19] employs a simple left-to-right traversal strategy.
A more restrictive upper bound for 4-tuples (k = 3) is achievable by permitting an exchange between the
two points in the middle.
To formalize, the refined upper bound for scenarios involving 5-tuples (k = 4) can also be detailed as
follows:
β2 ≤
1
4h
E min




Z1
h
2
(U0 − U1)

+

Z2
h
2
(U1 − U2)

+

Z3
h
2
(U2 − U3)

+

Z4
h
2
(U3 − U4)

Z1
h
2
(U0 − U1)

+

Z2 + Z3
h
2
(U1 − U3)

+

Z3
h
2
(U3 − U2)

+

Z3 + Z4
h
2
(U2 − U4)

Z1 + Z2
h
2
(U0 − U2)

+

Z2
h
2
(U2 − U1)

+

Z2 + Z3
h
2
(U1 − U3)

+

Z4
h
2
(U3 − U4)

Z1 + Z2
h
2
(U0 − U2)

+

Z3
h
2
(U2 − U3)

+

Z2 + Z3
h
2
(U3 − U1)

+

Z2 + Z3 + Z4
h
2
(U1 − U4)

Z1 + Z2 + Z3
h
2
(U0 − U3)

+

Z2 + Z3
h
2
(U3 − U1)

+

Z2
h
2
(U1 − U2)

+

Z3 + Z4
h
2
(U2 − U4)

Z1 + Z2 + Z3
h
2
(U0 − U3)

+

Z3
h
2
(U3 − U2)

+

Z2
h
2
(U2 − U1)

+

Z2 + Z3 + Z4
h
2
(U1 − U4)




The subsequent sections in this chapter are committed to the rigorous derivation of bounds for integrals
derived from the aforementioned expectations. The primary hurdles include the high-dimensional nature
42

of these integrals, as a (k + 1)-tuple requires 2k + 1 variables, and the absence of closed-form solutions
for these integrals. The next section will introduce techniques to streamline these integrals. The focus of
this dissertation will predominantly be on the analysis of 5-tuples, as this configuration offers an optimal
compromise between improving the upper bound of β2 and managing computational time and memory
resources efficiently.
5.2 Construction of high-dimensional integral
The subsequent analysis involves reverting variables to the original entries, Xi
, by establishing Xi =
Pi
j=1 Zi for i ∈ {1, . . . , k}, and designating X0 = 0 for convenience in notation. The bounds delineated
in the preceding section can now be reformulated for any (k + 1)-tuple through the following expression:
β2 ≤
1
khE min
π∈Πk
X
k
i=1

Xπ(i) − Xπ(i−1)
h
2
(Uπ(i) − Uπ(i−1))

, (5.1)
where Πk encompasses all permutations π of the set {0, . . . , k} with the conditions π(0) = 0 and π(k) = k.
To compute the expectation, it is necessary to convert it into an integral form, thereby necessitating the
joint probability density function for (X1, . . . , Xk). This derivation is fortunately direct.
Lemma 9. The joint density function of (X1, . . . , Xk) is f(x1, . . . , xk) = e
−xk for 0 ≤ x1 ≤ · · · ≤ xk.
Proof. Although this derivation is commonly encountered as an introductory exercise in the context of exponential random variables and their attributes, an inductive proof is provided here for thoroughness. It is
presupposed that the proposition is valid for(X1, . . . , Xk−1). The joint distribution for(X1, . . . , Xk−1, Xk)
is
f(x1, . . . , xk−1)g(xk|Xk−1 = xk−1),
43

where g represents the conditional pdf for xk. Given that Xk = Xk−1 + Zk,
g(xk|Xk−1 = xk−1) = e
−(xk−xk−1)
for xk ≥ xk−1 ,
from which the result follows.
(5.1) can consequently be articulated as the primary integral below:
β2 ≤
1
kh Z
D
e
−xk min
π∈Πk
X
k
i=1

xπ(i) − xπ(i−1)
h
2
(uπ(i) − uπ(i−1))

dV , (5.2)
where
D = {(x1, . . . , xk, u0, . . . , uk) : 0 ≤ x1 ≤ · · · ≤ xk , 0 ≤ ui ≤ 1} .
The ensuing sections of this chapter are dedicated to meticulously establishing upper bounds for the integral delineated in (5.2). For simplicity in notation, the integrand gπ = gπ(x1, . . . , xk, u0, . . . , uk) is defined
as
gπ =
e
−xk
kh
X
k
i=1

xπ(i) − xπ(i−1)
h
2
(uπ(i) − uπ(i−1))

,
thus rendering (5.2) synonymous with writing
β2 ≤
Z
D
min
π∈Πk
gπ dV .
It is pertinent to note that the term e
−xk can be placed either inside or outside the minπ∈Πk
expression, in
light of the assumption that π(k) = k.
To improve the original upper bound described in Section 3.1, consider a high-dimensional box (hereinafter referred to as a hyperrectangle) B ⊂ D, within which a certain permutation π
∗
B
, distinct from the
44

identity permutation, is favored. In the context of a (k+1)-tuple, B is designed to accommodate a (4k+2)
dimensional input (denoting the bounds of the hyperrectangle) and can be delineated as follows:
B([s1, t1], . . . , [sk, tk], [a0, b0], . . . , [ak, bk]) = {(x1, . . . , xk, u0, . . . , uk) : ai ≤ ui ≤ bi
, sj ≤ xj ≤ tj} ,
where
0 ≤ s1 ≤ t1 ≤ . . . ≤ sk−1 ≤ tk−1 ≤ sk ≤ tk, 0 ≤ ai ≤ bi ≤ 1 .
It is important to clarify that π
∗
B
does not have to be the superior permutation throughout every part of B:
a general preference within B for π
∗
B
over the identity permutation is adequate. This high-dimensional
box facilitates a small but positive advancement on the original upper bound established by [19] through
the following formulation:
β2 ≤
Z
D
min
π∈Πk
gπ dV
=
Z
D\B
min
π∈Πk
gπ dV +
Z
B
min
π∈Πk
gπ dV
≤
Z
D\B
gπId dV +
Z
B
gπ
∗
B
dV .
Here, π = πId represents the identity permutation. As a first step, it is assumed that the identity permutation is utilized for traversal outside B, with a shift to π
∗
B
for traversal within B. The domain D \ B is
cumbersome, yet simplification is attainable by representing the integral as the deduction of the integral
over B from the integral over D. Specifically,
Z
D\B
gπId dV =
Z
D
gπId dV
| {z }
(∗)
−
Z
B
gπId dV .
45

Moreover, the integral over D, denoted by (∗), aligns with the original upper bound introduced by [19],
which consistently applies the identity permutation throughout its evaluation. This is mathematically
depicted as:
Z
D
gπId dV =
1
kh Z
D
e
−xk X
k
i=1

xi − xi−1
h
2
(ui − ui−1)

dV
=
1
kh Z ∞
xk=0
· · · Z x2
x1=0
Z 1
uk=0
· · · Z 1
u0=0
e
−xk X
k
i=1

xi − xi−1
h
2
(ui − ui−1)

du0 · · · duk dx1 · · · dxk
=
1
3h
5
Z ∞
0
e
−z

3h
2
z
2
log(h
2
/z +
p
h
4/z2 + 1) + 2z
3 + (h
4 − 2z
2
)
p
h
4 + z
2

dz ,
which mirrors the expression found in (3.1), and is an established figure for any given h. Adopting h =
√
3,
as in the initial upper bound (i.e., 0.92116027) from [19], this methodology is poised to offer a modest
yet definitive improvement, based on the premise that π
∗
B
surpasses the identity permutation on average
within the hyperrectangle B.
Consequently, for a hyperrectangle B and the corresponding best permutation π
∗
B
, the refined upper
bound for the Euclidean travelling salesman constant β2 can be succinctly articulated as
β2 ≤
Z
D
gπId dV −
Z
B
gπId dV +
Z
B
gπ
∗
B
dV ,
from which the quantity of interest can be represented as
ϵ(B, π∗
B) = Z
B
gπId dV −
Z
B
gπ
∗
B
dV . (5.3)
46

In this context, π
∗
B
refers to the permutation that minimizes the average traversal over hyperrectangle B.
The expression ϵ(B, π∗
B
) is henceforth termed the contribution of hyperrectangle B, indicating a mandatory nonnegative value. It is important to recognize the possibility of selecting the identity permutation
as π
∗
B
, i.e., π
∗
B = πId. Under such circumstances, the contribution of hyperrectangle B equates to zero.
5.3 Separation of the integral
The discussion in the preceding section validates the feasibility of improving the upper bound established
by Beardwood et al. [19], as outlined in Section 3.1. This improvement is achieved in scenarios where
there exists a hyperrectangle B in which an alternative permutation π
∗
B
essentially yields a lower average
traversal cost compared to the identity permutation.
The incremental improvement of the upper bound documented in this dissertation results from the
aggregate impact of incorporating over a vast, non-overlapping assortment of such hyperrectangles and
systematically implementing the methodologies described earlier. For the establishment of guaranteed
bounds, it is critical to verify the manageability of integration across B.
Continuing with the premise that π
∗
B
denotes the optimal permutation while minimizing R
B
gπ dV ,
and leveraging the linearity property of integration, the expansion of terms facilitates the subsequent
representation:
Z
B
gπ dV =
1
kh Z
B
e
−xk X
k
i=1

xπ(i) − xπ(i−1)
h
2
(uπ(i) − uπ(i−1))

dV
=
1
kh
X
k
i=1
Z
B
e
−xk

xπ(i) − xπ(i−1)
h
2
(uπ(i) − uπ(i−1))

dV . (5.4)
This allows the integral over B to decompose into a sum of integrals, each encompassing no more than 5
variables, specifically xk, xπ(i)
, xπ(i−1), uπ(i)
, and uπ(i−1). For the last summand with i = k, the variable
47

count reduces to 4 as it is predetermined that π(k) = k. Contrarily, the first k − 1 summands, which
exclude xk from the | · | term, can be simplified as follows:
C
Z t1
t0
Z d2
c2
Z d1
c1
Z b2
a2
Z b1
a1
e
−s

x2 − x1
h
2
(u2 − u1)

du1 du2 dx1 dx2 ds
=C
′
Z d
′
2
c
′
2
Z d
′
1
c
′
1
Z b2
a2
Z b1
a1
p
(x2 − x1)
2 + (u2 − u1)
2 du1 du2 dx1 dx2 ,
where C denotes a constant derived by integrating out the variables absent in the integrand, and C
′
is
further obtained by integrating out the e
−s
term, removing h through a variable substitution in u1 and
u2, and letting c
′
i = h
2
ci
, d
′
i = h
2di
. Remarkably, it is discovered that this integral can be expressed in a
closed form:
Lemma 10. We have
Z d2
c2
Z b2
a2
Z d1
c1
Z b1
a1
p
(x2 − x1)
2 + (u2 − u1)
2 du1 dx1 du2 dx2 = F(x1, u1, x2, u2)|
b1
u1=a1

d1
x1=c1

b2
u2=a2

d2
x2=c2
where
F(x1, u1, x2, u2) = Z Z Z Z p
(x2 − x1)
2 + (u2 − u1)
2 du1 dx1 du2 dx2
= −
1
288
(3x
4
2 + 4x
2
2u
2
2 + 3x2u
3
2 − 12x
3
2x1 − 8x2u
2
2x1 − 3u
3
2x1 + 18x
2
2x
2
1 − 12x2x
3
1 − 12x
3
2u1
− 28x2u
2
2u1 + 12u
2
2x1u1 + 42x2u2u
2
1 − 18u2x1u
2
1 − 28x2u
3
1 + 12x1u
3
1
)
−
1
24
(u
2
2 − 2u2u1 + 2u
2
1
)x2(u2 − 2u1)u2 log 2
+
1
24
[(x1 − x2)
2 + (u1 − u2)
2 + x2u2 − u2x1 − x2u1 + x1u1](x2 − x1 + u1 − u2)
(x2 − x1)(u2 − u1) log(p
(x1 − x2)
2 + (u1 − u2)
2 − u1 + u2)
+
1
12
(x2 − x1)(u2 − u1)
4
log

p
(x1 − x2)
2 + (u1 − u2)
2 − x1 − u1 + x2 + u2

−
1
60
[(x1 − x2)
2 − (u1 − u2)
2 + x2u2 − u2x1 − x2u1 + x1u1]
[(x1 − x2)
2 − (u1 − u2)
2 − x2u2 + u2x1 + x2u1 − x1u1]
p
(x1 − x2)
2 + (u1 − u2)
2
48

Proof. This can be verified through the direct computation of the antiderivative function F, utilizing a
computer algebra system [44].
For the last summand of the integral mentioned in (5.4) where i = k, the absence of a closed-form
antiderivative is noted due to the exponential term. Nonetheless, we can still reduce it to a univariate
integral by performing integration over the other three variables, as demonstrated by the subsequent result:
Lemma 11. We have
Z dk
ck
Z bk
ak
Z b1
a1
Z d1
c1
e
−xk
p
(xk − x1)
2 + (uk − u1)
2 dx1 du1 duk dxk
=
Z dk
ck
e
−xk
Z bk
ak
Z b1
a1
Z d1
c1
p
(xk − x1)
2 + (uk − u1)
2 dx1 du1 duk

dxk =
Z dk
ck
e
−xkG(xk) dxk ,
where
G(xk) = F(x1, u1, xk, uk)|
bk
uk=ak

b1
u1=a1

d1
x1=c1
F(x1, u1, xk, uk) = Z Z Z p
(xk − x1)
2 + (uk − u1)
2 dx1 du1 duk
=
1
6
(xk − x1)
3u1 log p
(x1 − xk)
2 + (u1 − uk)
2 + u1 − uk

+
1
6
(7x
2
k − 2xkx1 + x
2
1
)(xk − x1)uk log p
(x1 − xk)
2 + (u1 − uk)
2 − u1 + uk

+
1
24
(12x
2
k + u
2
k
)u
2
k
log

p
(x1 − xk)
2 + (u1 − uk)
2 − x1 + xk − u1 + uk

−
1
24
(12x
2
k + u
2
k
)u
2
k
log

p
(x1 − xk)
2 + (u1 − uk)
2 + x1 − xk − u1 + uk

+
1
12
(2u
2
k − 3uku1 + 2u
2
1
)uku1 log p
(xk − x1)
2 + (u1 − uk)
2 − xk + x1
|u1 − uk|
!
−
1
24

12x
2
k − u
2
1

u
2
1
log

p
(x1 − xk)
2 + (u1 − uk)
2 − x1 + xk + u1 − uk

+
1
24

12x
2
k − u
2
1

u
2
1
log

p
(x1 − xk)
2 + (u1 − uk)
2 + x1 − xk + u1 − uk

−
1
24
(14x
2
k − 3u
2
k − 4xkx1 + 2x
2
1 + 6uku1 − 3u
2
1
)(xk − x1)
p
(x1 − xk)
2 + (u1 − uk)
2

Proof. Again, we used FriCAS [44] to compute the antiderivative F.
The absence of a closed-form solution for the preceding integral necessitates the use of numerical
methods, potentially introducing approximation errors that could compromise the reliability of our analysis. Fortunately, the inherent structure of the problem presents an opportunity to establish true rigorous
upper and lower bounds, thereby preserving the integrity of our findings.
Remark 12. The function G(xk) exhibits convexity, a characteristic derived from the integration of
p
(xk − x1)
2 + (uk − u1)
2 ,
which is itself convex. It is also monotonically increasing in xk, given that xk ≥ x1. Consequently, the integrand e
−xkG(xk) in Lemma 11 is the product of two convex functions. We can bound R dk
ck
e
−xkG(xk) dxk
from above by bounding the functions e
−xk and G(xk)from above individually with piecewise linear functions, namely the trapezoidal rule. Taking their product would result in a piecewise quadratic function.
Remark 13. In the context of evaluating a hyperrectangle B, the objective is to ascertain whether a specific
permutation π
∗
B
offers advantages over the identity permutation πId. Contrary to Remark 12, establishing a
lower bound for R dk
ck
e
−xkG(xk) dxk using the identity permutation becomes crucial. If the upper bounding
integral obtained from π
∗
B
outperforms the lower bounding integral obtained from πId, it confirms the
superiority of the permutation π
∗
B
. The lower bound of R dk
ck
e
−xkG(xk) dxk can be derived by bounding
e
−xk and G(xk) from below individually with piecewise linear functions, specifically tangent lines, and
their subsequent product yields a piecewise quadratic function.
Figure 5.5 graphically demonstrates the concepts introduced in Remark 12 and 13. The solid red, blue,
and purple curves depict the functions e
−xk , G(xk), and their product e
−xkG(xk), respectively. The red
and blue lines signify the piecewise linear bounds for e
−xk and G(xk), while the purple lines indicate
50

Figure 5.5: The piecewise quadratic bounding method is utilized to determine reliable lower and upper
bounds for the final summand in (5.4), applicable to any chosen valid hyperrectangle B. The trapezoidal
rule is employed for determining the upper bounds of G(xk) and e
−xk , and the midpoint rule is adopted
for calculating their lower bounds. Importantly, during computational experiments, the lower bounding
technique is exclusively applied to the identity permutation, while the upper bounding technique is reserved for assessing all alternative permutations. The capability to secure positive contributions under
such stringent criteria unequivocally substantiates the improvement of β2’s upper bound.
51

the piecewise quadratic bounds for their product. The dashed and dotted lines represent upper and lower
bounds, correspondingly.
An expanded explanation of how these bounds are constructed in our computational experiments is
provided as follows. The domain of G(xk), specified as [ck, dk], is divided evenly into M segments. The
endpoints of each segment are represented as sk = p0 ≤ p1 ≤ p2 ≤ . . . ≤ pM−1 ≤ pM = tk, where
pi = sk + i(tk − sk)/M. In Figure 5.5, M is intentionally kept small to enhance the visual distinction
between upper and lower bounds. In actual computations, M is set to 100.
Given the convex nature of G(xk) within its domain, a secant line connecting points (pi
, G(pi)) and
(pi+1, G(pi+1))serves as an overestimate for G(xk) across[pi
, pi+1]. Conversely, a tangent line that passes
through the midpoint of [pi
, pi+1] provides an underestimate for G(xk) on [pi
, pi+1]. These secant and
tangent lines are denoted as a
(1)
i
xk +b
(1)
i
and a
(2)
i
xk +b
(2)
i
, respectively. For e
−xk , the overestimation and
underestimation on the interval [pi
, pi+1] are represented by c
(1)
i
xk + d
(1)
i
and c
(2)
i
xk + d
(2)
i
, respectively.
A quadratic upper bound for e
−xkG(xk), expressed as (a
(1)
i
xk + b
(1)
i
)(c
(1)
i
xk + d
(1)
i
) for any xk ∈
[pi
, pi+1], is thereby established. Similarly, a quadratic lower bound is determined as(a
(2)
i
xk+b
(2)
i
)(c
(2)
i
xk+
d
(2)
i
). The quadratic bounds, illustrated in Figure 5.5 through purple dashed and dotted lines, effectively
demonstrate the method’s capability in delivering reliable, guaranteed bounds for the target integral. Furthermore, this quadratic bounding technique has been shown to yield minimal empirical relative differences in comparison to Gaussian quadrature estimates, with discrepancies reaching magnitudes of no more
than 10−8
.
To summarize, this chapter lays down the foundational elements necessary for refining the upper
bound of β2 through systematic integration across a series of non-overlapping hyperrectangles B. The
next section will delve into the process of identifying these hyperrectangles and resolving any challenges
that might impede the advancement of the upper bound refinement.
52

5.4 Algorithms
This section provides an in-depth explanation of the functions contained in the Python script helper.py.
As discussed earlier, k = 4 is set universally in order to balance computational efficiency with prospective
solution quality. The primary purpose of the procedures outlined in helper.py is to compute the guaranteed contribution ϵ(B, πB) once the boundaries of the hyperrectangle B and the permutation index
πB are known. This framework is essential for accurately assessing the impact of specific hyperrectangle
configurations on the overall upper bound improvement of β2.
Algorithm 1 describes the process for calculating integrals starting from antiderivatives, with the antiderivative computation facilitated by FriCAS [44]. Note that the procedure F4 addresses various limiting
cases of Lemma 10: when x1 = x2, the logarithmic term within the given expression becomes undefined,
prompting the need for specialized considerations for scenarios where u1 < u2, u1 > u2, or u1 = u2.
To address these, the capabilities of SageMath [83] are employed to articulate the antiderivatives corresponding to these limiting cases. For the general case where x1 ̸= x2, the expression from Lemma 10 is
applicable without modifications. Subsequently, the procedure x4 is tasked with computing the definite
integral within specified bounds, leveraging the antiderivative obtained from F4. This method is tailored
to address the first, second, and third summands of (5.4) as k = 4, enabling the direct computation of their
closed-form expressions.
Algorithm 2 adopts a methodology akin to that of Algorithm 1, albeit focusing on the evaluation of
the final (k
th) integral summand as depicted in (5.4). The procedure F4_x4out addresses different limiting cases of Lemma 11, including situations where x1 = x4 with u1 and u4 either greater than, less
than, or equal to each other. In these limiting scenarios, SageMath [83] is again utilized to articulate the
antiderivatives, ensuring precise representations. For the general case where x1 ̸= x4 and u1 ̸= u4, the
antiderivative remains unaltered. Consequently, the procedure x4_x4out capitalizes on the antiderivative
53

Algorithm 1 Calculation of integral using FriCAS (1st/2nd/3rd summand)
1: procedure F4 (x1, x2, u1, u2: float) ▷ See Lemma 10
2: if x1 = x2 and u1 < u2 then
3: return F(x1, u1, x2, u2; h) for limiting case 1
4: else if x1 = x2 and u1 > u2 then
5: return F(x1, u1, x2, u2; h) for limiting case 2
6: else if x1 = x2 and u1 = u2 then
7: return F(x1, u1, x2, u2; h) for limiting case 3
8: else
9: return F(x1, u1, x2, u2; h) for the general case ▷ Representation of antiderivative
10: procedure x4 (x1, x2, u1, u2: List[float]) ▷ Lower and upper bounds for the integral
11: call F4, Fxyuv_4 ▷ Antiderivative F4 into integral
12: return Fxyuv_4() ▷ 1st/2nd/3rd summand of (5.4) as float
Algorithm 2 Calculation of integral using FriCAS (last summand)
1: procedure F4_x4out (x1, x4, u1, u4: float) ▷ See Lemma 11
2: if x1 = x4 and u1 < u4 then
3: return F(x1, u1, x4, u4; h) for limiting case 1
4: else if x1 = x4 and u1 > u4 then
5: return F(x1, u1, x4, u4; h) for limiting case 2
6: else if u1 = u4 then
7: return F(x1, u1, x4, u4; h) for limiting case 3
8: else
9: return F(x1, u1, x4, u4; h) for the general case ▷ Representation of antiderivative
10: procedure x4_x4out(x1, x4, u1, u4: List[float]) ▷ Lower and upper bounds for the integral
11: call F4_x4out, Fxuv_4out ▷ Antiderivative F4_x4out into integral
12: return Fxuv_4out(x4) ▷ Last summand of (5.4) as a function of x4
13: procedure x4_tangent(x1, x4, u1, u4: List[float]) ▷ Lower and upper bounds for the integral
14: call x4_x4out ▷ Apply the midpoint rule
15: return underestimation ▷ See Figure 5.5 in Section 5.3
16: procedure x4_secant(x1, x4, u1, u4: List[float]) ▷ Lower and upper bounds for the integral
17: call x4_x4out ▷ Apply the trapezoidal rule
18: return overestimation ▷ See Figure 5.5 in Section 5.3
54

provided by F4_x4out to execute the computation of definite integrals within the designated bounds of
B. This step distinctly defines the univariate function G(x4), in line with Lemma 11. It is crucial to note
that the function G(x4) only forms the primary component of the integrand for the final summand of (5.4),
complemented by an additional element e
−x4
.
In the concluding phase, the procedures x4_tangent and x4_secant enrich the analysis by calculating
the guaranteed lower and upper bounds for the k
th and last integral summand. Specifically, x4_tangent
applies the midpoint rule, while x4_secant implements the trapezoidal rule for numerical analysis. The
calculation of numerical bounds is fundamental for resolving the integral mentioned in Lemma 11, especially when closed-form solutions are unattainable, and numerical methods offer a viable alternative. The
practicality of these numerical approaches is visually demonstrated in Figure 5.5 of Section 5.3.
Algorithm 3 presents a structured approach to calculating the expected path lengths for all feasible
permutations of 5-tuples, given any valid hyperrectangle B. Each procedure of the algorithm, such as
x4_0101, x4_0202, x4_0303, and so forth, is tasked with computing the expected length of a specific segment, leveraging the boundaries of B. These computations utilize scalars CB,·,·
, detailed in the Appendix,
which adjust the integral calculations to account for the unique geometry of B.
The algorithm effectively distinguishes between straightforward segment length computations (x4_0101
through x4_2323) and those necessitating bound estimations (x4_UB_1414 through x4_UB_3434, as
well as x4_LB_3434). The latter employs numerical methods, such as the trapezoidal rule for overestimation and the midpoint rule for underestimation, to approximate the expected path lengths, particularly
when the last summand lacks a closed-form solution.
The final set of procedures (calc_01234_LB, as well as calc_01324_UB through calc_03214_UB)
aggregate these calculated expected segment lengths to assess the overall path length across all six feasible permutations. Notably, underestimation applies to the identity permutation, whereas overestimation
pertains to all other permutations, providing a reliable upper bound for the Euclidean travelling salesman
55

Algorithm 3 Calculation of expected path lengths for all possible permutations, given B
1: procedure x4_0101(B: 2DArray) ▷ Input all the boundaries of hyperrectangle B
2: return CB,0,1 x4(x0, x1, u0, u1) ▷ See Appendix for CB,·,·
3: procedure x4_0202(B: 2DArray)
4: return CB,0,2 x4(x0, x2, u0, u2)
5: procedure x4_0303(B: 2DArray)
6: return CB,0,3 x4(x0, x3, u0, u3)
7: procedure x4_1212(B: 2DArray)
8: return CB,1,2 x4(x1, x2, u1, u2)
9: procedure x4_1313(B: 2DArray)
10: return CB,1,3 x4(x1, x3, u1, u3)
11: procedure x4_2323(B: 2DArray)
12: return CB,2,3 x4(x2, x3, u2, u3)
13: procedure x4_UB_1414(B: 2DArray) ▷ Calculate overestimation
14: return CB,1,4 x4_secant(x1, x4, u1, u4)
15: procedure x4_UB_2424(B: 2DArray) ▷ Calculate overestimation
16: return CB,2,4 x4_secant(x2, x4, u2, u4)
17: procedure x4_UB_3434(B: 2DArray) ▷ Calculate overestimation
18: return CB,3,4 x4_secant(x3, x4, u3, u4)
19: procedure x4_LB_3434(B: 2DArray) ▷ Calculate underestimation
20: return CB,3,4 x4_tangent(x3, x4, u3, u4)
21: procedure calc_01234_LB(B: 2DArray) ▷ Underestimation of path length 0-1-2-3-4
22: return x4_0101(B) + x4_1212(B) + x4_2323(B) + x4_LB_3434(B)
23: procedure calc_01324_UB(B: 2DArray) ▷ Overestimation of path length 0-1-3-2-4
24: return x4_0101(B) + x4_1313(B) + x4_2323(B) + x4_UB_2424(B)
25: procedure calc_02134_UB(B: 2DArray) ▷ Overestimation of path length 0-2-1-3-4
26: return x4_0202(B) + x4_1212(B) + x4_1313(B) + x4_UB_3434(B)
27: procedure calc_02314_UB(B: 2DArray) ▷ Overestimation of path length 0-2-3-1-4
28: return x4_0202(B) + x4_2323(B) + x4_1313(B) + x4_UB_1414(B)
29: procedure calc_03124_UB(B: 2DArray) ▷ Overestimation of path length 0-3-1-2-4
30: return x4_0303(B) + x4_1313(B) + x4_1212(B) + x4_UB_2424(B)
31: procedure calc_03214_UB(B: 2DArray) ▷ Overestimation of path length 0-3-2-1-4
32: return x4_0303(B) + x4_2323(B) + x4_1212(B) + x4_UB_1414(B)
56

constant β2.
Algorithm 4 Calculation of contribution ϵ(B, πB), given hyperrectangle B and permutation index πB
1: procedure contribution(B: 2DArray, label: int) ▷ Calculate the contribution of B, given πB
2: if label = 0 then
3: return 0
4: identity ← iv.mpf(calc_01234_LB(B)) ▷ Apply interval arithmetic
5: if label = 1 then
6: alternative ← iv.mpf(calc_01324_UB(B))
7: else if label = 2 then
8: alternative ← iv.mpf(calc_02134_UB(B))
9: else if label = 3 then
10: alternative ← iv.mpf(calc_02314_UB(B))
11: else if label = 4 then
12: alternative ← iv.mpf(calc_03124_UB(B))
13: else if label = 5 then
14: alternative ← iv.mpf(calc_03214_UB(B))
15: return max(0, identity − alternative)
16: ▷ Contribution becomes zero if identity permutation outperforms all alternatives
Algorithm 4 consolidates the methodologies discussed earlier to evaluate the contribution ϵ(B, πB) for
a specified hyperrectangle B and permutation index πB. It juxtaposes the path lengths generated by the
identity permutation with that from a designated alternative permutation, assigning a zero contribution
when the identity permutation is optimal. Utilizing interval arithmetic [62] (iv.mpf) on the computations
derived from FriCAS [44] and SageMath [83], the algorithm ensures the precision of bound improvements
for β2. These improvements are quantified by the lower bound within the interval arithmetic. The subsequent section delves into the computational techniques and outcomes, showcasing the systematic process
employed to refine the upper bound of β2, thereby encapsulating the algorithm’s strategic value to handle
complex combinatorial optimization challenges.
57

Chapter 6
Computational Experiments
6.1 A new upper bound of the travelling salesman constant
The refinement of the upper bound of β2 predominantly aims at constraining the right-hand side of (5.2).
Mathematical proofs presented in Chapter 5 introduce a methodology for bounding the right-hand side
of (5.2) through the accumulation of positive contribution ϵ (5.3), derived from an extensive collection of
disjoint hyperrectangles. This chapter will elaborate on the construction of these non-overlapping hyperrectangles, alongside determining the most cost-effective permutation π
∗
B
and calculating the contribution
ϵ(B, π∗
B
) for any selected hyperrectangle B.
6.1.1 Selection of hyperparameters
Initiating this discussion, it is pivotal to decide on specific values for h and k. Echoing the seminal work of
[19] and the insights from [94], the value h =
√
3 is adopted due to its proven effectiveness in delivering
an optimal default upper bound ≈ 0.92116027. Opting for k = 4 considers the significant computational
demand associated with higher k values. Nonetheless, exploring alternative settings for the hyperparameters h and k is considered advantageous, as it holds the potential to further refine the guaranteed upper
bound of β2.
58

h
2 Upper Bound of (3.1) 95% CI of (5.1) with 5-tuples 95% CI of (5.1) with 6-tuples
2 0.94892720 0.934055 ± 5.3 × 10−5 0.931268 ± 3.0 × 10−5
3 0.92116027 0.890195 ± 4.1 × 10−5 0.884025 ± 3.7 × 10−5
4 0.93471669 0.884960 ± 3.2 × 10−5 0.874237 ± 2.5 × 10−5
5 0.96432637 0.894581 ± 3.0 × 10−5 0.878494 ± 3.3 × 10−5
6 1.00102758 0.910974 ± 2.8 × 10−5 0.888951 ± 2.9 × 10−5
7 1.04094797 0.930837 ± 3.4 × 10−5 0.902545 ± 2.6 × 10−5
8 1.08221001 0.952621 ± 2.7 × 10−5 0.917852 ± 4.4 × 10−5
Table 6.1: For zigzag configurations involving 5-tuples and 6-tuples, a prioritized task is to determine the
improvement limit for the upper bound of β2 across different settings of hyperparameter h. The righthand side of (5.1) is simulated for each h and k pair, with the 95% confidence interval obtained through
bootstrapping.
Recall that Section 3.1 details a process for simulating (k+1)-tuples and evaluating path lengths across
all (k − 1)! permutations, for any given positive integer k and positive number h. An investigation into
varying h values that might offer a more refined upper bound for β2, beyond the conventional choice of
h =
√
3, has been conducted for k = 4, 5 and h
2 values ranging from 2 to 8, using 30 million samples of
(k + 1)-tuples each.
For each pair of h and k, the evaluation of the right hand side of (5.1), with the 95% confidence interval determined through bootstrapping, is documented in Table 6.1. Notably, for k = 4 and h =
√
3,
the simulated upper bound is approximately 0.8902, corroborating the simulation findings of [94]. This
result indicates a practical limit to the potential for improving the upper bound of β2 with the default
hyperparameter setting k = 4 and h =
√
3.
According to Table 6.1, the scenario with k = 5 (involving 6-tuples) consistently yields lower (and
therefore better) values for the quantity of interest described in (5.1) compared to when k = 4. Additionally,
opting for a suitable h value shows promise for achieving a shorter expected travelling salesman tour than
applying the default choice of h =
√
3. Specifically, for k = 4, 5, selecting h = 2 appears to be the
most advantageous theoretically, assuming an immaculate mechanism for targeting hyperrectangles B and
bountiful computational resources. With k = 5 and h = 2, the methodology presented in Chapter 5 could
at its best enable an upper bound of β2 of approximately 0.87424, taking into account all 24 permutations
and selecting B in a perfect manner.
59

Figure 6.1: The computational experiments for h
2 values of 3, 3.5, 4, 4.5, 5, and 5.5, across a range of k from
3 to 9, showcase the potential improvement limits for the upper bound of β2. This investigation reveals that
the optimal choice of h tends to increase with larger k, suggesting the inadequacy of the default h =
√
3
when computational capabilities permit the exploration of higher k values.
The requirement for more computational resources grows significantly with an increase in k, attributed
to the superexponential surge in permutation possibilities. Nonetheless, this increase in computational
demand is compensated by the increased opportunity for refining the upper bound, as illustrated in Table
6.1 and Figure 6.1.
Figure 6.1 elucidates how higher k values can lead to improved simulation outcomes and demonstrates
the impact of different h values on these results. The best improvement limit (5.1) for the upper bound of
β2 is an estimated 0.852, attainable at k = 9 and h
2 = 5. A notable trend is the inclination for the optimal
h to increase with k, suggesting that a larger h could result in a more substantial cumulative contribution
from the selected hyperrectangles. This aspect will be expanded upon in the following section.
60

Finally, as highlighted in Section 5.3, the adoption of the partition hyperparameter M = 100 for the
one-dimensional numerical integration of the last summand in (5.4) represents a strategic compromise
aimed at increasing computational efficiency while preserving the solution quality.
Considering the current computational constraints, a focused investigation on 5-tuples is deemed most
practical for the selection of hyperrectangles. Future studies could extend this methodology to examine the
asymptotic behaviors and potential convergence of the upper bound of β2, or more broadly βd, as k continues to increase. Such research endeavors hold the promise of unveiling deeper insights into the intricate
dynamics of combinatorial optimization problems in the context of evolving computational capabilities.
6.1.2 Construction of hyperrectangles with a decision tree
This section delves into the procedure for identifying hyperrectangles within the multidimensional space
D, encapsulating the strategies employed to navigate and partition this hyperspace for the purposes of upper bound refinement. To effectively refine the upper bound of β2, it is paramount to delineate a substantial
array of disjoint hyperrectangles B, along with their most efficacious permutation π
∗
B
. The definitions of
B, π
∗
B
, and the contribution metric ϵ(B, π∗
B
) are comprehensively outlined in Section 5.2.
For a specific hyperrectangle B ⊂ D, π
∗
B
is designated as the permutation that yields the maximum
improvement to the upper bound of β2. It is pivotal to acknowledge, however, that permutations not
achieving this optimal status can still offer positive contributions to the upper bound of β2, provided
that they surpass the baseline performance established by the identity permutation πId. Consequently,
pinpointing the absolute optimal permutation π
∗
B
for each B is not always necessary. The emphasis should
indeed be placed on formulating an effective partition strategy for the hyperspace D, alongside a method
for transforming each partitioned segment into a disjoint hyperrectangle. For every hyperrectangle, there
lies the potential for alternative permutations to outperform the identity permutation, and the crucial task
is to accurately identify and categorize these instances.
61

Algorithm 5 Implementation of node and tree classes
1: procedure generate_data( ) ▷ Generate 2DArray with feature and label columns
2: generate data ▷ Features: xi
, ui
; label: permutation index
3: procedure calc_cross_entropy(data: 2DArray)
4: return cross_entropy ▷ Calculate cross entropy
5: procedure calc_info_gain(feature: int, threshold: float, data: 2DArray)
6: call calc_cross_entropy
7: return info_gain ▷ Calculate information gain
8: procedure calc_threshold(data: 2DArray)
9: call calc_info_gain
10: return (feature, threshold) ▷ Determine the optimal split
11: class Node
12: attribute (lef t: node, right: node, label: int, feature: int, threshold: float, ancestor: List[tuple])
13: procedure insert(lef tchild: node, rightchild: node)
14: lef t ← lef tchild; right ← rightchild ▷ Link parent node to children nodes
15:
16: class Tree
17: attribute (ancestor: List[tuple])
18: procedure create_decision_tree(data: 2DArray, leaf_nodes: List[node])
19: if data is homogeneously labeled by label then
20: leaf ← Node(None, None, label, None, None, ancestor)
21: leaf_nodes.append(leaf)
22: return leaf
23: else
24: feature, threshold ← calc_threshold(data)
25: parent ← Node(None, None, None, feature, threshold, ancestor)
26: ancestor.append((feature, threshold, “L”))
27: lef tchild ← create_decision_tree(data[data[:, feature] < threshold], leaf_nodes)
28: ancestor.pop()
29: ancestor.append((feature, threshold, “R”))
30: rightchild ← create_decision_tree(data[data[:, feature] ≥ threshold], leaf_nodes)
31: ancestor.pop()
32: parent.insert(lef tchild, rightchild)
33: return parent ▷ Construct decision tree using recursion
34:
Algorithm 5 efficiently constructs a decision tree with two primary classes: Node and Tree. The
Node class holds crucial attributes like pointers to child nodes, labels, and the criteria for splitting at
the current node. The Tree class, through its recursive method create_decision_tree, orchestrates the
construction of the tree. It starts with the entire dataset and iteratively divides it based on predetermined split criteria, thereby generating subtrees for each branch. This methodology allows for the
dynamic construction of the tree, with each node and its corresponding subtree tailored to the subset
of data it represents.
62

Figure 6.2: Illustration of a decision tree in an experiment with a small sample of N = 1000 hyperpoints.
Each node within the tree represents a subset of rows in the original dataset, with non-leaf nodes carrying
decisions based on specific features and threshold values. Leaf nodes, which can emerge at any depth beyond the root, are defined by a homogeneity in their labels (permutation indices). Additionally, every node
is associated with an ancestor attribute, a list that traces its lineage back to the root, with its length corresponding to the node’s depth in the tree. The root node’s list of ancestors is notably empty, highlighting
its origin status in the tree structure.
63

Utilizing a decision tree for binary space partitioning, in conjunction with a Monte Carlo simulation,
emerges as an effective strategy for achieving this goal, as illustrated in Figure 6.2. This methodology
enables the systematic identification and categorization of a multitude of non-overlapping hyperrectangles
B ⊂ D. The initial phase of this methodology involves generating a considerable volume, denoted by N,
of hyperpoints p1, p2, . . . , pN ∈ D, where
D = {(x1, . . . , xk, u0, . . . , uk) : 0 ≤ x1 ≤ · · · ≤ xk , 0 ≤ ui ≤ 1} ⊂ R
2k+1
.
In this process, ui values are sampled independently and uniformly, while xi values are sampled in accordance with the distribution specified in Lemma 9. Following the completion of this procedure, N hyperpoints, each representing a (k+ 1)-tuple in the unit square, are independently generated. Each hyperpoint
is then assigned a label corresponding to one of the (k − 1)! categorical outcomes, indicating its true
optimal permutation. This demonstrates the data generation process within Algorithm 5.
With the array of N hyperpoints, each characterized by (2k + 1) dimensions, along with their designated categorical outcomes, the training of a decision tree becomes a viable approach. This decision tree
is tasked with partitioning the hyperspace R
2k+1 into segments through the implementation of separating
hyperplanes. By employing cross entropy to ascertain the splitting thresholds, the training of the tree
results in a segmentation of R
2k+1 according to the delineation provided by the tree’s leaves.
It is noteworthy that the branching mechanism persists recursively until the subsets at the leaf nodes
exhibit uniform categorical outcomes, i.e., permutation indices. At the end of this branching process, a
series of leaf nodes emerges, categorized either under the identity permutation or an alternative permutation. The next phase of the algorithm entails isolating those subsets at the leaf nodes associated with
alternative permutations, and then refining them into precisely defined hyperrectangles.
64

Figure 6.3: Given k = 4, an example hyperrectangle B ⊂ D ⊂ R
9
is depicted as a combination of
one line segment and four rectangular regions. Additionally, the representation of a 9-dimensional hyperpoint as five different points on the 2-dimensional coordinate plane bridges the gap between abstract
high-dimensional data and its tangible graphical interpretation. For this specific hyperpoint, the visual
contrast between the identity permutation (0,1,2,3,4) and the optimal alternative permutation (0,1,3,2,4) is
highlighted. For a specified hyperrectangle B, the goal is to discern an alternative permutation that, on
average, outperforms the identity permutation across the infinite set of hyperpoints within this hyperrectangle.
The decision tree’s architecture intrinsically ensures that the subsets represented by its leaf nodes
are non-overlapping, facilitating a systematic division of the search space into segments bounded by clear
linear boundaries, as dictated by the decision criteria at each node of the tree. This characteristic makes the
leaf nodes particularly conducive to forming a collection of hyperrectangles. These valid hyperrectangles
are then employed to contribute towards the refinement of the upper bound of β2.
65

Revisiting the concept of the hyperrectangle, it is noted that for fixed positive integer k, any selected
hyperrectangle B ⊂ D ⊂ R
2k+1 is distinctly characterized by 4k + 2 boundary parameters. It can be
defined as
B([s1, t1], . . . , [sk, tk], [a0, b0], . . . , [ak, bk]) = {(x1, . . . , xk, u0, . . . , uk) : ai ≤ ui ≤ bi
, sj ≤ xj ≤ tj} ,
where the parameters ai
, bi
, sj , tj are constrained to fulfill the conditions
0 ≤ s1 ≤ t1 ≤ . . . ≤ sk−1 ≤ tk−1 ≤ sk ≤ tk, 0 ≤ ai ≤ bi ≤ 1 .
Figure 6.3 serves as a detailed visual guide to understanding a 9-dimensional hyperrectangle for k = 4.
This hyperrectangle B is visually translated onto a 2-dimensional coordinate plane, represented as a mix
of one line segment and four rectangular areas. It is precisely defined by
B([0, 0.4], [0.5, 1], [1, 1.6], [1.7, 2], [0, 0.18], [0.9, 1], [0, 0.15], [0.88, 1], [0, 0.08]) ⊂ D ⊂ R
9
.
For each hyperrectangle B depicted on the coordinate plane, the enhanced capability of its hyperpoints
(i.e., (k + 1)-tuples) to identify favorable zigzags over standard left-to-right paths would increase the
likelihood of this hyperrectangle to accompany an alternative permutation π
∗
B
that yields a significant
contribution ϵ(B, π∗
B
) (5.3).
Dividing a larger hyperrectangle into smaller subsets and allocating the optimal permutation to each
ensures that the cumulative contribution is at least as substantial as in the initial configuration. Hence,
66

further subdivision is always encouraged, considering the minimization problem related to further segmentation as a relaxation of the original minimization problem tied to the larger hyperrectangle. Therefore, augmenting the number of sampled hyperpoints prior to beginning the construction of the decision
tree is anticipated to facilitate a further refinement of the upper bound.
Last but not the least, in the process of partitioning the hyperspace R
2k+1 into hyperrectangles, it
is often observed that the resulting hyperrectangles can extend beyond the boundary of D, rendering
them unsuitable for assessing contribution to the upper bound of β2. Note that D can be conceptualized
as a simplex-hypercube, combining a k-dimensional simplex defined by the xi variables with a (k + 1)-
dimensional hypercube defined by the ui variables.
For any hyperrectangle B
′
⊂ R
2k+1 located at the leaf node of the decision tree, if there exists a nonempty difference set B
′
\ D, it is imperative to identify a valid hyperrectangle B such that B ⊂ B
′
∩ D.
The natural approach involves selecting the largest-volume hyperrectangle B that meets the boundary
conditions of D via convex optimization.
Considering computational efficiency, various alternative methods exist to locate valid hyperrectangles
B ⊂ D. This section delves into the strategies for pruning the hyperrectangular sets B
′
located at the leaf
nodes of the decision tree, with the objective of generating valid hyperrectangular subsets B within the
simplex-hypercube domain D. This pruning process is critical for generating hyperrectangles that are not
only valid within the boundaries of D but also capable of facilitating a substantial improvement to the
upper bound of the Euclidean travelling salesman constant, β2.
Here, an illustrative example is provided in Figure 6.2, focusing on the decision tree’s leaf node marked
“LRRRL”. By tracing the ancestry of this hyperrectangular set, the constraints defining B
′
are specified as
x4 ≤ 3.695, x3 ≥ 1.749, x1 ≥ 1.013, x2 ≥ 2.679, u1 ≤ 0.382,
67

which reveal the initial boundaries of B
′
within R
9
as
B
′
([1.013, ∞], [2.679, ∞], [1.749, ∞], [0, 3.695], [0, 1], [0, 0.382], [0, 1], [0, 1], [0, 1]) ̸⊂ D.
Notice the absence of explicit constraints for u0, u2, u3, u4. By imposing the necessary ordering conditions
s1 ≤ . . . ≤ sk and t1 ≤ . . . ≤ tk on the simplex coordinates, B
′
is further pruned to
B
′
([1.013, 3.695], [2.679, 3.695], [2.679, 3.695], [2.679, 3.695], [0, 1], [0, 0.382], [0, 1], [0, 1], [0, 1]) ̸⊂ D .
Adjustments to the bounds of B
′
are carried out to ensure that the resulting hyperrectangle B fits within
the simplex-hypercube D, and it brings potential for a meaningful contribution ϵ(B, πB) towards improving the upper bound of β2. The algorithm achieves this by refining B
′
to
B([1.013, 2.679], [2.679, 3.018], [3.018, 3.356], [3.356, 3.695], [0, 1], [0, 0.382], [0, 1], [0, 1], [0, 1]) ⊂ D ,
where B is meticulously crafted to ensure maximum volume within the boundaries of D ⊂ R
9
. It’s critical
to note, however, that this method does not aim to identify the global optimal ϵ(B, πB) for B
′
emanating
from the leaf nodes of the decision tree. Rather, the goal is to attain a substantial, albeit potentially suboptimal, contribution ϵ(B, πB). This reflects a deliberate strategic compromise, balancing the pursuit of
global optimal solutions with computational practicality. Such a balance is sometimes essential within the
complex landscape of combinatorial optimization problems.
As the number of sampled hyperpoints N increases, the occurrence of scenarios lacking explicit constraints for certain dimensions in the hyperrectangular set B
′
becomes less frequent. This reduction in
ambiguity directly correlates with a diminished necessity for extensive pruning from B
′
to B, streamlining
the process of refining hyperrectangles for analysis.
68

The entire process, from receiving a leaf node defined by the hyperrectangular set B
′
, through the
pruning phase to derive a valid hyperrectangle B, to the final calculation of contribution ϵ(B, πB), is encapsulated in Algorithm 6.
Algorithm 6 Calculation of contribution for a leaf node in the decision tree
1: procedure necessary_condition(B
′
: 2DArray) ▷ Refine B
′
based on simplex constraints
2: B
′
[1, 0], B′
[2, 0] ← max(B
′
[0, 0], B′
[1, 0]), max(B
′
[0, 0], B′
[1, 0], B′
[2, 0])
3: B
′
[3, 0] ← max(B
′
[0, 0], B′
[1, 0], B′
[2, 0], B′
[3, 0])
4: B
′
[1, 1], B′
[2, 1] ← min(B
′
[1, 1], B′
[2, 1], B′
[3, 1]), min(B
′
[2, 1], B′
[3, 1])
5: B
′
[0, 1] ← min(B
′
[0, 1], B′
[1, 1], B′
[2, 1], B′
[3, 1])
6: return B
′
7: procedure prune(B
′
: 2DArray) ▷ Ensure B fits within D
8: apply pruning rules to adjust B
′
into a valid hyperrectangle B, maximizing its volume
9: return B
10: procedure leaf2box(leaf: node) ▷ Construct hyperrectangle from a leaf node
11: B
′ ← 2DArray([0, ∞], [0, ∞], [0, ∞], [0, ∞], [0, 1], [0, 1], [0, 1], [0, 1], [0, 1])
12: for constraint in leaf.ancestor do
13: apply constraint to refine B
′
14: B
′ ← necessary_condition(B
′
) ▷ Impose necessary conditions
15: B ← prune(B
′
) ▷ Maximize volume of B ⊂ D
16: return B
17: procedure check_monotone(B: 2DArray) ▷ Confirm the validity of B
18: return boolean
19: procedure leaf_contribution(leaf: node) ▷ Derive ϵ(B, πB) from a leaf node
20: B ← leaf2box(leaf)
21: if check_monotone(B) then
22: return contribution(B, leaf.label) ▷ See Algorithm 4
6.1.3 Computational results
Building upon the methodologies outlined in the previous section for constructing hyperrectangles, this
section transitions into the presentation of computational results. The experiments underpinning these
results were conducted on an Intel Xeon E5-2640 v3, 2.60 GHz CPU with 59 GB memory per node, utilizing
the computational resources provided by the USC Center for Advanced Research Computing (CARC).
69

Figure 6.4: Comparative analysis of the upper bound of β2 against sample size N, for hyperparameter
values h
2 = 3, 3.25, 3.5. The solid lines depict the new guaranteed upper bounds of β2, which decrease
with increasing N and significantly outperform Steinerberger’s benchmark. The bars indicate computation
time, revealing an expected increase with larger sample sizes.
This section keeps the selection of k = 4, maintaining the balance between solution quality and computational efficiency, whereas our framework is applicable for any given k. The best guaranteed outcomes
for the upper bound of the Euclidean travelling salesman constant β2 were systematically explored over
a range of sample sizes N and hyperparameter settings h. For every combination of (h, N), 10 separate
decision trees were constructed to ensure the robustness of the results, with h
2 ∈ {3, 3.25, 3.5, 3.75, 4},
and N ∈ {104
, 104.5
, 105
, 105.5
, 106
, 106.5
, 107
, 107.5
, 108}. The findings from Figure 6.4 underscore that
larger sample size N consistently amplifies the aggregate contribution across different hyperparameter
settings, which is a testament to the underlying relaxation properties of the binary hyperspace partitioning strategy employed in this research.
When examining the influence of the hyperparameters h on the aggregate contribution, another interesting trend emerges in Figure 6.4: smaller h values lead to better results at lower sample sizes (N), while
larger h values begin to outperform them as N increases. This trend, consistently observed across the
70

N Upper bound of β2
104 0.92104
104.5 0.91919
105 0.91786
105.5 0.91652
106 0.91457
106.5 0.91258
107 0.91031
107.5 0.90818
108 0.90604
108.5 0.90403
3.6 · 108 0.90380
Table 6.2: The best upper bounds of β2 obtained by our method, for increasing values of N. The overall
best upper bound β2 < 0.9038 is highlighted.
three curves for h
2 values of 3, 3.25, and 3.5, resonates with the earlier findings depicted in Figure 6.1. Additionally, the computation time, while increasing with N, exhibits only marginal changes across different
h values when plotted on a logarithmic scale. This characteristic is advantageous for computational experiments as it allows for flexibility in hyperparameter tuning without incurring significant computational
inefficiencies.
To mitigate any potential floating-point errors, the dissertation utilized the mpmath Python package,
known for its variable-precision interval arithmetic. This choice aligns with the methodological rigor, as
demonstrated in Algorithm 4, ensuring the precision of the computational analysis. For the purpose of
transparency and to promote the reproducibility of the research, the source code and the comprehensive
computational results have been made publicly accessible via the referenced GitHub repository [101].
Due to the substantial memory requirements of managing exceedingly large decision trees and the oneweek computation time limit imposed by CARC, the current capacity to process hyperpoints was capped
around N = 3.6 · 108
. Within this constraint, Figure 6.4 indicates that the hyperparameter setting h
2 =
3.25 is expected to yield the most effective results for the new upper bounds of β2 at the given magnitude
of sample size. Table 6.2 succinctly captures the improvements in the best guaranteed upper bounds of β2
as N escalates. A single computational run with the hyperparameter h
2 = 3.25 and the sample size of
N = 3.6·108 hyperpoints was performed to test the upper limits of the available computational resources,
71

producing the most refined improvement result of β2 < 0.9038 and thus validating the claim of Theorem
2.
6.2 A lossless splitting rule for split-delivery routing problems
Recall that an a priori splitting rule, elaborated in Chapter 4 of this dissertation, enables the reduction of
the split-delivery vehicle routing problem (SDVRP) into the capacitated vehicle routing problem (CVRP)
through the simple division of customer nodes into several nodes with reduced demands. This process
increases the problem size yet maintains the lossless property, allowing the utilization of existing CVRP
solvers for the SDVRP. This dissertation section undertakes computational experiments on the benchmarks of the SDVRP and its time window variant (SDVRPTW), leveraging the enumeration schemes that
optimally solve the minimum size coalescing partition problem as described in Chapter 4.
Compared to earlier solvers [50, 51, 78], the LKH-3 solver [59, 61] has exhibited effectiveness in addressing constrained vehicle routing problems including the CVRP and CVRPTW, making it the solver
of choice for the non-split CVRP instances generated by the lossless a priori splitting rule. The results
obtained are competitive with the current state-of-the-art for both the SDVRP and SDVRPTW, achieved
without significant adjustments to the metaheuristics. Notably, when comparing the computational outcomes in this section with those derived from other a priori splitting rules, such as the strategy developed
by [29] as mentioned in Section 3.2, the method in this dissertation exhibits a marked advantage.
The performance evaluation of our proposed method spans two sets of SDVRP benchmarks and one
set of SDVRPTW benchmarks. The first set, presented by [20], includes 25 instances. This set features
the exact Christofides and Eilon instances from TSPLIB [82] and additional instances with the same coordinates as eil51, eil76, eil101, and a fixed vehicle capacity of Q = 160. Customer demands across these
72

instances are randomly generated according to various ranges, proportionally scaled to the vehicle’s capacity: D1: [⌈0.01Q⌉, ⌊0.1Q⌋]; D2: [⌈0.1Q⌉, ⌊0.3Q⌋]; D3: [⌈0.1Q⌉, ⌊0.5Q⌋]; D4: [⌈0.1Q⌉, ⌊0.9Q⌋]; D5:
[⌈0.3Q⌉, ⌊0.7Q⌋]; D6: [⌈0.7Q⌉, ⌊0.9Q⌋].
The second set of benchmarks, devised by [13], comprises 42 instances, each applying varied demand
generation strategies across six different customer maps, distinguished by the diversity in customer count
and location. The first and second sets of SDVRP benchmark instances have been extensively examined
in a wide range of researches [2, 11, 20, 29, 54, 76, 88], with [54, 88] achieving most of the best known
solutions.
The third set introduces SDVRPTW benchmark instances, initially proposed by [90] and categorized
into groups C, RC, and R. The computational experiments detailed in this section focus on the C and RC
groups, featuring 58 instances notable for lacking feasible non-split solutions. The best known SDVRPTW
solutions for these instances have been recognized and documented by [12, 38, 69].
Preliminary computational experiments were conducted on the classical benchmark Set 1, as introduced by [20]. These experiments aim to evaluate the impact of constraining the sizes of non-split CVRP
instances generated by the a priorisplitting rule on both the computational time and the quality of solutions
for the SDVRP. Furthermore, the effectiveness of the splitting rule in tackling more complex split-delivery
routing challenges was further confirmed through a comparative analysis with the leading solutions for
the SDVRP problems in Set 2 and the SDVRPTW problems in Set 3.
The LKH-3 solver was executed 10 times for each problem instance, employing non-default hyperparameter settings labeled “SPECIAL” and modifying the maximum number of trials in each run through the
hyperparameter “MAX_TRIALS”. Details can be found in the description of hyperparameters for LKH-3
[61]. By standard configuration, LKH-3 addresses the split-delivery vehicle routing problem with a limited fleet (SDVRP-LF) rather than the variant with an unlimited fleet (SDVRP-UF), setting the fleet size
to the smallest feasible number of vehicles, denoted as kmin = ⌈
Pn
i=1 di/Q⌉. Consequently, for the first
73

and second instance sets, the computational outcomes are benchmarked against the best known SDVRPLF solutions. For the third instance set, concerning split-delivery routing problems with time windows,
the derived SDVRPTW-LF solutions are evaluated in comparison to the best known SDVRPTW solutions
documented by [12, 38, 69].
All experiments were conducted using Intel Xeon E5-2640 v3, 2.60 GHz CPUs at the USC Center for
Advanced Research Computing (CARC), equipped with 59 GB of memory per node. The computational
time was limited to three hours for instances from [20], and extended to six hours for instances from [13,
90].
6.2.1 Limiting the non-split problem size
In the tables presented from now on, several key terms are defined for clarity. Instance refers to the name
of the problem instance. k represents the minimum feasible number of vehicles, denoted before as kmin.
BKS stands for the best known solution for the SDVRP (or SDVRPTW) benchmark instance. Best sol.
indicates the best solution identified by the a priori splitting rule for the SDVRP (or SDVRPTW) across
all executions. Gap (%) is the percentage difference between the best solution and the BKS, calculated as
100 × (Best sol. − BKS)/BKS. Time (s) denotes the time in seconds to achieve the best solution. Solutions
that equal or surpass the BKS will be highlighted in bold.
Building on the observations from Chapter 4, it is postulated that the extent of demand splitting across
trucks is likely not excessively large in practical scenarios. It is generally expected that some ¯k, significantly
lower than k, exists, ensuring that no more than ¯k trucks are required to service any specific customer. This
premise facilitates a significant reduction in the sizes of the generated non-split CVRP instances requiring
solutions.
Meanwhile, the hyperparameter q dictates that only customer demands greater than qQ should be
considered for division using the minimum size coalescing partition strategy. This approach stems from
74

the understanding that segmenting small customer demands into smaller units does not typically enhance
solution quality but tends to adversely affect computational time. When q = 1 (or more precisely when
q ≥ dmax/Q), the scenario simplifies to that of a non-split capacitated vehicle routing problem instance.
Conversely, when q = 0 (or more precisely when q < dmin/Q), the scenario of a fully relaxed split-delivery
vehicle routing problem needs to be addressed.
In the forthcoming Tables 6.3 and 6.4, the columns for Best ¯k and Best q are included to showcase the
effects of ¯k and q on both the solution quality and computation time of the generated non-split CVRP
instances. Specifically, ¯k and q are defined within the domains ¯k ∈ {2, . . . , k} and q ∈ [0, 1], respectively.
A lesser value of ¯k yields a reduced size for the generated non-split problem, and similarly, an increased
q value also diminishes the generated non-split problem size. Overall, this measure demonstrates how
¯k and q can be effectively utilized to reduce the non-split problem size, in comparison to completing the
minimum size coalescing partitioning procedure, while maintaining or enhancing the solutions for SDVRP
benchmark instances.
Lastly, the Full size column in Tables 6.3 and 6.4 refers to the comprehensive magnitude of the generated non-split problem that guarantees a lossless solution to each split-delivery routing problem instance.
Setting ¯k = k and q = 0, it represents the aggregate cardinality of the minimum size partitions µi that
coalesce to all k-partitions of each customer demand di
. Prob size signifies the size of the non-split problem
created with the specified ¯k and q, which facilitate the best split-delivery solution identified by the LKH-3
solver. Ratio (%) denotes the percentage ratio of the problem size that culminates in the best outcome,
relative to the comprehensive non-split problem size, computed as 100 × Prob size/Full size.
Upon reviewing all 25 instances from Set 1, Figures 6.5 and 6.6 depict how different choices of ¯k and q
influence the size of the non-split problem and the gap of the LKH-3 solution from the best known solution, respectively. The heat map in Figure 6.6 reveals that the solution gap tends to be higher in the bottom
right corner (where ¯k is large and q is small) and at the top (where q is large). This outcome is logical, as
75

Figure 6.5: The a priori splitting rule, when applied to Instance Set 1 from [20], results in varying sizes of
non-split capacitated vehicle routing problem instances for different (¯k, q) combinations. The generated
problem size increases with larger ¯k or smaller q. It is important to note that when q ≥ dmax/Q, the
resulting problems effectively disallow splitting, hence such instances are excluded from Figure 6.5.
76

Figure 6.6: The gap between the solution obtained using the a priori splitting rule and the best known
solution, for varying selections of (
¯k, q), analyzed across all instances of Set 1 [20].
77

an excessively large problem size may lead to diminished solution quality due to the constraints of limited
computation time. In contrast, a higher q limits the permissible splitting, producing a capacitated vehicle
routing problem that is more akin to the problem instance which prohibits splitting. Thus, the ideal selections for q appear to be neither excessively small (approaching 0) nor excessively large (approaching
dmax/Q), but rather an intermediate value.
Instance k BKS Best sol. Gap (%) Timea
(s) Prob size Full size Ratio (%) Best k¯ Best q
eil22 4 375b
375 0.000 0.01 33 455 7.25 2 0.40
eil23 3 569c
569 0.000 0.01 42 287 14.63 3 0.90
eil30 3 510b
510 0.000 0.02 64 381 16.80 3 0.30
eil33 4 835c
835 0.000 0.07 44 654 6.73 2 0.45
eil51 5 521b
521 0.000 0.43 59 417 14.15 3 0.25
eilA76 10 818c
818 0.000 308.30 465 989 47.02 4 0.10
eilB76 14 1002c
1003 0.100 30.83 187 1137 16.45 2 0.20
eilC76 8 733c
733 0.000 86.44 322 887 36.30 2 0.05
eilD76 7 681c
682 0.147 0.11 90 814 11.06 3 0.15
eilA101 8 815c
815 0.000 0.59 137 1005 13.63 3 0.15
eilB101 14 1061c
1062 0.094 269.18 427 1246 34.27 3 0.10
S51D1 3 458c
458 0.000 0.23 167 220 75.91 2 0.00
S51D2 9 703c
705 0.284 34.99 381 788 48.35 3 0.10
S51D3 15 942e
945 0.318 10.31 265 1288 20.57 3 0.25
S51D4 27 1551e
1560 0.580 10584.38 842 2381 35.36 6 0.20
S51D5 23 1328c
1334 0.452 528.29 745 2092 35.61 8 0.35
S51D6 41 2153d
2183 1.393 2635.05 1035 3804 27.21 6 0.70
S76D1 4 592c
592 0.000 17.67 292 400 73.00 4 0.05
S76D2 15 1081c
1086 0.463 848.87 413 1646 25.09 2 0.05
S76D3 23 1419c
1427 0.564 1595.84 611 2446 24.98 4 0.20
S76D4 37 2071c
2085 0.676 1726.02 624 3852 16.20 3 0.20
S101D1 5 716c
716 0.000 40.42 326 567 57.50 2 0.00
S101D2 20 1360e
1375 1.103 473.32 542 2468 21.96 2 0.10
S101D3 31 1858e
1874 0.861 2001.65 1091 3772 28.92 7 0.25
S101D5 48 2767e
2819 1.879 5832.25 563 6168 9.13 3 0.45
Average 0.357 1081.01 28.72 3.4 0.24
a
Intel Xeon E5-2640 v3, 2.6 GHz
b Belenguer et al. (2000) [20]
c
Silva et al. (2015) [88]
d Ozbaygin et al. (2018) [76]
e He and Hao (2023) [54]
Table 6.3: LKH-3 results for SDVRP with rounded costs and limited fleet on the instances of Belenguer et
al. [20]
Table 6.3 focuses on the 25 Set 1 instances introduced by [20], showcasing the results of applying the
lossless a priori splitting rule, in conjunction with the LKH-3 solver. This analysis considers all practical q
values, initiating at 0 and increasing in increments of 0.05, up to 1. For instances where different q values
result in the same non-split problem, the splitting rule utilizes the highest q value. The range of ¯k values
78

analyzed is {2, 3, . . . , 8}, based on the rationale that allocating more than this number of vehicles to service
the same customer is unlikely in practical scenarios.
Findings on the selection of ¯k and q are also presented in Table 6.3. The evaluation of different (¯k, q)
combinations prioritizes solution quality, followed by computational runtime. For each instance, Table 6.3
details the (¯k, q) pair leading to the best solution as determined by the LKH-3 solver, along with the full
and reduced non-split problem sizes, as well as the reduction ratio. It is observed that in most scenarios,
¯k = 2 and ¯k = 3 yield the ideal outcomes, surpassing those obtained with larger values of ¯k, which
correspond to increased non-split problem sizes. This observation aligns with the rationale that a small
number of vehicles should suffice for serving each customer. Moreover, in most instances, a nonzero
value of q is chosen to limit the number of demand nodes that require splitting, thereby diminishing the
overall problem size. The average selections for ¯k and q are found to be 3.4 and 0.24, respectively. The
average reduction in problem size across the 25 instances is 28.72%, signifying that the problem size can
be reduced to 28.72% relative to the solution to the minimum size coalescing partition (MSCP), without
adversely affecting the solution quality of the SDVRP.
The solutions reported in Table 6.3 demonstrate a match with 11 of the 25 best known solutions for the
SDVRP, achieved through the application of the splitting rule alone, without any advanced adjustments
to the LKH-3 metaheuristics. The average gap between these solutions and the BKS is a modest 0.357%.
Note that the focus of this dissertation is not to modify the metaheuristics but to demonstrate the efficacy
of the a priori splitting rule. The LKH-3 solver, not being able to divide non-split instances further into
separate subtasks for multi-core processor utilization, is optimized for single-thread execution. Consequently, all computational experiments were conducted using a single thread. Despite these constraints,
computational experiments in this section have achieved solution quality comparable to that reported in
[54, 88] within a feasible computation time.
79

Comprehensive computational experiments were also conducted on the 42 instances of Set 2, as introduced by [13], with findings summarized in Table 6.4. It is noted that the minimum coalescing partition sizes for this set are considerably larger on average compared to those in Set 1. Considering the
excessive sizes of non-split problems, evaluations were limited to configurations with ¯k ∈ {2, 3} and
q ∈ {0, 0.05, 0.2, 0.3, 0.5, 0.8}. The average gap from the BKS identified by [54, 88] is 1.122%. Note that
the contributions to metaheuristic development for split-delivery routing problems by [54, 88] are widely
recognized, whereas this dissertation underscores the remarkable efficacy of enumeration schemes in addressing the SDVRP, facilitated by a priori splitting rules. Achieving solution quality on par with the
state-of-the-art is deemed satisfactory.
The result of experiments further indicates a slight preference for ¯k = 2 over ¯k = 3. The average
reduction in problem size for the 42 instances in Set 2 is 17.17%, demonstrating that by permitting two or
three trucks to serve each customer, the problem size can be effectively reduced to 17.17% of the minimum
coalescing partition size.
6.2.2 Adding time windows
This section of the dissertation delves into the split-delivery vehicle routing problem with time windows
(SDVRPTW), as delineated in Definition 5. The investigation employs the a priori splitting rule on a collection of classic instances introduced by [90], specifically focusing on 58 instances from groups C and RC,
with sizes n = 26, 51 and vehicle capacity Q = 30. Instances with the same prefix have uniform node coordinates and demands, but exhibit substantial variations in time windows at the individual instance level.
Notably, all instances feature at least one node with a demand exceeding the vehicle’s capacity, precluding
any feasible non-split solutions and necessitating splitting in the split-delivery context. The best known
solutions for these instances were achieved by [12, 38, 69].
80

In Table 6.5, the column titled No. split reflects the number of splits occurring within the chosen tour
generated by LKH-3. For the 58 instances of Set 3, the average gap between the outcomes produced by the
a priori splitting rule in conjunction with the LKH-3 algorithm and the BKS is a negligible 0.084%, with
the methodology of this dissertation successfully matching 32 of the BKS within Group RC. For the 26
instances categorized under Group C, the average gap slightly increases to 0.112%.
For each instance in Set 3, binary decisions were made regarding ¯k ∈ {2, 3} and q ∈ {0, 0.5}. The
mean value of ¯k calculated is 2.14, indicating a predominant preference for ¯k = 2 over ¯k = 3. The average
value of q stands at 0.44, suggesting a general inclination towards q = 0.5 rather than q = 0. The mean
number of splits in the solutions that were realized for split-delivery is 4.54, confirming that splits do
occur. The average reduction in problem size is 21.72%, showcasing the effective decrease relative to the
minimum coalescing partition size, while still maintaining solution quality.
81

Instance k BKS Best sol. Gap (%) Timea
(s) Prob size Full size Ratio (%) Best k¯ Best q
p01-50 5 524.61b
524.61 0.000 0.02 56 417 13.43 2 0.20
p02-75 10 823.89b
831.79 0.959 1413.20 478 989 48.33 3 0.00
p03-100 8 826.14b
826.14 0.000 0.51 109 1005 10.85 3 0.20
p04-150 12 1023.23c
1029.56 0.619 23.20 167 1817 9.19 3 0.20
p05-199 16 1285.79b
1295.96 0.791 17493.23 216 2814 7.68 3 0.20
p11-120 7 1037.88b
1042.12 0.409 0.42 121 1029 11.76 2 0.20
p01-50D1 3 459.50b
459.50 0.000 0.21 126 228 55.26 2 0.05
p02-75D1 5 617.85b
621.60 0.607 1.28 193 425 45.41 2 0.05
p03-100D1 6 760.00b
760.19 1.006 3161.28 313 803 38.98 2 0.05
p04-150D1 9 921.20c
925.99 0.520 5315.51 466 1440 32.36 2 0.05
p05-199D1 12 1073.57c
1086.14 1.171 10042.39 769 2021 38.05 3 0.05
p11-120D1 8 1042.80c
1042.89 0.009 1819.58 475 1103 43.06 2 0.00
p01-50D2 11 756.71b
768.79 1.596 216.22 278 952 29.20 2 0.00
p02-75D2 16 1109.24c
1117.34 0.730 370.51 265 1628 16.28 2 0.20
p03-100D2 22 1458.46b
1466.20 0.531 2593.03 565 3116 18.13 3 0.20
p04-150D2 32 2016.93c
2051.91 1.734 14050.17 596 5442 10.95 2 0.20
p05-199D2 41 2478.37c
2508.80 1.228 20282.88 1148 7656 14.99 2 0.00
p11-120D2 26 2898.25c
2906.88 0.298 16325.68 681 4029 16.90 3 0.20
p01-50D3 16 1005.75b
1014.91 0.911 5323.18 304 1410 21.56 2 0.00
p02-75D3 24 1502.05b
1517.66 1.039 17693.39 637 2437 26.14 3 0.00
p03-100D3 33 1996.76c
2021.21 1.224 4413.62 638 4663 13.68 3 0.30
p04-150D3 49 2849.59c
2901.45 1.820 16784.15 688 8234 8.36 2 0.30
p05-199D3 63 3469.90c
3525.17 1.593 20120.86 855 11471 7.45 2 0.30
p11-120D3 40 4215.98c
4236.70 0.491 21332.67 770 6103 12.62 3 0.30
p01-50D4 26 1487.18c
1496.13 0.602 17425.97 332 2303 14.42 2 0.00
p02-75D4 40 2301.61b
2319.73 0.787 5149.12 488 4030 12.11 2 0.00
p03-100D4 56 3085.69c
3120.80 1.138 20732.48 721 7815 9.23 2 0.00
p04-150D4 83 4533.82c
4660.13 2.786 21381.17 1109 13814 8.03 3 0.50
p05-199D4 105 5521.61c
5654.87 2.413 21437.06 1368 19020 7.19 3 0.50
p11-120D4 67 6849.73c
6974.97 1.828 15613.14 884 10150 8.71 3 0.50
p01-50D5 26 1481.71b
1487.81 0.412 992.34 506 2371 21.34 3 0.00
p02-75D5 39 2219.11c
2260.89 1.883 21538.35 732 4083 17.93 3 0.00
p03-100D5 53 2986.27c
3025.42 1.311 16483.32 473 7700 6.14 2 0.50
p04-150D5 79 4332.75c
4400.41 1.562 18464.32 1051 13637 7.71 3 0.50
p05-199D5 102 5398.15c
5521.30 2.281 20779.62 890 19384 4.59 2 0.50
p11-120D5 64 6639.95c
6712.05 1.086 10921.86 831 10095 8.23 3 0.50
p01-50D6 41 2155.80c
2172.10 0.756 20473.01 378 3793 9.97 2 0.00
p02-75D6 61 3217.51c
3257.50 1.243 18018.49 526 6493 8.10 2 0.00
p03-100D6 82 4378.33c
4453.35 1.713 19439.10 1201 12116 9.91 3 0.00
p04-150D6 122 6378.28c
6538.92 2.519 20736.82 1201 21278 5.64 2 0.00
p05-199D6 161 8181.44c
8389.22 2.540 20880.95 1593 31536 5.05 2 0.00
p11-120D6 98 10192.00c
10291.36 0.975 21157.62 961 15588 6.16 2 0.00
Average 1.122 11676.24 17.17 2.43 0.16
a
Intel Xeon E5-2640 v3, 2.6 GHz
b
Silva et al. [88]
c He and Hao [54]
Table 6.4: LKH-3 results for SDVRP with limited fleet on the instances of Archetti et al. [13]
82

Instance k BKS Best sol. Gap (%) Timea
(s) Prob size Ratio (%) No. split Best k¯ Best q
C101.25.30 16 821.1b
821.3 0.024 1.77 109 28.02 3 3 0.50
C102.25.30 16 821.1c
821.3 0.024 0.40 81 20.82 3 2 0.50
C103.25.30 16 821.1c
821.3 0.024 0.31 81 20.82 3 2 0.50
C104.25.30 16 821.1c
821.3 0.024 0.18 81 20.82 3 2 0.50
C105.25.30 16 821.1b
821.3 0.024 82.69 81 20.82 3 2 0.50
C106.25.30 16 821.1b
821.3 0.024 3.48 117 30.08 3 2 0.00
C107.25.30 16 821.1c
821.3 0.024 0.70 81 20.82 3 2 0.50
C108.25.30 16 821.1c
821.3 0.024 0.24 81 20.82 3 2 0.50
C109.25.30 16 821.1c
821.3 0.024 0.30 81 20.82 3 2 0.50
C201.25.30 16 909.8b
912.8 0.330 2.01 81 20.82 3 2 0.50
C202.25.30 16 909.8b
912.8 0.330 0.17 81 20.82 3 2 0.50
C203.25.30 16 909.8c
912.8 0.330 0.44 81 20.82 3 2 0.50
C204.25.30 16 909.8c
912.8 0.330 0.20 81 20.82 3 2 0.50
C205.25.30 16 909.8b
912.8 0.330 0.38 81 20.82 3 2 0.50
C206.25.30 16 909.8b
912.8 0.330 0.17 81 20.82 3 2 0.50
C207.25.30 16 909.8c
912.8 0.330 5.65 109 28.02 3 3 0.50
C208.25.30 16 909.8c
912.8 0.330 0.62 81 20.82 3 2 0.50
C101.50.30 29 1599.5d
1599.6 0.006 9.27 147 17.67 4 2 0.50
C102.50.30 29 1599.5d
1599.6 0.006 1.41 147 17.67 4 2 0.50
C103.50.30 29 1598.3c
1598.4 0.006 2.66 147 17.67 4 2 0.50
C104.50.30 29 1598.3c
1598.4 0.006 1.09 147 17.67 4 2 0.50
C105.50.30 29 1599.5d
1599.6 0.006 12.53 198 23.80 4 3 0.50
C106.50.30 29 1599.5d
1599.6 0.006 4.00 147 17.67 4 2 0.50
C107.50.30 29 1599.5d
1599.6 0.006 4.59 147 17.67 4 2 0.50
C108.50.30 29 1598.3d
1598.4 0.006 52.71 147 17.67 4 2 0.50
C109.50.30 29 1598.3c
1598.4 0.006 6.97 147 17.67 4 2 0.50
RC101.25.30 18 1419.8b
1419.8 0.000 4.87 102 21.98 5 2 0.50
RC102.25.30 18 1419.8b
1419.8 0.000 3.03 102 21.98 5 2 0.50
RC103.25.30 18 1419.8b
1419.8 0.000 3.00 102 21.98 5 2 0.50
RC104.25.30 18 1419.8b
1419.8 0.000 0.68 102 21.98 5 2 0.50
RC105.25.30 18 1419.8b
1419.8 0.000 4.86 141 30.39 5 3 0.50
RC106.25.30 18 1419.8b
1419.8 0.000 6.62 123 26.51 5 2 0.00
RC107.25.30 18 1419.8b
1419.8 0.000 1.60 102 21.98 5 2 0.50
RC108.25.30 18 1419.8b
1419.8 0.000 1.15 102 21.98 5 2 0.50
RC201.25.30 18 1419.8b
1419.8 0.000 13.64 141 30.39 5 3 0.50
RC202.25.30 18 1419.8b
1419.8 0.000 3.13 123 26.51 5 2 0.00
RC203.25.30 18 1419.8b
1419.8 0.000 5.16 102 21.98 5 2 0.50
RC204.25.30 18 1419.8b
1419.8 0.000 3788.68 102 21.98 5 2 0.50
RC205.25.30 18 1419.8b
1419.8 0.000 4.44 169 36.42 5 3 0.00
RC206.25.30 18 1419.8b
1419.8 0.000 7.73 102 21.98 5 2 0.50
RC207.25.30 18 1419.8b
1419.8 0.000 16.66 102 21.98 5 2 0.50
RC208.25.30 18 1419.8c
1419.8 0.000 2.80 141 30.39 5 3 0.50
RC101.50.30 33 2739.5c
2739.5 0.000 26.65 172 18.09 6 2 0.50
RC102.50.30 33 2739.5c
2739.5 0.000 16.68 172 18.09 6 2 0.50
RC103.50.30 33 2739.5c
2739.5 0.000 67.46 238 25.03 6 3 0.50
RC104.50.30 33 2739.5c
2739.5 0.000 110.02 172 18.09 6 2 0.50
RC105.50.30 33 2739.6c
2739.6 0.000 25.38 172 18.09 6 2 0.50
RC106.50.30 33 2739.5c
2739.5 0.000 42.26 172 18.09 6 2 0.50
RC107.50.30 33 2739.5c
2739.5 0.000 7.73 172 18.09 6 2 0.50
RC108.50.30 33 2739.5c
2739.5 0.000 3.95 172 18.09 6 2 0.50
RC201.50.30 33 2739.5c
2739.5 0.000 136.74 235 24.71 6 2 0.00
RC202.50.30 33 2739.5c
2739.5 0.000 18.04 172 18.09 6 2 0.50
RC203.50.30 33 2739.5c
2739.5 0.000 62.54 172 18.09 6 2 0.50
RC204.50.30 33 2739.5c
2739.5 0.000 15.94 172 18.09 6 2 0.50
RC205.50.30 33 2739.5c
2739.5 0.000 96.10 235 24.71 6 2 0.00
RC206.50.30 33 2739.5c
2739.5 0.000 27.41 235 24.71 6 2 0.00
RC207.50.30 33 2739.5c
2739.5 0.000 16.11 172 18.09 6 2 0.50
RC208.50.30 33 2739.5c
2739.5 0.000 3.00 172 18.09 6 2 0.50
Average 0.084 81.71 21.72 4.54 2.14 0.44
a
Intel Xeon E5-2640 v3, 2.6 GHz
b Desaulniers et al. [38]
c Archetti et al. [12]
d
Luo et al. [69]
Table 6.5: LKH-3 results for SDVRPTW with limited fleet on the Group C and RC instances of Solomon
[90]
83

Chapter 7
Conclusion and Future Work
This dissertation explores the applications of explicit enumeration schemes in combinatorial optimization,
with one of its emphasis on introducing a heuristic for split-delivery routing problems, leveraging the
foundational work of [29]. A demand splitting rule has been developed, yielding the minimum-size reduction of a split-delivery problem instance to an equivalent unsplittable instance of the capacitated vehicle
routing problem (CVRP). Inherent to the method is the capability to incorporate an assumed upper bound
on the number of vehicles servicing a given customer (¯k), thereby further reducing the size of the resulting non-split problem while preserving solution quality. Computational experiments for the split-delivery
vehicle routing problem (SDVRP) reveal that the a priori splitting rule proposed in this dissertation, when
utilized in conjunction with the LKH-3 solver without modifications to its metaheuristics, approaches the
performance levels of the best existing approaches [54, 88].
The real value of the enumeration schemes is further highlighted by their immediate applicability to
other extensions of the vehicle routing problem with split delivery. Computational experiments for the
split-delivery vehicle routing problem with time windows (SDVRPTW) show performance comparable to
that of state-of-the-art solutions. An interesting future direction for this work would involve considering
the full Pareto frontier defined by the trade-off between the size of a partition µ and the family of partitions
84

that it coalesces to. While this paper requires µ to be lossless, it might be sensible to use partitions that
are not entirely lossless, but very close to being so, if their cardinality is sufficiently small.
The major contribution of this dissertation is the provision of a tighter guaranteed upper bound of
the Euclidean travelling salesman constant, a figure that has not seen significant improvement for more
than six decades. It is explicitly proved that β2 < 0.9038, highlighting that the bound could be further
improved with additional computing power. The methodology employed in this dissertation also suggests
potential extensions to improve estimates for the corresponding constants βd for travelling salesman problems in higher dimensions d > 2, despite these constants appearing less frequently in the literature. The
enumeration and high-dimensional integration approach outlined is further applicable, with much simpler
antiderivative expressions, to bound constants for metrics such as the Manhattan metric, although these
topics also have not garnered much interest.
7.1 Lower bound of the travelling salesman constant
In contemplating future avenues of research, it is noteworthy that [46] achieved an improvement in the
lower bound of β2, leveraging the arguments put forth in [94]. This section proposes an innovative, albeit
challenging, enumeration-based approach aiming at further refining the lower bound.
First, define Lk as a random variable that denotes the length of a travelling salesman tour consisting
of k points uniformly distributed within a unit square, but with a special condition whereby any distance
covered along the boundary incurs no cost, making Lk a lower bound for the length of the true travelling
salesman tour:
ℓk := E[Lk] = Z 1
0
· · · Z 1
0
min
π∈Πk
X
k
i=1
d(xπ(i)
, xπ(i−1)) dV ,
where Πk represents the collection of all permutations π of the set {1, . . . , k}, with the condition that
π(0) = π(k). The distance metric d(xi
, xj ) indicates the shortest path between points xi and xj , which
85

(a) (b)
Figure 7.1: Figure 7.1a demonstrates that a unit square containing n points, each independently and
uniformly sampled, can be divided into n/k square cells. These cells are designed to permit cost-free
traversal along their boundaries, and each has a side length calculated as p
k/n. The quantity of points
per cell adheres to a Poisson distribution with mean k. Likewise, the unit square can be divided into n/k
hexagonal regions, facilitating more cost-free movements along their boundaries, as depicted in Figure
7.1b. In these hexagons, the distribution of points remains Poisson with mean k. The side length for each
hexagon is determined by 2
0.5
· 3
−0.75p
k/n.
may include traversal along the boundary. The integral, defined as the minimum among convex functions,
can be bounded below by considering a very large collection of tangent planes. By employing the methodologies outlined in Chapter 5, it is possible to identify hyperrectangles within which specific permutations
are universally optimal. Subsequently, the antiderivatives computed by FriCAS can be utilized to calculate
the lower bound improvement based on these identifications.
To understand the benefit of this approach, now consider n points X = {X1, . . . , Xn} that are independent and uniformly distributed within a unit square. The square is partitioned into n/k square cells,
denoted as Ci
. The distance matrix is adapted among all the points Xi to eliminate any cost associated
with travelling along the boundaries of the square cells, in contrast to traversing the perimeter of the entire
unit square as previously discussed. A lower bound on the length of this adjusted travelling salesman tour
86

Figure 7.2: The simulated lower bound of β2 was explored for both square and hexagonal configurations.
For each value of k, spanning from 5 to 39, a total of 10, 000 instances were simulated. The sample averages
and standard deviations derived from these simulations were then utilized to approximate the ℓi values for
each scenario. The yellow region highlights the 95% confidence interval for the average lower bound of
β2, revealing that achieving a value of k = 29 is necessary to exceed the best known lower bound reported
in [46].
87

involving all n points is termed the rooted dual [91], also known as the boundary functional [102]. This
modification is symbolized as TSP0, thereby establishing
TSP(X ) ≥
X
n/k
i=1
TSP0(X ∩ Ci).
Observe that when |X ∩ Ci
| = k, the entries TSP0(X ∩ Ci) follow the same distribution as the Lk’s,
scaled by a factor of p
k/n. Moreover, with a sufficiently large n, the quantity of points in each cell Ci
aligns with a Poisson distribution with mean k, and hence
β2 ≥
1
√
n
E
X
n/k
i=1
TSP0(X ∩ Ci) = 1
√
n
·
n
k
X∞
i=1
r
k
n
Pr(Zk = i)ℓi =
e
−k
√
k
X∞
i=1
k
i
i!
ℓi
,
where Zk denotes a Poisson random variable with mean k. Consequently, estimates of β2 can be refined by
lower bounding as many ℓi terms as the computational resources allow, while for i → ∞, simpler analytic
lower bounds, such as those based on nearest-neighbor distances, can be applied.
Figure 7.2 demonstrates how changes in the Poisson mean k influence the expected lower bound of
β2, referred to as β
sim
LB . These values are Monte Carlo estimates rather than actual lower bounds. The
analysis reveals that with an increase in k, there is a corresponding rise in β
sim
LB , which is consistent with
theoretical expectations. Although reaching these true lower bound values might exceed the limits of
present computational resources, this information is recorded for the sake of future exploration.
88

Bibliography
[1] Emile HL Aarts, Jan HM Korst, and Peter JM van Laarhoven. “A quantitative analysis of the
simulated annealing algorithm: A case study for the traveling salesman problem”. In: Journal of
Statistical Physics 50 (1988), pp. 187–206.
[2] Rafael E Aleman, Xinhui Zhang, and Raymond R Hill. “An adaptive memory algorithm for the
split delivery vehicle routing problem”. In: Journal of Heuristics 16 (2010), pp. 441–473.
[3] Kemal Altinkemer and Bezalel Gavish. “Heuristics for unequal weight delivery problems with a
fixed error guarantee”. In: Operations Research Letters 6.4 (1987), pp. 149–158.
[4] İ K Altınel and Temel Öncan. “A new enhancement of the Clarke and Wright savings heuristic for
the capacitated vehicle routing problem”. In: Journal of the Operational Research Society 56 (2005),
pp. 954–961.
[5] Andris Ambainis, Kaspars Balodis, Janis Iraids, Martins Kokainis, Krišj ¯ anis Pr ¯ usis, and ¯
Jevgenijs Vihrovs. “Quantum speedups for exponential-time dynamic programming algorithms”. ¯
In: Proceedings of the Thirtieth Annual ACM-SIAM Symposium on Discrete Algorithms. SIAM. 2019,
pp. 1783–1793.
[6] Sina Ansari, Mehmet Başdere, Xiaopeng Li, Yanfeng Ouyang, and Karen Smilowitz.
“Advancements in continuous approximation models for logistics and transportation systems:
1996–2016”. In: Transportation Research Part B: Methodological 107 (2018), pp. 229–252.
[7] David Applegate. “Concorde: A code for solving traveling salesman problems”. In: http://www.
tsp. gatech. edu/concorde. html (2003).
[8] David L Applegate, Robert E Bixby, Vašek Chvátal, William Cook, Daniel G Espinoza,
Marcos Goycoolea, and Keld Helsgaun. “Certification of an optimal TSP tour through 85,900
cities”. In: Operations Research Letters 37.1 (2009), pp. 11–15.
[9] David L Applegate, William J Cook, Vašek Chvátal, and Robert E Bixby. Concorde TSP Solver.
https://www.math.uwaterloo.ca/tsp/concorde/. 2015.
[10] Claudia Archetti, Nicola Bianchessi, and M Grazia Speranza. “Branch-and-cut algorithms for the
split delivery vehicle routing problem”. In: European Journal of Operational Research 238.3 (2014),
pp. 685–698.
89

[11] Claudia Archetti, Nicola Bianchessi, and Maria Grazia Speranza. “A column generation approach
for the split delivery vehicle routing problem”. In: Networks 58.4 (2011), pp. 241–254.
[12] Claudia Archetti, Mathieu Bouchard, and Guy Desaulniers. “Enhanced branch and price and cut
for vehicle routing with split deliveries and time windows”. In: Transportation Science 45.3 (2011),
pp. 285–298.
[13] Claudia Archetti, M Grazia Speranza, and Martin WP Savelsbergh. “An optimization-based
heuristic for the split delivery vehicle routing problem”. In: Transportation Science 42.1 (2008),
pp. 22–31.
[14] Claudia Archetti and Maria Grazia Speranza. “The split delivery vehicle routing problem: A
survey”. In: The vehicle routing problem: Latest advances and new challenges (2008), pp. 103–122.
[15] Claudia Archetti, Maria Grazia Speranza, and Alain Hertz. “A tabu search algorithm for the split
delivery vehicle routing problem”. In: Transportation science 40.1 (2006), pp. 64–73.
[16] Alper Atamtürk, Gemma Berenguer, and Zuo-Jun Shen. “A conic integer programming approach
to stochastic joint location-inventory problems”. In: Operations Research 60.2 (2012), pp. 366–381.
[17] Florin Avram and Dimitris Bertsimas. “The minimum spanning tree constant in geometrical
probability and under the independent model: a unified approach”. In: The Annals of Applied
Probability (1992), pp. 113–130.
[18] Andrew D Barbour, Lars Holst, and Svante Janson. Poisson approximation. Oxford University
Press, 1992.
[19] Jillian Beardwood, John H Halton, and John Michael Hammersley. “The shortest path through
many points”. In: Mathematical Proceedings of the Cambridge Philosophical Society. Vol. 55.
Cambridge University Press. 1959, pp. 299–327.
[20] José-Manuel Belenguer, MC Martinez, and E Mota. “A lower bound for the split delivery vehicle
routing problem”. In: Operations research 48.5 (2000), pp. 801–810.
[21] Luca Bertazzi, Bruce Golden, and Xingyin Wang. “The bin packing problem with item
fragmentation: a worst-case analysis”. In: Discrete Applied Mathematics 261 (2019), pp. 63–77.
[22] Jannis Blauth, Vera Traub, and Jens Vygen. “Improving the approximation ratio for capacitated
vehicle routing”. In: Mathematical Programming 197.2 (2023), pp. 451–497.
[23] Mourad Boudia, Christian Prins, and Mohamed Reghioui. “An effective memetic algorithm with
population management for the split delivery vehicle routing problem”. In: Hybrid Metaheuristics:
4th International Workshop, HM 2007, Dortmund, Germany, October 8-9, 2007. Proceedings 4.
Springer. 2007, pp. 16–30.
[24] Sylvia Boyd, René Sitters, Suzanne van der Ster, and Leen Stougie. “The traveling salesman
problem on cubic and subcubic graphs”. In: Mathematical Programming 144.1-2 (2014),
pp. 227–245.
90

[25] Vicente Campos, Angel Corberán, and Enrique Mota. “A scatter search algorithm for the split
delivery vehicle routing problem”. In: Advances in computational intelligence in transport, logistics,
and supply chain management (2008), pp. 137–152.
[26] Emmanuel J Candes, Michael B Wakin, and Stephen P Boyd. “Enhancing sparsity by reweighted
ℓ1 minimization”. In: Journal of Fourier analysis and applications 14.5 (2008), pp. 877–905.
[27] Marco Casazza and Alberto Ceselli. “Exactly solving packing problems with fragmentation”. In:
Computers & Operations Research 75 (2016), pp. 202–213.
[28] Sangit Chatterjee, Cecilia Carrera, and Lucy A Lynch. “Genetic algorithms and traveling
salesman problems”. In: European journal of operational research 93.3 (1996), pp. 490–510.
[29] Ping Chen, Bruce Golden, Xingyin Wang, and Edward Wasil. “A novel approach to solve the split
delivery vehicle routing problem”. In: International Transactions in Operational Research 24.1-2
(2017), pp. 27–41.
[30] Si Chen, Bruce Golden, and Edward Wasil. “The split delivery vehicle routing problem:
Applications, algorithms, test problems, and computational results”. In: Networks: An
International Journal 49.4 (2007), pp. 318–329.
[31] Nicos Christofides. “Worst-case analysis of a new heuristic for the travelling salesman problem”.
In: (1976).
[32] Geoff Clarke and John W Wright. “Scheduling of vehicles from a central depot to a number of
delivery points”. In: Operations research 12.4 (1964), pp. 568–581.
[33] William J Cook. “In pursuit of the traveling salesman”. In: In Pursuit of the Traveling Salesman.
Princeton University Press, 2011.
[34] William J Cook, David L Applegate, Robert E Bixby, and Vasek Chvatal. The traveling salesman
problem: a computational study. Princeton university press, 2011.
[35] Carlos Daganzo. Logistics systems analysis. Springer Science & Business Media, 2005.
[36] I Deif and L Bodin. “Extension of the Clarke and Wright algorithm for solving the vehicle routing
problem with backhauling”. In: Proceedings of the Babson conference on software uses in
transportation and logistics management. Babson Park, MA. 1984, pp. 75–96.
[37] Ulrich Derigs, B Li, and Ulrich Vogel. “Local search-based metaheuristics for the split delivery
vehicle routing problem”. In: Journal of the Operational Research Society 61.9 (2010),
pp. 1356–1364.
[38] Guy Desaulniers. “Branch-and-price-and-cut for the split-delivery vehicle routing problem with
time windows”. In: Operations research 58.1 (2010), pp. 179–192.
[39] Moshe Dror and Pierre Trudeau. “Stochastic vehicle routing with modified savings algorithm”. In:
European Journal of Operational Research 23.2 (1986), pp. 228–235.
91

[40] Leonard Few. “The shortest path and the shortest road through n points”. In: Mathematika 2.2
(1955), pp. 141–144.
[41] C-N Fiechter. “A parallel tabu search algorithm for large traveling salesman problems”. In:
Discrete Applied Mathematics 51.3 (1994), pp. 243–267.
[42] Steven R Finch. Mathematical constants. Cambridge university press, 2003.
[43] Anna Franceschetti, Ola Jabali, and Gilbert Laporte. “Continuous approximation models in freight
distribution management”. In: Top 25 (2017), pp. 413–433.
[44] FriCAS team. FriCAS—an advanced computer algebra system. http://fricas.github.io/. 2021.
[45] Alan Frieze and Wesley Pegden. “Separating subadditive Euclidean functionals”. In: Proceedings of
the forty-eighth annual ACM symposium on Theory of Computing. 2016, pp. 22–35.
[46] Julia Gaudio and Patrick Jaillet. “An improved lower bound for the Traveling Salesman constant”.
In: Operations Research Letters 48.1 (2020), pp. 67–70.
[47] Shayan Oveis Gharan, Amin Saberi, and Mohit Singh. “A randomized rounding approach to the
traveling salesman problem”. In: 2011 IEEE 52nd Annual Symposium on Foundations of Computer
Science. IEEE. 2011, pp. 550–559.
[48] Fred Glover. “Tabu search—part I”. In: ORSA Journal on computing 1.3 (1989), pp. 190–206.
[49] John Grefenstette, Rajeev Gopal, Brian Rosmaita, and Dirk Van Gucht. “Genetic algorithms for
the traveling salesman problem”. In: Proceedings of the first International Conference on Genetic
Algorithms and their Applications. Psychology Press. 2014, pp. 160–168.
[50] Christopher Groër. Parallel and serial algorithms for vehicle routing problems. University of
Maryland, College Park, 2008.
[51] Christopher Groër, Bruce Golden, and Edward Wasil. “A library of local search heuristics for the
vehicle routing problem”. In: Mathematical Programming Computation 2 (2010), pp. 79–101.
[52] Gregory Gutin and Abraham P Punnen. The traveling salesman problem and its variations. Vol. 12.
Springer Science & Business Media, 2006.
[53] Mordecai Haimovich and Alexander HG Rinnooy Kan. “Bounds and heuristics for capacitated
routing problems”. In: Mathematics of operations Research 10.4 (1985), pp. 527–542.
[54] Pengfei He and Jin-Kao Hao. “General edge assembly crossover-driven memetic search for split
delivery vehicle routing”. In: Transportation Science 57.2 (2023), pp. 482–511.
[55] Michael Held and Richard M Karp. “A dynamic programming approach to sequencing problems”.
In: Journal of the Society for Industrial and Applied mathematics 10.1 (1962), pp. 196–210.
[56] Michael Held and Richard M Karp. “The traveling-salesman problem and minimum spanning
trees”. In: Operations Research 18.6 (1970), pp. 1138–1162.
92

[57] Michael Held and Richard M Karp. “The traveling-salesman problem and minimum spanning
trees: Part II”. In: Mathematical programming 1.1 (1971), pp. 6–25.
[58] Keld Helsgaun. “An effective implementation of the Lin–Kernighan traveling salesman heuristic”.
In: European journal of operational research 126.1 (2000), pp. 106–130.
[59] Keld Helsgaun. “An extension of the Lin-Kernighan-Helsgaun TSP solver for constrained
traveling salesman and vehicle routing problems”. In:
http://akira.ruc.dk/~keld/research/LKH-3/LKH-3_REPORT.pdf (2017).
[60] Keld Helsgaun. “General k-opt submoves for the Lin–Kernighan TSP heuristic”. In: Mathematical
Programming Computation 1 (2009), pp. 119–163.
[61] Keld Helsgaun. LKH-3. http://akira.ruc.dk/~keld/research/LKH-3/. 2019.
[62] Fredrik Johansson, Vinzent Steinberg, Sergey B Kirpichev, Kristopher L Kuhlman, Aaron Meurer,
Ondřej Čertık, Case Van Horsen, Paul Masson, Juan Arias de Reyna, Timo Hartmann, et al.
“mpmath: a Python library for arbitrary-precision floating-point arithmetic”. In: Zenodo (2013).
[63] David S Johnson, Lyle A McGeoch, and Edward E Rothberg. “Asymptotic experimental analysis
for the Held-Karp traveling salesman bound”. In: Proceedings of the 7th Annual ACM-SIAM
Symposium on Discrete Algorithms. Vol. 341. ACM Press San Francisco. 1996, p. 350.
[64] Bo Jones and John Gunnar Carlsson. “Minimum size generating partitions and their application
to demand fulfillment optimization problems”. In: arXiv preprint arXiv:1909.09363 (2019).
[65] Anna R Karlin, Nathan Klein, and Shayan Oveis Gharan. “A (slightly) improved approximation
algorithm for metric TSP”. In: Proceedings of the 53rd Annual ACM SIGACT Symposium on Theory
of Computing. 2021, pp. 32–45.
[66] Jon M Kleinberg. “Single-source unsplittable flow”. In: Proceedings of 37th Conference on
Foundations of Computer Science. IEEE. 1996, pp. 68–77.
[67] Jooyoung Lee and MooYoung Choi. “Optimization by multicanonical annealing and the traveling
salesman problem”. In: Physical Review E 50.2 (1994), R651.
[68] Shen Lin and Brian W Kernighan. “An effective heuristic algorithm for the traveling-salesman
problem”. In: Operations research 21.2 (1973), pp. 498–516.
[69] Zhixing Luo, Hu Qin, Wenbin Zhu, and Andrew Lim. “Branch and price and cut for the
split-delivery vehicle routing problem with time windows and linear weight-related cost”. In:
Transportation Science 51.2 (2017), pp. 668–687.
[70] Miroslaw Malek, Mohan Guruswamy, Mihir Pandya, and Howard Owens. “Serial and parallel
simulated annealing and tabu search algorithms for the traveling salesman problem”. In: Annals
of Operations Research 21 (1989), pp. 59–84.
[71] Tobias Mömke and Ola Svensson. “Approximating graphic TSP by matchings”. In: 2011 IEEE 52nd
Annual Symposium on Foundations of Computer Science. IEEE. 2011, pp. 560–569.
93

[72] Marcin Mucha. “-Approximation for Graphic TSP”. In: Theory of computing systems 55.4 (2014),
pp. 640–657.
[73] Paul Abraham Mullaseril. “Capacitated rural postman problem with time windows and split
delivery”. In: The University of Arizona PhD Dissertation (1997).
[74] Yuichi Nagata. “New EAX crossover for large TSP instances”. In: International Conference on
Parallel Problem Solving from Nature. Springer. 2006, pp. 372–381.
[75] HL Ong and HC Huang. “Asymptotic expected performance of some TSP heuristics: an empirical
evaluation”. In: European Journal of Operational Research 43.2 (1989), pp. 231–238.
[76] Gizem Ozbaygin, Oya Karasan, and Hande Yaman. “New exact solution approaches for the split
delivery vehicle routing problem”. In: EURO Journal on Computational Optimization 6.1 (2018),
pp. 85–115.
[77] Allon G Percus and Olivier C Martin. “Finite size and dimensional dependence in the Euclidean
traveling salesman problem”. In: Physical Review Letters 76.8 (1996), p. 1188.
[78] Laurent Perron and Vincent Furnon. OR-Tools. https://developers.google.com/optimization/.
2020.
[79] Tantikorn Pichpibul, Ruengsak Kawtummachai, et al. “A heuristic approach based on
clarke-wright algorithm for open vehicle routing problem”. In: The scientific world journal 2013
(2013).
[80] Tantikorn Pichpibul and Ruengsak Kawtummachai. “An improved Clarke and Wright savings
algorithm for the capacitated vehicle routing problem”. In: ScienceAsia 38.3 (2012), pp. 307–318.
[81] C Redmond and JE Yukich. “Limit theorems and rates of convergence for Euclidean functionals”.
In: The Annals of Applied Probability (1994), pp. 1057–1073.
[82] Gerhard Reinelt. “TSPLIB—A traveling salesman problem library”. In: ORSA journal on computing
3.4 (1991), pp. 376–384.
[83] William Stein, David Joyner, David Kohel, John Cremona, and Burçin Eröcal. SageMath, the Sage
Mathematics Software System. https://www.sagemath.org/. 2022.
[84] Matteo Salani and Ilaria Vacca. “Branch and price for the vehicle routing problem with discrete
split deliveries and time windows”. In: European Journal of Operational Research 213.3 (2011),
pp. 470–477.
[85] Stefan Schröder. jsprit. https://github.com/graphhopper/jsprit. 2020.
[86] Anatolii Ivanovich Serdyukov. “Some extremal bypasses in graphs”. In: Diskretnyi Analiz i
Issledovanie Operatsii 17 (1978), pp. 76–79.
[87] David B Shmoys and David P Williamson. “Analyzing the Held-Karp TSP bound: A monotonicity
property with application”. In: Information Processing Letters 35.6 (1990), pp. 281–285.
94

[88] Marcos Melo Silva, Anand Subramanian, and Luiz Satoru Ochi. “An iterated local search heuristic
for the split delivery vehicle routing problem”. In: Computers & Operations Research 53 (2015),
pp. 234–249.
[89] Geoffrey De Smet. Optaplanner. http://www.optaplanner.org. 2020.
[90] Marius M Solomon. “Algorithms for the vehicle routing and scheduling problems with time
window constraints”. In: Operations research 35.2 (1987), pp. 254–265.
[91] J Michael Steele. Probability theory and combinatorial optimization. SIAM, 1997.
[92] J Michael Steele. “Subadditive Euclidean functionals and nonlinear growth in geometric
probability”. In: The Annals of Probability (1981), pp. 365–376.
[93] David M Stein. Scheduling Dial-A-Ride Transportation Systems: An Asymptotic Approach. Tech. rep.
HARVARD UNIV CAMBRIDGE MASS DIV OF ENGINEERING and APPLIED PHYSICS, 1977.
[94] Stefan Steinerberger. “New bounds for the traveling salesman constant”. In: Advances in Applied
Probability 47.1 (2015), pp. 27–36.
[95] L Fejes Tóth. “Über einen geometrischen Satz”. In: Mathematische Zeitschrift 46.83-85 (1940),
p. 463.
[96] Christine L Valenzuela and Antonia J Jones. “Estimating the Held-Karp lower bound for the
geometric TSP”. In: European journal of operational research 102.1 (1997), pp. 157–175.
[97] Santosh S Venkatesh. The theory of probability: Explorations and applications. Cambridge
University Press, 2013.
[98] Samuel Verblunsky. “On the shortest path through a number of points”. In: Proceedings of the
American Mathematical Society 2.6 (1951), pp. 904–913.
[99] Weijun Xie and Yanfeng Ouyang. “Optimal layout of transshipment facility locations on an
infinite homogeneous plane”. In: Transportation Research Part B: Methodological 75 (2015),
pp. 74–88.
[100] Zhengguang Xie and Jianping Hu. “Rewighted l1-minimization for sparse solutions to
underdetermined linear systems”. In: 2013 6th International Congress on Image and Signal
Processing (CISP). Vol. 3. IEEE. 2013, pp. 1660–1664.
[101] Julien Yu. BHH upper bound. https://github.com/yyu56253/BHH_upper_bound. 2023.
[102] Joseph E Yukich. Probability theory of classical Euclidean optimization problems. Springer, 2006.
95

Chapter 8
Appendix
This appendix details the process of computing CB,·,·
, as outlined in Algorithm 3. For the initial k − 1
summands, assign π(i) = v, π(i + 1) = w. The calculation of CB,v,w proceeds as follows:
Z
B
e
−xk

xv − xw
h
2
(uv − uw)

dV
=
Z t1
s1
· · · Z tk
sk
Z b0
a0
· · · Z bk
ak
e
−xk

xv − xw
h
2
(uv − uw)

duk . . . du0 dxk . . . dx1
= C
Z t1
s1
· · · Z tk
sk
Z bw
aw
Z bv
av
e
−xk

xv − xw
h
2
(uv − uw)

duv duw dxk . . . dx1
= C
Z t1
s1
· · · Z tk−1
sk−1
(e
−sk − e
−tk )
Z bw
aw
Z bv
av

xv − xw
h
2
(uv − uw)

duv duw dxk−1 . . . dx1
= CB,v,w Z tw
sw
Z tv
sv
Z bw
aw
Z bv
av

xv − xw
h
2
(uv − uw)

duv duw dxv dxw.
Here, CB,v,w = (e
−sk − e
−tk )Πk−1
i=1,i̸=v,w
(ti − si)Πk
i=0,i̸=v,w(bi − ai), depending on the values of v and w,
as well as the bounds of hyperrectangle B.
96

For the final summand where i = k, since π(k) = k, the variable xk is included within the norm. Let
π(k − 1) = v. The detailed calculation for CB,v,k is as follows:
Z
B
e
−xk

xv − xk
h
2
(uv − uk)

dV
=
Z t1
s1
· · · Z tk
sk
Z b0
a0
· · · Z bk
ak
e
−xk

xv − xk
h
2
(uv − uk)

duk . . . du0 dxk . . . dx1
= C
Z t1
s1
· · · Z tk
sk
Z bk
ak
Z bv
av
e
−xk

xv − xk
h
2
(uv − uk)

duv duk dxk . . . dx1
= CB,v,k Z tk
sk
Z tv
sv
Z bk
ak
Z bv
av
e
−xk

xv − xk
h
2
(uv − uk)

duv duk dxv

dxk.
In this scenario, CB,v,k = Πk−1
i=1,i̸=v
(ti − si)Πk−1
i=0,i̸=v
(bi − ai), which is contingent upon v and the bounds
of hyperrectangle B.
97

Abstract (if available)

Abstract Enumeration schemes play a pivotal role in addressing the challenge of complex optimization problems. This dissertation embarks on a detailed exploration of two vital topics in combinatorial optimization, refining the upper bound of the Euclidean Travelling Salesman Problem (TSP), and providing novel insights into efficient resource allocation in the Split-Delivery Vehicle Routing Problem (SDVRP) and the Split-Delivery Vehicle Routing Problem with Time Windows (SDVRPTW). The integration of these two themes offers a novel perspective in the field of combinatorial optimization, combining theoretical innovation with practical problem-solving approaches.

At the outset, this dissertation focuses on a fundamental challenge in combinatorial optimization: improving the upper bound of the Euclidean TSP constant. Let X1, X2, . . . , Xn be n independent and uniformly distributed random points in a compact region R ⊂ R2 of area 1. Let TSP(X1, . . . , Xn) denote the length of the optimal Euclidean travelling salesman tour that traverses all these points. The classical Beardwood-Halton-Hammersley theorem proves the existence of a universal constant β2 such that TSP(X1, . . . , Xn)/√n → β2 almost surely, which satisfies 0.625 ≤ β2 ≤ 0.92116027. This research presents a computer-aided proof using numerical quadrature and tree-based machine learning that β2 ≤ 0.9038, marking a significant stride in optimizing the Euclidean TSP and underscoring the critical role of enumeration schemes in unraveling combinatorial optimization problems. This advancement also opens the door for further improvements with the evolution of computational hardware over time. In addition to the explicit upper bound improvement, this research also suggests a novel enumeration-based direction for enhancing the lower bound of the Euclidean TSP constant, although it appears computationally intractable at present.

Simultaneously, the dissertation ventures into the intricate territory of splittable demand in resource allocation problems, a challenging aspect of combinatorial optimization. It presents an innovative a priori splitting rule, minimizing the number of pieces while maintaining the feasibility of all possible demand splitting patterns. This strategy effectively transforms splittable problems into unsplittable ones, leveraging the advanced heuristics developed for the latter. This part of the research demonstrates the versatility and efficacy of enumeration schemes in devising practical solutions, as validated by computational experiments on large-scale instances, including those involving vehicle routing problems with time-windows.

On the whole, this dissertation underscores the indispensable role of enumeration schemes in tackling combinatorial optimization problems, paving the way for innovative solutions.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Dynamic programming-based algorithms and heuristics for routing problems

PDF

The warehouse traveling salesman problem and its application

PDF

Improving decision-making in search algorithms for combinatorial optimization with machine learning

PDF

A continuous approximation model for the parallel drone scheduling traveling salesman problem

PDF

Continuous approximation for selection routing problems

PDF

Asymptotic analysis of the generalized traveling salesman problem and its application

PDF

Package delivery with trucks and UAVs

PDF

Applications of topological data analysis to operational research problems

PDF

Continuous approximation formulas for cumulative routing optimization problems

PDF

A stochastic conjugate subgradient framework for large-scale stochastic optimization problems

PDF

Landscape analysis and algorithms for large scale non-convex optimization

PDF

Artificial Decision Intelligence: integrating deep learning and combinatorial optimization

PDF

Applications of Wasserstein distance in distributionally robust optimization

PDF

Train routing and timetabling algorithms for general networks

PDF

Modern nonconvex optimization: parametric extensions and stochastic programming

PDF

Mixed-integer nonlinear programming with binary variables

PDF

Advancing distributed computing and graph representation learning with AI-enabled schemes

PDF

Congestion reduction via private cooperation of new mobility services

PDF

Integer optimization for analytics in high stakes domain

PDF

Scalable optimization for trustworthy AI: robust and fair machine learning

Asset Metadata

Creator Yu, Yue (author)

Core Title Applications of explicit enumeration schemes in combinatorial optimization

School Viterbi School of Engineering

Degree Doctor of Philosophy

Degree Program Industrial and Systems Engineering

Degree Conferral Date 2024-05

Publication Date 04/17/2024

Defense Date 04/10/2024

Publisher Los Angeles, California (original), University of Southern California (original), University of Southern California. Libraries (digital)

Tag combinatorial optimization,enumeration,high dimension integral,OAI-PMH Harvest,splitting,Travelling Salesman problem

Format theses (aat)

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Carlsson, John (committee member), Razaviyayn, Meisam (committee member), Savla, Ketan (committee member)

Creator Email julienyue@berkeley.edu,yyu56253@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-oUC113877822

Unique identifier UC113877822

Identifier etd-YuYue-12823.pdf (filename)

Legacy Identifier etd-YuYue-12823

Document Type Dissertation

Format theses (aat)

Rights Yu, Yue

Internet Media Type application/pdf

Type texts

Source 20240417-usctheses-batch-1141 (batch), University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the author, as the original true and official version of the work, but does not grant the reader permission to use the work if the desired use is covered by copyright. It is the author, as rights holder, who must provide use permission if such use is covered by copyright.

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA

Repository Email cisadmin@lib.usc.edu