 A Data-Driven Methodology for Last-Mile Routing

Learn Global and Optimize Local: A Data-Driven Methodology for Last-Mile
Routing

Mayukh Ghosh1, Alex Kuiper1, Roshan Mahes1,2, Donato Maragno1 1 Amsterdam
Business School, University of Amsterdam, 1018 TV Amsterdam, the Netherlands
2 Korteweg-de Vries Institute for Mathematics, University of Amsterdam, 1098
XG Amsterdam, the Netherlands
<EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS><EMAIL_ADDRESS>
In last-mile routing, the task of finding a route is often framed as a
Traveling Salesman Problem to minimize travel time and associated cost.
However, solutions stemming from this approach do not match the realized paths
as drivers deviate due to navigational considerations and preferences. To
prescribe routes that incorporate this tacit knowledge, a data-driven model is
proposed that aligns well with the hierarchical structure of delivery data
wherein each stop belongs to a zone—a geographical area. First, on the global
level, a zone sequence is established as a result of a minimization over a
cost matrix which is a weighted combination of historical information and
distances (travel times) between zones. Subsequently, within zones, sequences
of stops are determined, such that, integrated with the predetermined zone
sequence, a full solution is obtained.

The methodology is particularly promising as it propels itself within the top-
tier of submissions to the Last-Mile Routing Research Challenge, while it
maintains an elegant decomposition that ensures a feasible implementation into
practice. The concurrence between prescribed and realized routes underpins the
adequateness of a hierarchical breakdown of the problem, and the fact that
drivers make a series of locally optimal decisions when navigating.
Furthermore, experimenting with the balance between historical information and
distance exposes that historic information is pivotal in deciding a starting
zone of a route. The experiments also reveal that at the end of a route,
historical information can best be discarded, making the time it takes to
return to the station the primary concern.

last-mile logistics, route prediction, navigation, traveling salesman problem,
data-driven optimization

##  1 Introduction

Last-mile routing is often the shortest, yet most complex and expensive leg in
distribution systems. It accounts for about 28% of the total costs of
transportation from the distribution center to the final customer (Goodman
2005, Gevaers, Van de Voorde, and Vanelslander 2011). Last-mile operations
also have a large impact on the livability, because of traffic congestion and
pollution it creates (Taniguchi and Thompson 2002, Chopra 2003). This pressure
is set to increase, because the trend of urbanization continues; around 68% of
the global population will reside in urban areas by 2050, according to a
report by the United Nations (2018). Currently e-commerce sales are estimated
at about $709.78 billion, which is projected to increase as the coronavirus
has accelerated a channel shift to e-commerce, see eMarketer (2020). The same
trends are echoed by a recent report of the World Economic Forum (Deloison et
al. 2020), which further points out that also the number of same-day
deliveries will grow due to customers ever-increasing expectations, amplifying
the importance of last-mile operations and its impact on the ecosystem.

Considering parcel delivery in last-mile operations, the main goal is to
ensure that products arrive in good shape. In addition, it should be on time
(Gunasekaran, Patel, and McGaughey 2004), i.e., the time windows that have
been communicated will be met. As routes require a multitude of stops to be
visited, it is critical that a proposed route needs to be followed up as
deviations can incur redeliveries or dissatisfied customers (by not meeting
time slots), unforeseen time losses along the route and additional operational
expenses (Boyer, Prud’homme, and Chung 2009, Gevaers, Van de Voorde, and
Vanelslander 2014).

Despite the benefits of adhering to a proposed ‘optimal’ route are clear, in
practice—reported by managers in last-mile delivery (e.g., Amazon)—drivers
often deviate from proposed routes. This discrepancy between theoretical route
planning and real-life route execution is an important gap in the literature
on last-mile routing, and has remained an open research question, because of
the typical inaccessibility of real-life data. Recently, Amazon in combination
with MIT’s Center for Transportation & Logistics hosted a challenge (Amazon
and MIT Center for Transportation & Logistics 2021a), wherein by means of
offering a vast routing dataset, researchers around the world have been
invited to bridge this gap. In this work we answer that call by providing a
complete methodology that proposes routes that adequately reflect drivers’
tacit knowledge and preferences by learning from historic route data.

In the US, Amazon steadily holds the number one position as the largest
e-commerce retailer having a market share of 38.0% (eMarketer 2020), and
inherently is responsible for the majority of traffic associated to last-mile
delivery. Therefore, by applying the methodology on the data associated to
Amazon will set an example for many other organizations facing the bewildering
setting of last-mile delivery. Although we consider it as a predictive
problem, which helps uncovering on a conceptual level the driver’s decision
process when determining a route, we emphasize that the approach can also be
used to prescribe routes, which are more feasible from a driver’s perspective.

Figure 1: Examples of two realized routes by using OpenStreetMap (2021). The
numbers represent the order in which the stops are visited and the colors
represent the different zones within the routes.

To consider the problem at a more detailed level, two delivery drivers’ routes
are given in Figure 1 as an increasing sequence of numbers; each number
represents a stop. Note that not all stops are shown as for readability
purposes. Along the route the driver traverses through different zones, given
by their zone ids, which are (preset) geographical areas. Comparing the routes
in the left (more rectangular grid) and the right map, it is observed that the
size of a zone ranges from a part of a street to a set of multiple streets.

Leveraging historical routes to learn recurrent patterns is feasible on a zone
level rather than at a stop level. This is due to the fact that a zone appears
in multiple routes while a stop does not (or rarely). The intrinsic importance
of a zone sequence depends on how the zone is established. Ideally, a zone
should be constructed such that every stop in it is visited before moving to
the next zone (as for example in Figure 1). However, with this presumption,
classical heuristics and optimization models fall short to exploit this
information. Two classical examples are illustrated in Figure 2. On the left,
we find the route obtained by applying the nearest neighbor heuristic, which
is sometimes argued to mimic the driver’s behavior, while on the right the
route that minimizes distance, as the solution of the corresponding Traveling
Salesman Problem (TSP).

On top of the condition that drivers visit all stops in a zone before moving
to another, the challenge is to determine the sequence of zones a driver will
follow. In Figure 2, there are six possible permutations of zone sequences,
but in real-life settings (Figure 1) there are easily 30 zones leading to 30!
unique zone sequences. The sequence might be based on information at ground
level, navigational considerations, or merely a preference, but most
importantly an approach should be able to learn this tacit information from
realized routes.

(a) Nearest Neighbor solution.

(b) TSP solution.

Figure 2: An example of a heuristic (a) and the solution when using an
optimization model (b).

The contributed methodology divides the problem into two parts: on the zone
level, it aims to infer the sequence of zones after which, within each zone,
the sequence of stops is determined. It stems from the idea that first a
routing problem is abstractified into zones, for which a sequence is
established by solving a TSP formulation that balances historical zone
transition information, resembling the tacit knowledge, and distances. So the
methodology learns on a global level, which aligns with the intuition that
humans congregate route sequences on a higher level than remembering exact
stop sequences. Secondly, the problem in each zone is of relatively low
complexity and thus a driver is capable to determine an optimal sequence out
of the set of stops that minimizes distance or travel time. In the end, the
performance of the methodology confirms the competency of the approach.

In the next section, we provide the literature review, which outlines
important elements to our methodology. Section 3 describes the methodology,
and is split into learning and prediction. In Section 4, we introduce the
dataset on which the methodology will be applied and provide the route
prediction metrics aligned with the challenge. Using the case study data we
carry out various experiments and demonstrate the methodology. Furthermore, we
compare the performance to the aforementioned benchmarks as well as to the
submissions of the competition. In Section 5, we present our concluding
remarks and implications to practice, followed by a discussion and some
suggestions for future research.

##  2 Literature Review

Solving a last-mile routing problem can be translated to a TSP, i.e., a set of
stops should be visited—starting from and ending at the station. However, in
this research the goal is to predict a driver’s route using historic data.
Therefore, in this pursuit, elements of different disciplines are combined.
More specifically, we consider classical TSP optimization frameworks,
hierarchical structures as an approximation (such as cluster first, routing
second), human approaches to solving TSPs, and lastly, learning from route
data. In relation to our work, we elaborate further on these themes below.

###  2.1 Traveling Salesman Problem

The TSP with Euclidean distance is notorious for being hard and is proven to
be NP-complete (Papadimitriou 1977), but for finding the optimal or near-
optimal tour there is a rich literature (Lawler 1985, Applegate et al. 2011).
Also, several variations of the standard TSP have been introduced in
literature, to name a few: an open TSP (OTSP) where the driver does not
necessarily return to the depot; TSP with time windows to meet (Gendreau et
al. 1998); additionally allowing pick-ups such that the capacity of the truck
is not exceeded (Gendreau, Laporte, and Vigo 1999); or the covering Salesman
problem (CSP) (Current and Schilling 1989), where a minimal tour is found such
that all stops are in the vicinity of a stop on the tour (e.g., dropping goods
at neighbors).

Acquainted to the TSP is its generalization known as Vehicle Routing Problem
(VRP); see Dantzig and Ramser (1959) for its formulation and Bräysy and
Gendreau (2005a), Montoya-Torres et al. (2015) for recent surveys on multiple
depots VRPs with variants that include time windows, split delivery,
heterogeneous fleet, periodic deliveries, and pick-up. Also, various studies
have considered the ecological footprint as presented in the literature review
given by Lin et al. (2014). Finally, also for the more general VRP, numerous
heuristics and metaheuristics have been proposed over the last decades (Bräysy
and Gendreau 2005b, Speranza and Archetti 2014, Purkayastha et al. 2020).

Both TSP and VRP have seen many different integer programming formulations
(Orman and Williams 2007). In our model, we adopt the Miller-Tucker-Zemlin
(MTZ) formulation to solve the TSP, introduced for the first time by Miller,
Tucker, and Zemlin (1960). The MTZ formulation has the advantage of being
intuitive and straightforward to implement and works well when the number of
stops (variables) to consider is relatively small, fitting our instances. As a
downside, it has been proven to lead to a weak relaxation, leading to higher
computation times, see e.g., Campuzano, Obreque, and Aguayo (2020) for a
discussion. In our case, as we break down the full problem into smaller ones,
this is not an issue, but as desired any other TSP formulation can be adopted.

###  2.2 Hierarchical Structure

To reduce the computational complexity, the problem can be broken down into
smaller instances, which are computationally much less involved, and connect
these together. For example, in the seminal work by Karp (1977), the proposed
partitioning algorithms are asymptotically optimal heuristics in the case of
uniformly distributed stop locations, i.e., the error tends to zero with
probability one as the length of a route (i.e., the number of stops)
increases.

Both Liao and Liu (2018) and Jiang et al. (2014) decompose the problem with
small-scale nodes by relying on clustering algorithms. Each subproblem is
subsequently optimized and the center nodes of the subproblems constitute a
TSP in itself. Connecting all local tours in the order of the upper layer
problem generates approximative solutions with significantly reduced
computation times. Such an approach is also found from an empirical point of
view, see Vickers et al. (2003).

Moreover, evidence to support the use of a hierarchical structure in practice
is found in Graham, Joshi, and Pizlo (2000). They show that solution methods
originating from artificial intelligence or operations research algorithms are
insufficiently capable to mimic the human approach of solving TSPs. They
introduce a hierarchical model by means of a pyramid algorithm on the visual
representation. This algorithm is capable to render human-alike solutions to
the various TSP instances where classical algorithms fail in this task, see
also Pizlo et al. (2006) for a more refined algorithm. More importantly, the
aforementioned works hint at the use of a hierarchical breakdown, because of
lower computational complexity and being concurrent with practice.

###  2.3 Human Navigation

The TSP lends itself to a broad range of experimental studies as its goal of
minimization of the route is easy to understand and visualize. Many
experimental studies have been devoted to visual versions of the TSP
(MacGregor and Ormerod 1996, Van Rooij, Stege, and Schactman 2003, Vickers et
al. 2003). MacGregor and Ormerod (1996) show that humans outperform well-known
TSP heuristics and are in small problem sizes capable to be in 1% from
optimality.

Humans rely on various tactics to generate near-optimal solutions for the TSP.
The potential of these tactics is further demonstrated by the fact that the
time needed increases only in a linear fashion compared to the problem size
(Pizlo et al. 2006, Dry et al. 2006). One of these tactics is to consider the
problem first globally, considering the tour that visits all ‘exterior’
points, to which interior points are inserted; this resembles the so-called
convex hull approach (MacGregor and Ormerod 1996, Vickers et al. 2003). Argued
as the underlying motivation for the convex hull approach, is the avoidance of
crossings in a tour (Van Rooij, Stege, and Schactman 2003). The motivation
behind this tactic is the intuition that a cross in a tour is suboptimal,
which is even a fact for metric TSPs. In an OTSP, where the solution is not a
tour, the starting point can have a profound impact on the performance, see
Sengupta, Mariescu-Istodor, and Fränti (2018). However, they also report that
humans are quite capable to select a ‘right’ starting point.

Wiener, Ehbauer, and Mallot (2009) find in a series of experiments that, when
navigating, a coarse route is stipulated first. This route visits a set of
‘regions’ after which it is ‘optimized’ on a detailed level along the way.
Furthermore, if the problem size increases, the problem is divided into more
regions to make the problem manageable and approachable. Such a global-to-
local approach echoes the hierarchical decomposition of first determining the
zone sequence on a global level, and subsequently a series of local problems
to determine the sequence in which stops are to be visited in each zone.
Additionally, another experiment from Wiener and Mallot (2003) reveals that a
segmentation into zones affects the route planning and navigational behavior
as it primes a driver to approach the problem from such a structure. A similar
effect can be expected when drivers are presented a route, for example a
solution from routing software as in the Amazon challenge, from which we take
the data to demonstrate our methodology (Amazon and MIT Center for
Transportation & Logistics 2021a).

###  2.4 Route Prediction

Typically, the problem of finding a sequence of stops is positioned in the
field of combinatorial optimization, but in this research the goal of find an
optimal sequence to minimize an objective function does not take center stage.
Here, the goal is how route data can be used to predict routes that drivers
(will) follow; for an example where learning takes place in order to find the
optimal route, see Kool, van Hoof, and Welling (2019).

The capability to learn and predict (part of) the routes is demonstrated in
Krumm (2008). The model is trained from drivers’ long-term trip history using
GPS data. It uses a Markovian approach, i.e., on the basis of the last road
segment, which can also be adjusted to incorporate more history, it predicts
the next segment the driver will take. In the same vein, Ye et al. (2015)
propose a route prediction method based on a hidden Markov model that can
accurately predict an entire route early in the trip. Wang et al. (2015) also
employ a Markov model; their algorithm relies on a probability transition
matrix that is developed to represent the knowledge of the driver’s preferred
links and routes. For the VRP, Canoy and Guns (2019) show the potential of
using a Markov model in an optimization framework. By constructing a
transition probability matrix, based on historical data, and by exploiting the
VRP structure, they render solutions that resemble actual route plans much
better than relying on a distance metric.

The previous works demonstrate that a Markov model is a powerful tool to learn
from historical data and to use these in predictions. We adopt such a model by
using historical information, which is weighted against a distance metric. The
full methodology is detailed in the next section.

##  3 Methodology

The methodology that we propose relies on the fourfold of elements reviewed in
Section 2. We impose a hierarchical structure on two levels: finding the zone
sequence on a global level and finding the stop sequence within each zone
locally. Such an approach is in line with evidence from human navigation, and
more importantly matches the structure of typical routing data—each stop
belongs to a zone. We learn on the global level the preference to traverse
from one zone to another by adopting a Markov model, which we combine with the
distances between zones to form the cost matrix for the global problem.
Relying on standard methods, we solve the TSP on the global level to determine
the zone sequence, and a series of TSPs on the local level, which mimics the
driver capabilities to find near-optimal solutions in small instances.
Finally, by patching the local solutions, we deliver a full prediction of the
route.

To learn the drivers’ preference over zones we use the fact that each route
starts at a station $s$—also commonly known as depots, where $s$ belongs to
the set of stations $\mathcal{S}$. Then for each stop in a route, the data
consists of the GPS coordinates, which can be used to compute distances
between stops, but the methodology can also be used when the costs of
traversing from one stop to another are given, e.g., (expected) travel times.
Furthermore, each stop within the route should be given an identifier, which
relates the stop to belong to a geographical area, i.e., a zone. The zones
should not overlap and therefore they form a segmentation. As a side note, the
zone ids should be consistent within the set of routes in the dataset.
Specifically in relation to zones there are two phases, which form the basis
of our approach:

  * •

Learning. Extracting zone sequences from historical data enables the counting
and collection of zone transitions in a count matrix.

  * •

Prediction. For new routes, a combination of a distance and count matrix is
used in a TSP formulation which solution provides the zone sequence; after
which within zones the stop sequences are iteratively determined by means of a
series of OTSPs.

###  3.1 Learning Zone Preferences

For learning, we exploit the information contained in route data to extract
information about drivers’ preferences at the level of zones. Each route
consists of an ordered series of $n+1$ stops, namely $n$ delivery nodes (this
number differs per route) and the station. Therefore, a route is represented
as a sequence of tuples,
$\texttt{route}=((x_{0},y_{0},z_{0}),\ldots,(x_{n},y_{n},z_{n}))$, with
$(x_{i},y_{i})\in\mathbb{R}^{2}$ indicating the stop’s geographical location,
and $z_{i}$ the corresponding zone id, with
$z_{i}\in\left\\{Z_{1},\ldots,Z_{M}\right\\}$ for $i\in\\{0,\ldots,n\\}$. As
the data provides the stop sequence per route, we have to convert it to the
higher-level sequence of zones.

Given the observations in the introduction and a navigational perspective that
all stops within a zone are completed before a driver moves to another zone
(Section 2.3), we determine for each route the sequence of unique zones. To
illustrate the procedure let us consider a route with $n$ delivery nodes.
First, we consider the sequence of zones for each stop $(z_{1},\dots,z_{n})$.
Then, the route’s zone sequence is reduced to the sequence of pairs
$((\zeta_{1},\tau_{1})\ldots,(\zeta_{\ell},\tau_{\ell}))$, where each pair
($\zeta_{i},\tau_{i}$) represents how many nodes ($\tau_{i}$) belonging to the
same zone ($\zeta_{i}$) are visited in a row. When the $\zeta_{i}$ are
pairwise distinct, we consider $(\zeta_{1},\dots,\zeta_{\ell})$ to be the
desired sequence of unique zones.

###### Example 3.1

The following route consisting of six stops over three zones is reduced into a
sequence of three unique zones:

$(z_{1},z_{2},z_{3},z_{4},z_{5},z_{6})=(Z_{1},Z_{1},Z_{1},Z_{2},Z_{3},Z_{3})\rightarrow((\zeta_{1},3),(\zeta_{2},1),(\zeta_{3},2))\rightarrow(\zeta_{1},\zeta_{2},\zeta_{3})=(Z_{1},Z_{2},Z_{3}).$

In case of a zone’s reoccurrence ($\zeta_{i}=\zeta_{j}$ for an $i<j$), we only
keep the occurrence with the highest number of stops by comparing $\tau_{i}$
to $\tau_{j}$, or keep the earliest appearance in case of a tie, here that is
$\zeta_{i}$. This procedure is repeated until we are left with a sequence of
singly occurring zones. Two extreme, yet illustrative, examples of this
procedure are given below.

###### Example 3.2

The reduction of two routes consisting of six stops over three zones under the
procedure outlined above is illustrated below. First, zone $Z_{3}$ is put last
as $\tau_{1}=1$ and $\tau_{4}=2$:

$\displaystyle(Z_{3},Z_{1},Z_{1},Z_{2},Z_{3},Z_{3})$
$\displaystyle\rightarrow((\zeta_{1},1),(\zeta_{2},2),(\zeta_{3},1),(\zeta_{4},2))\rightarrow(\zeta_{2},\zeta_{3},\zeta_{4})=(Z_{1},Z_{2},Z_{3}).$

In the next case, zone $Z_{1}$ occurs thrice, but breaking the tie twice puts
zone $Z_{1}$ at its first occurrence in the sequence:

$\displaystyle(Z_{1},Z_{3},Z_{1},Z_{2},Z_{2},Z_{1})$
$\displaystyle\rightarrow((\zeta_{1},1),(\zeta_{2},1),(\zeta_{3},1),(\zeta_{4},2),(\zeta_{5},1))\rightarrow(\zeta_{1},\zeta_{2},\zeta_{4})=(Z_{1},Z_{3},Z_{2}).$

Evidently, not all zones are visited within each route, therefore we keep
track of all $M$ zones that have been visited at least once in any route.
Having reduced each route to its corresponding zone sequence, we compute an
asymmetric count matrix $N$ where each entry $N_{ij}$ represents the number of
times a driver went from the $i$-th to the $j$-th zone, so that $i$ and $j$
range from $1$ to $M+S$, where $S$ is the number of stations in the dataset.
The pseudocode related to the procedures outlined is provided in the Appendix.

###  3.2 Predicting the Zone Sequence

In the next phase, we predict the route that a driver will take. For that
purpose, we translate the count matrix, created in the previous phase, to a
transition matrix $P$ via $P_{ij}=\dfrac{N_{ij}}{\sum_{j=1}^{m}N_{ij}}$, so
that $P_{ij}$ reflects the Markovian probability of transitioning from zone
$i$ to zone $j$.

Besides the zone sequences, the locations of the station and zones are
necessary ingredients in the suspected trade-off a driver makes. For this
purpose, we define the center of a zone, if not given, by computing the mean
latitude and longitude of all stops that were visited in that zone. Having the
zone centers’ location, we can compute the distances $\delta_{ij}$ between
zone $i$ and zone $j$. We normalize the distances by dividing each element
through the largest, i.e., $\max_{i,j}\delta_{ij}$, and store these distances
in a matrix $D\in\mathbb{R}^{(M{+S})\times(M{+S})}$.

As an alternative to the Euclidean distances in $D$, a sensible choice is to
use the travel times (or another measure for costs between stops) as a link to
measure ‘distance’ between zones. However, in the case when we are only given
route data, we can proceed with a two-step procedure to have an approximation
of the travel times between zones. Since stops are typically assigned to
zones, we first identify for each zone the closest stop within the route to
the zone center by using the Euclidean distance metric. Second, once each zone
has a stop assigned, we extract the travel times between these stops to
compose a matrix $T$ which consists of elements $T_{ij}$ corresponding to
normalized travel times from station/zone $i$ to $j$, i.e.,
$T_{ij}=\frac{t_{ij}}{\max_{i,j}t_{ij}}$ where $t_{ij}$ are the provided
travel times. We will compare the performances of both options.

Since we are to minimize a cost matrix, we reverse the transition matrix to
$1-P_{ij}$ in the cost matrix to arrive at the following linear combination
with $\omega\in[0,1]$ being a weight parameter:

$C_{ij}=\omega D_{ij}+(1-\omega)(1-P_{ij}),$ (1)

where $D_{ij}$ can be interchanged with $T_{ij}$, and additionally we set the
diagonal entries $C_{ii}=0$. In fact this problem can be considered as a TSP
with unknown costs as we do not a priori know how distance is evaluated
against likelihood of zone-to-zone transitions, which contain the implicit
tacit knowledge and drivers’ preferences. So a good value of $\omega$ is still
to be determined. Next, we are able to use this cost matrix in a TSP
formulation, which optimal solution constitutes a tour through all zones in
which one or more stops are located. For details about the implementation we
refer the reader to the Appendix.

###  3.3 Predicting the Stop Sequence

At this point, we have a sequence of zones that yet has to be transformed into
a sequence of stops, e.g., we are in Figure 3(a) where the zone sequence is
$(Z_{A},Z_{B},Z_{C}).$ As the problem size in each zone is considerably
smaller, drivers are able to make (near-)optimal navigational decisions that
minimize travel time or distance, see Section 2.3. To mimic that behavior, we
employ again the framework for solving TSPs, but with the additional
flexibility to be an open TSP (OTSP) as we do not need to close the loop
within a zone, but rather traverse from one zone to another. Although there
are several options available, the procedure detailed in Figure 3 is in
accordance with practice. The approach builds on two elements, which are
applied iteratively in each zone:

  * •

The last stop of a previous zone (or station, Figure 3(b)) is used as the
starting point of the OTSP for the current zone.

  * •

An additional stop of the next zone is added to each OTSP to add sense of
direction to find a good stop to finish its zone.

As seen in Figure 3(a), the extra stop (not filled) is selected as the one
closest (in terms of Euclidean distance) to the computed zone centers
$(\times)$. Naturally, after the solution for an OTSP is found with the
additional stop added, the last stop is removed from the solution such that
for the subsequent zone the starting point is again a stop in its previous
zone. The procedure continues until the last zone is reached for which the
station is used as the additional ‘stop’.

(a) Overview of stops with zone centers $(\times)$.

(b) Solution of the OTSP for zone A.

(c) Solution of the OTSP for zone B.

(d) Solution of the OTSP for zone C.

Figure 3: An example of the iterative procedure that is applied on (a), which
each time ((b), (c) & (d)) solves an open traveling salesman problem locally.

Because of clarity of the exposition, we considered Euclidean distances in the
figure, but the implementation can also be fed with travel times. The
resulting solution thus forms a logical sequence from a navigational point of
view as it avoids path crossings and mimics human behavior when solving
(O)TSPs, see Section 2.3. A full implementation of the approach is given in
the Appendix and is tailored such that it fits the dataset corresponding to
the Amazon Last Mile Research Challenge. Also, we use this dataset to show the
adequacy of the approach, as the achieved performance puts the approach among
the top-tier of submissions for this challenge.

##  4 Performance Evaluation

To assess the performance of the proposed methodology we employ it on a
routing dataset, which was part of the Amazon Last Mile Research Challenge
(Amazon and MIT Center for Transportation & Logistics 2021a). Also the metrics
that we use are adopted from this challenge and are provided in Section 4.1.

The dataset counts 6,112 routes from $17$ different stations. We split the
data such that 5,112 routes are used for learning and $1,000$ routes are
reserved as an out-of-sample test set. Furthermore, we impute missing zone ids
by taking over the id of the closest stop in terms of travel time—a sensible
choice from the driver’s perspective. The experiments related to our
methodology are run using an Intel Core i7-10850H 2.70 GHz CPU and 32 GB RAM.
The learning takes approximately 1 minute and with a test set of $1,000$
routes, it takes on average 30 minutes to render all routes. However, for
choosing a right value of $\omega$ for Eq. (1) and reconsidering it in 4.3, we
rely on 5-fold cross validation over the training set of $5,112$ routes.

###  4.1 Route Prediction Metrics

The route prediction performance can be evaluated in two dimensions, that is,
how often and by how much a prediction $B$ deviates from the realized stop
sequence $A:=(0,1,\ldots,n)$, and thus $B$ is in essence a permutation of $A$.
Below, we describe these different components and how they together form an
overall route performance metric.

  * •

The sequence deviation SD${}_{\text{stop}}(A,B)$ measures the difference in
stop sequences $A$ and $B.$ Given that in both cases all $n+1$ stops have been
visited, we take the stop sequence of $B$, and for each position we trace back
when it has been visited in $A$ to create a vector
$(a_{0},a_{1}\ldots,a_{n})=\pi(A)$ so that

$\textbf{SD${}_{\text{stop}}(A,B)$}:=\frac{2}{n(n-1)}\sum_{i=1}^{n}\left(|a_{i}-a_{i-1}|-1\right).$

Note that indeed if $A\equiv\pi(A)=B$, the deviation becomes $0$. Lastly, due
to the fact that we deduce the zone sequence from a route by running the
procedure outlined in Section 3.1, we can also compute the zone sequence
deviation SD${}_{\text{zone}}(A,B)$ that on the higher level compares sequence
deviations between two zone sequences $A$ and $B$.

  * •

Besides considering the position in the sequence, one can also count the
number of edit operations needed to come from the predicted route to the
realized route, familiarly known as the Levenshtein distance as introduced by
Levenshtein et al. (1966). This distance metric, renamed to
ERP${}_{\text{edit}}(A,B)$, counts the number of deletions and insertions to
get the same sequence.

  * •

Apart from considering the concurrence and concordance of stop sequences, the
size of deviation between two sequences $A$ and $B$ should be evaluated as
well. To this end, normalized travel time is introduced to evaluate the
difference between two routes on a stop basis; let $A[i]$ and $B[j]$ be the
$i$-th stop of $A$ and the $j$-th stop of $B$, with $i,j=0,1,\dots,n$, then it
is defined as

$\displaystyle\textbf{time${}_{\text{norm}}$}(A[i],B[j])$
$\displaystyle:=Y_{A[i],B[j]}-\min_{i,j}Y_{A[i],B[j]},$ $\displaystyle
Y_{A[i],B[j]}$ $\displaystyle:=\frac{t_{A[i],B[j]}-\bar{t}}{\text{std}(t)},$

where $t_{A[i],B[j]}$ is the travel time from stop $A[i]$ to stop $B[j]$, and
$\bar{t}$ and $\text{std}(t)$ are the mean and standard deviation over all
travel times between the stops in the route. Considering per stop the
difference between two sequences, we compute the so-called Edit Distance with
Real Penalty component (ERP${}_{\textbf{norm}}$). Formally, it is computed as

$\textbf{ERP${}_{\text{norm}}(A,B)$}:=\sum_{i=0}^{n}\textbf{time${}_{\text{norm}}$}(A[i],B[i]).$

As an aside, dividing ERP${}_{\text{norm}}(A,B)$ by ERP${}_{\text{edit}}(A,B)$
translates to the average additional travel time incurred by deviating, and is
called ERP${}_{\text{ratio}}(A,B)$.

Besides these separate metrics, they are combined to create a comprehensive
score

$\textbf{route\\_score}=\frac{\textbf{SD${}_{\text{stop}}(A,B)$}\cdot\textbf{ERP${}_{\text{norm}}(A,B)$}}{\textbf{ERP${}_{\text{edit}}(A,B)$}},$

so that a performance score is obtained by averaging over all $I$ routes to be
predicted:

$\displaystyle\textbf{Performance}=\frac{1}{I}\sum_{i=1}^{I}\textbf{route\\_score}_{i}.$
(2)

Evidently, a lower route score means that the predicted sequence coincides
more with the realized sequence; $0$ means they are exactly the same.

###  4.2 Choosing the Weight

As a first experiment, we consider varying the weight parameter $\omega$ in
Eq. (1). Furthermore, we can replace distance by travel times $T$ as these are
provided in Amazon’s dataset. The performances when varying $\omega$ in both
versions are displayed in Figure 4.

$0$$0.1$$0.2$$0.3$$0.4$$0.5$$0.6$$0.7$$0.8$$0.9$$1$0.030.040.050.06$\omega$Performance$\omega
D+(1-\omega)(1-P)$$\omega T+(1-\omega)(1-P)$ Figure 4: Comparison of
performances, as defined in Eq (2), when varying $\omega$ in Eq. (1), where
Euclidean distance is dashed and (expected) travel time solid. The curves are
obtained averaging over the 5-fold cross-validation scores.

As expected, we find that only relying on historical information ($\omega=0$),
or only using distance or travel times ($\omega=1$) does not perform as well
as having a trade-off between the two. Overall, we observe that an $\omega$
value around $0.9$ achieves the best trade-off, but actually any value between
0.1 and 0.9 in Eq. (1) (dashed and solid lines) will result in a performance
of around $0.035.$ This observation demonstrates the robustness of the model
as being insensitive to misspecification of $\omega$. In addition, we conclude
that the difference between travel times or Euclidean distance is small and
thus distance can be used as a good proxy for travel time.

###  4.3 Adjusting the Station’s Weights

Reconsidering Eq. (1), wherein we define an $\Omega$ matrix to weigh
historical information against travel times, we further study the role of the
station as the starting and ending point of a route. For the sake of
readability, we consider a matrix $\Omega^{(r)}$ which contains the weights in
$\Omega$ corresponding to the station and zones that are visited in a specific
route $r$ where the station is related to the first row and column. So as to
refine the role of the station, we define $\Omega^{(r)}$ matrix as

$\Omega^{(r)}=\left(\begin{array}[]{c|cccc}0&\omega_{F}&\dots&\omega_{F}\\\
\hline\cr\omega_{L}&\omega_{Z}&\dots&\omega_{Z}\\\
\vdots&\vdots&\ddots&\vdots\\\
\omega_{L}&\omega_{Z}&\dots&\omega_{Z}\end{array}\right),$

where $\omega_{F}\in[0,1]$ is the weight in the cost function from the station
to the first zone. Analogously, $\omega_{Z}\in[0,1]$ sets the balance for
zone-to-zone transitions, and lastly $\omega_{L}\in[0,1]$ is the weight
corresponding to the return to the station. The elements $\Omega_{ij}$ from
this matrix are put in the following refinement of Eq. (1) with travel times:

$C_{ij}=\Omega_{ij}\,T_{ij}+(1-\Omega_{ij})\,(1-P_{ij}),$ (3)

where, similarly to Eq. (1), the diagonal values are set to zero. With this
generalized cost matrix we first set $\omega_{L}=1$ as we anticipate that at
the end of a route heading back to the station should be the only concern, and
thus the number of times that a specific zone acted as the finish point of a
route can be disregarded. This resonates with the thought that a driver is
only focused at the end on the travel time it takes to return to station.
Having this parameter restrained, we study combinations of $\omega_{F}$ and
$\omega_{Z}$ to obtain the contour plot in Figure 5. In the plot, we observe
that the optimal point should lie somewhere in the region where
$\omega_{F}\in[0.1,0.2]$ and $\omega_{Z}\in[0.7,0.8]$. So when starting a
route in the first transition from the station, historical information should
be weighted more important than regular zone-to-zone transitions. This
implicates that when determining the route globally preference and
navigational considerations have a greater impact than the standard trade-off
with travel times.

Figure 5: A contour plot showing the performance across different weightings
of $\omega_{F}$ (station-to-zone) and $\omega_{Z}$ (zone-to-zone) transitions,
while keeping $\omega_{L}=1$ (zone-to-station) in Eq. (3). The values are
obtained averaging over the 5-fold cross-validation scores.

The choice of $\omega_{L}=1$ is justified by the degradation in performance
for values less than 1. This is shown in Figure 6 by taking four combinations
of $(\omega_{F},\omega_{Z})$ that lie around the inner contour—close to the
optimum. In all, when comparing the points in detail, we conclude that the
best performance in our cross-validation is obtained when
$(\omega_{F},\omega_{Z},\omega_{L})=(0.2,0.8,1)$, which results in the
reported performance of the final model in Figure 7.

$0$$0.1$$0.2$$0.3$$0.4$$0.5$$0.6$$0.7$$0.8$$0.9$$1$0.0300.0350.040$\omega_{L}$Performance$\omega_{F}=0.1,\,\omega_{Z}=0.7$$\omega_{F}=0.1,\,\omega_{Z}=0.8$$\omega_{F}=0.2,\,\omega_{Z}=0.7$$\omega_{F}=0.2,\,\omega_{Z}=0.8$
Figure 6: A sensitivity study on $\omega_{L}$, (zone-to-station), for selected
combinations of $\omega_{F}$ and $\omega_{Z}$ that lie around the inner
contour of Figure 5. The performances, as defined in Eq. (2), are obtained by
averaging over the 5-fold cross-validation scores.

###  4.4 Comparison against Benchmarks

In Table 4.4, we provide the performance of our approach abbreviated to LG-OL
and compare it to two benchmarks to shed some light onto the range of
performance scores that can be attained. Besides the metrics outlined in
Section 4 we also added the mean travel time in seconds in the Time column.
First, we apply the Nearest Neighbor algorithm (Figure 2(a)) on the out-of-
sample test set of $1,000$ routes. This algorithm sequentially adds a
connection from the last stop to the nearest unvisited stop; a strategy that
is sometimes considered to model human route planning behavior (Graham, Joshi,
and Pizlo 2000, Wiener, Ehbauer, and Mallot 2009). As another idea, consider
the route which minimizes the travel times as the route a driver follows
(Figure 2(b)). Such route is generally hard to compute, but within the range
of state-of-the-art solvers (e.g., Gurobi). Surprising is that drivers average
travel time (9496 seconds) is about 18% more than when travel time
minimization was the only objective (Full TSP); indeed showing that drivers’
navigational decisions cannot be approached as a standard optimization
problem.

Our methodology obtains a prediction performance (Performance) of $0.0299$—a
more than $64\%$ reduction compared to the aforementioned benchmarks.
Interestingly, the Nearest Neighbor and full TSP models have similar
ERP${}_{\text{ratio}}$ scores, but the full TSP renders better sequences
(lower SD${}_{\text{zone}}$ and SD${}_{\text{stop}}$ scores). When considering
our individual route scores of the test set in Figure 7 against the
benchmarks, we find that our model succeeds to have $40\%$ of the test routes
to be near-perfectly predicted, as they are below 0.01; and even gets $80\%$
of them below the $0.05$ threshold. In retrospect, reconsidering Figure 4 we
find that even for any $\omega$ the unadjusted approach outperforms these
benchmarks, supporting the fundamental idea of our learn global, optimize
local approach.

Evaluating our performance to the submissions for the challenge (Amazon and
MIT Center for Transportation & Logistics 2021b), we find that the approach
propels itself in the top-tier of submissions, which performances range from
0.025 to 0.050. Finally, to verify the correctness and potential of the
approach, we did an additional hypothetical experiment on the test set of
which the results are provided in Table 4.4 as LG-OL (hypothetical). We
extracted the zone sequences using Section 3.1 on the test set so as to
consider the situation that we are fully capable to predict the zone sequence
right. Thus remaining are the inconsistencies between the design choices made
in Section 3, and drivers’ behavior that is not reflected in the methodology.
In the experiment, we find that the performance becomes 0.0094, which is well
below all submissions for this challenge.

Comparison of the methodology against benchmarks on the out-of-sample test
set.  Model SD${}_{\text{zone}}$ SD${}_{\text{stop}}$ ERP${}_{\text{ratio}}$
Time (s) Performance Nearest Neighbor 0.2441 0.0873 1.3382 9865 0.1119 Full
TSP 0.1912 0.0650 1.3330 8036 0.0826 LG-OL 0.0940 0.0335 0.7389 8994 0.0299
LG-OL (hypothetical) 0 0.0238 0.3831 9535 0.0094 Driver 0 0 0 9496 0

Nearest NeighborFull TSPLG-OL0.000.050.100.150.200.25Route Score Figure 7:
Boxplot representing the route score distribution on the out-of-sample test
set of 1,000 routes for our model, and the two benchmarks of Table 4.4.

##  5 Conclusion and Discussion

The costs associated with last-mile delivery surmount any other shipping
costs. Hence, good predictions about the route taken are crucial from an
operations management point of view. The proposed approach leads to several
broadly applicable insights about delivery driving.

The approach first determines a zone sequence by solving a traveling salesman
problem (TSP) over a cost matrix which is a weighted combination of historical
information, which has to be learnt from realized routes, and distance (or
travel time) between zones. Next, on the local level, adhering to the
established zone sequence on the global level, the order of stops to be
visited in a zone is found by solving an open TSP (OTSP), which generates
locally optimal solutions. Repeating this for all subsequent zones and
patching them together generates a feasible route. Therefore this approach can
be considered to learn global (over zones) and optimize local (within zones).
Besides that the approach can be argued to follow the logic of human
navigation, the approach is convenient, because it only requires learning on
the zone level as by tracking the zone transitions. The breakdown makes the
approach relatively fast, as the zone-sequence TSP and within-zone OTSPs are
not so computationally involved compared to solving a single TSP for the
entire route. This also renders the approach amenable to real-time training
using newly realized routes of the same day.

The key element in the approach is the computation of a cost matrix that
combines distance and historical information when establishing the zone
sequence as a TSP. The experiments show that a weighted combination with
travel times provides a slight performance improvement over using Euclidean
distance. This fact underpins that when determining the zone sequence an
implicit trade-off is made between the driver’s zone-to-zone preferences and a
measure of distance. Studying the weight value by changing it from $0$ to $1,$
where the setting of 1 corresponds to only considering travel time, uncovers
that the exact value is not critical. Overall, the best performance is
obtained when the value is close to $0.9$. However, the weighting of station-
to-zone and zone-to-station transitions can be further adapted in the cost
matrix to better mimic driver’s behavior resulting in another performance
improvement. In fact, these adaptations are motivated by the logic that at the
beginning of a route the driver is more concerned with where to start, whereas
at the end, when all stops are visited, the only consideration is the return
to the station. Indeed, the experiments confirm these hypotheses and we find
that the travel time matrix should be weighted less important in the station-
to-zone transitions than in the zone-to-zone transitions, and from zone-to-
station it is the only factor. These insights come with the caveat that
drivers were suggested a route beforehand (see Amazon and MIT Center for
Transportation & Logistics 2021a), which might have tainted the driver’s
decisions.

The approach as-is only utilizes historical stop data (with zone ids) in
making predictions. Using additional data which are part of the case study,
such as package size, start time, quality of the route, did not provide any
clear starting points for improvement. Since the zone ids were given, it is
worthwhile to further study the underlying segmentation, because it determines
the breakdown of the problem into a global and a series of local problems and
as such it is pivotal to the performance of the approach. We have seen that
occurrence of zone revisitations are rare. However, in 2.6% of the routes
there are more than five revisitations, which hints that in these routes it
might be valuable to reconsider the zone segmentation.

The test set also provides us an opportunity to quantify the remaining
learning gap. By first extracting the zone sequences in the test set after
which we apply the approach, we arrive at a score just below $0.01.$ This
shows that the learn global and optimize local approach can yield a score as
low as 0.009 (as a reference we get it to 0.0299, as shown in Section 4). The
learning took place on only 5,112 routes, which results in a sparse transition
matrix as the zone transition probabilities, when available, are learnt based
on the relatively few routes available for each region. Increasing the number
of routes for learning, for example by daily tracking of realized routes and
updating the count matrix accordingly, the approach can likely improve its
accuracy on zone sequence prediction, which have a profound impact on the
overall performance.

Another endeavor is to further study the cost matrix structure. The
methodology readily incorporates other functional forms, and also different
weightings between historic information and distance-based measures. These
weightings might depend on the size of the dataset available and even on the
length of the route itself. Nevertheless, we find that a linear structure
between distance or travel time and historical information is robust to
misspecification of this weight parameter and is practically sound.

Finally, the methodology can be employed in other (last-mile) delivery
concepts, see for example Boysen, Fedtke, and Schwerdfeger (2021), or
generalized to other operational optimization problems, which are susceptible
to human interference in the decision making process—a logical starting point
is the allied Vehicle Routing Problem, e.g., Bräysy and Gendreau (2005a). In
all, this research demonstrates that integrating learning in an optimization
pipeline leverages its connection and value to practice.

## References

  * Amazon and MIT Center for Transportation & Logistics (2021a) Amazon and MIT Center for Transportation & Logistics, 2021a _Amazon Last-Mile Routing Research Challenge_. URL https://routingchallenge.mit.edu/about-the-challenge/, (accessed on 30 April 2022).
  * Amazon and MIT Center for Transportation & Logistics (2021b) Amazon and MIT Center for Transportation & Logistics, 2021b _Amazon Last-Mile Routing Research Challenge_. URL https://routingchallenge.mit.edu/last-mile-routing-challenge-team-performance-and-leaderboard/, (accessed on 30 April 2022).
  * Applegate et al. (2011) Applegate DL, Bixby RE, Chvátal V, Cook WJ, 2011 _The Traveling Salesman Problem_ (Princeton University Press).
  * Boyer, Prud’homme, and Chung (2009) Boyer KK, Prud’homme AM, Chung W, 2009 _The last mile challenge: Evaluating the effects of customer density and delivery window patterns_. _Journal of Business Logistics_ 30(1):185–201.
  * Boysen, Fedtke, and Schwerdfeger (2021) Boysen N, Fedtke S, Schwerdfeger S, 2021 _Last-mile delivery concepts: a survey from an operational research perspective_. _OR Spectrum_ 43(1):1–58.
  * Bräysy and Gendreau (2005a) Bräysy O, Gendreau M, 2005a _Vehicle routing problem with time windows, part I: Route construction and local search algorithms_. _Transportation Science_ 39(1):104–118.
  * Bräysy and Gendreau (2005b) Bräysy O, Gendreau M, 2005b _Vehicle routing problem with time windows, part II: Metaheuristics_. _Transportation Science_ 39(1):119–139.
  * Campuzano, Obreque, and Aguayo (2020) Campuzano G, Obreque C, Aguayo MM, 2020 _Accelerating the Miller-Tucker-Zemlin model for the asymmetric traveling salesman problem_. _Expert Systems with Applications_ 148:113229.
  * Canoy and Guns (2019) Canoy R, Guns T, 2019 _Vehicle routing by learning from historical solutions_. _International Conference on Principles and Practice of Constraint Programming_ , 54–70 (Springer).
  * Chopra (2003) Chopra S, 2003 _Designing the distribution network in a supply chain_. _Transportation Research Part E: Logistics and Transportation Review_ 39(2):123–140.
  * Current and Schilling (1989) Current JR, Schilling DA, 1989 _The covering salesman problem_. _Transportation Science_ 23(3):208–213.
  * Dantzig and Ramser (1959) Dantzig GB, Ramser JH, 1959 _The truck dispatching problem_. _Management Science_ 6(1):80–91.
  * Deloison et al. (2020) Deloison T, Hannon E, Huber A, Heid B, Klink C, Sahay R, Wolff C, 2020 _The future of the last-mile ecosystem: Transition roadmaps for public-and private-sector players_ (World Economic Forum).
  * Dry et al. (2006) Dry M, Lee MD, Vickers D, Hughes P, 2006 _Human performance on visually presented traveling salesperson problems with varying numbers of nodes_. _The Journal of Problem Solving_ 1(1):4.
  * eMarketer (2020) eMarketer, 2020 _US ecommerce 2020, coronavirus boosts ecommerce forecast and will accelerate channel-shift_. URL https://www.emarketer.com/content/us-ecommerce-2020, (accessed on 30 April 2022).
  * Gendreau et al. (1998) Gendreau M, Hertz A, Laporte G, Stan M, 1998 _A generalized insertion heuristic for the traveling salesman problem with time windows_. _Operations Research_ 46(3):330–335.
  * Gendreau, Laporte, and Vigo (1999) Gendreau M, Laporte G, Vigo D, 1999 _Heuristics for the traveling salesman problem with pickup and delivery_. _Computers & Operations Research_ 26(7):699–714.
  * Gevaers, Van de Voorde, and Vanelslander (2011) Gevaers R, Van de Voorde E, Vanelslander T, 2011 _Characteristics and typology of last-mile logistics from an innovation perspective in an urban context_. _City Distribution and Urban Freight Transport_ (Edward Elgar Publishing).
  * Gevaers, Van de Voorde, and Vanelslander (2014) Gevaers R, Van de Voorde E, Vanelslander T, 2014 _Cost modelling and simulation of last-mile characteristics in an innovative b2c supply chain environment with implications on urban areas and cities_. _Procedia - Social and Behavioral Sciences_ 125:398–411, eighth International Conference on City Logistics 17-19 June 2013, Bali, Indonesia.
  * Goodman (2005) Goodman RW, 2005 _Whatever you call it, just don’t think of last-mile logistics, last_. _Global Logistics & Supply Chain Strategies_ 9(12).
  * Graham, Joshi, and Pizlo (2000) Graham SM, Joshi A, Pizlo Z, 2000 _The traveling salesman problem: A hierarchical model_. _Memory & Cognition_ 28(7):1191–1204.
  * Gunasekaran, Patel, and McGaughey (2004) Gunasekaran A, Patel C, McGaughey RE, 2004 _A framework for supply chain performance measurement_. _International Journal of Production Economics_ 87(3):333–347.
  * Jiang et al. (2014) Jiang J, Gao J, Li G, Wu C, Pei Z, 2014 _Hierarchical solving method for large scale tsp problems_. _International symposium on neural networks_ , 252–261 (Springer).
  * Karp (1977) Karp RM, 1977 _Probabilistic analysis of partitioning algorithms for the traveling-salesman problem in the plane_. _Mathematics of Operations Research_ 2(3):209–224.
  * Kool, van Hoof, and Welling (2019) Kool W, van Hoof H, Welling M, 2019 _Attention, learn to solve routing problems!_ _International Conference on Learning Representations_ , URL https://openreview.net/forum?id=ByxBFsRqYm.
  * Krumm (2008) Krumm J, 2008 _A Markov model for driver turn prediction_. _Society of Automotive Engineers (SAE) World Congress, April 2008_.
  * Lawler (1985) Lawler EL, 1985 _The Traveling Salesman Problem: A Guided Tour of Combinatorial Optimization_ (Wiley-Interscience Series in Discrete Mathematics).
  * Levenshtein et al. (1966) Levenshtein VI, et al., 1966 _Binary codes capable of correcting deletions, insertions, and reversals_ 10:707–710.
  * Liao and Liu (2018) Liao E, Liu C, 2018 _A hierarchical algorithm based on density peaks clustering and ant colony optimization for traveling salesman problem_. _IEEE Access_ 6:38921–38933.
  * Lin et al. (2014) Lin C, Choy K, Ho G, Chung SH, Lam H, 2014 _Survey of green vehicle routing problem: Past and future trends_. _Expert Systems with Applications_ 41:1118–1138.
  * MacGregor and Ormerod (1996) MacGregor JN, Ormerod T, 1996 _Human performance on the traveling salesman problem_. _Perception & Psychophysics_ 58(4):527–539.
  * Miller, Tucker, and Zemlin (1960) Miller CE, Tucker AW, Zemlin RA, 1960 _Integer programming formulation of traveling salesman problems_. _Journal of the ACM_ 7(4):326–329.
  * Montoya-Torres et al. (2015) Montoya-Torres JR, López Franco J, Nieto Isaza S, Felizzola Jiménez H, Herazo-Padilla N, 2015 _A literature review on the vehicle routing problem with multiple depots_. _Computers & Industrial Engineering_ 79:115–129.
  * OpenStreetMap (2021) OpenStreetMap, 2021 _Planet dump retrieved from https://planet.osm.org_. URL https://www.openstreetmap.org, (accessed on 2021-20-11).
  * Orman and Williams (2007) Orman A, Williams HP, 2007 _A survey of different integer programming formulations of the travelling salesman problem_. _Advances in Computational Management Science_ , volume 9, 91–104 (Springer Berlin Heidelberg).
  * Papadimitriou (1977) Papadimitriou CH, 1977 _The Euclidean traveling salesman problem is NP-complete_. _Theoretical Computer Science_ 4(3):237–244.
  * Pizlo et al. (2006) Pizlo Z, Stefanov E, Saalweachter J, Li Z, Haxhimusa Y, Kropatsch WG, 2006 _Traveling salesman problem: A foveating pyramid model_. _The Journal of Problem Solving_ 1(1):8.
  * Purkayastha et al. (2020) Purkayastha R, Chakraborty T, Saha A, Mukhopadhyay D, 2020 _Study and analysis of various heuristic algorithms for solving travelling salesman problem–a survey_. Mandal JK, Mukhopadhyay S, eds., _Proceedings of the Global AI Congress 2019_ , 61–70 (Springer Singapore).
  * Sengupta, Mariescu-Istodor, and Fränti (2018) Sengupta L, Mariescu-Istodor R, Fränti P, 2018 _Planning your route: Where to start?_ _Computational Brain & Behavior_ 1(3-4):252–265.
  * Speranza and Archetti (2014) Speranza M, Archetti C, 2014 _A survey on matheuristics for routing problems_. _EURO Journal on Computational Optimization_ 2:223–246.
  * Taniguchi and Thompson (2002) Taniguchi E, Thompson RG, 2002 _Modeling city logistics_. _Transportation Research Record_ 1790(1):45–51.
  * United Nations (2018) United Nations, 2018 _68% of the world population projected to live in urban areas by 2050, says UN._ URL https://www.un.org/development/desa/en/news/population/2018-revision-of-world-urbanization-prospects.html, (accessed on 2021-20-11).
  * Van Rooij, Stege, and Schactman (2003) Van Rooij I, Stege U, Schactman A, 2003 _Convex hull and tour crossings in the Euclidean traveling salesperson problem: Implications for human performance studies_. _Memory & Cognition_ 31(2):215–220.
  * Vickers et al. (2003) Vickers D, Lee MD, Dry M, Hughes P, 2003 _The roles of the convex hull and the number of potential intersections in performance on visually presented traveling salesperson problems_. _Memory & Cognition_ 31(7):1094–1104.
  * Wang et al. (2015) Wang X, Ma Y, Di J, Murphey YL, Qiu S, Kristinsson J, Meyer J, Tseng F, Feldkamp T, 2015 _Building efficient probability transition matrix using machine learning from big data for personalized route prediction_. _Procedia Computer Science_ 53:284–291.
  * Wiener, Ehbauer, and Mallot (2009) Wiener J, Ehbauer N, Mallot H, 2009 _Planning paths to multiple targets: Memory involvement and planning heuristics in spatial problem solving_. _Psychological Research PRPF_ 73(5):644–658.
  * Wiener and Mallot (2003) Wiener JM, Mallot HA, 2003 _’Fine-to-coarse’ route planning and navigation in regionalized environments_. _Spatial cognition and computation_ 3(4):331–358.
  * Ye et al. (2015) Ye N, Wang Zq, Malekian R, Lin Q, Wang Rc, 2015 _A method for driving route predictions based on hidden Markov model_. _Mathematical Problems in Engineering_ 2015\.

Pseudocode

In this Appendix, we provide the pseudocode to the methodology presented in
Section 3. Although the pseudocode has been designed on the basis of the
Amazon Last-Mile Routing Research Challenge dataset, it can be easily adapted
to different data structures. As a side note, the zone ids are not consistent
for different stations, therefore, we run our code for each $s\in\mathcal{S}$,
with $\mathcal{S}$ the set of stations. The pseudocode used for the learning
phase is reported in Algorithm 1, where the function ToZoneSeq converts the
realized sequence of stops into a sequence of zones (see Section 3.1), while
function ComputeCountMatrix computes an asymmetric matrix where the
$({i,j})$-th entry represents the number of times a driver went from the
$i$-th to the $j$-th zone. Algorithm 2 illustrates the pseudocode used in the
prediction phase. Here, in order to compute the cost matrix, as thoroughly
described in Section 3.2, we compute the distance matrix where each entry
refers to the travel time between two different zone centers (computed with
EstimateZoneCenter).
The code used to run the experiments presented in Section 4 is open source and
will be made available.

Algorithm 1 Learning Zone Preferences and Preprocessing

1:$\mathcal{D}_{s}$: route data set corresponding to station $s$

2:StopSeqs: actual sequence of stops for each route in $\mathcal{D}_{s}$

3:$\texttt{CountMatrix}\in\mathbb{R}^{(M+1)\times(M+1)}$: counts of
station/zone transitions

4:ZoneSequences $\leftarrow\emptyset$

5:AllZones $\leftarrow\\{s\\}$

6:for route in $\mathcal{D}_{s}$ do

7: ZoneSequences[route] $\leftarrow$ ToZoneSeq(StopSeqs[route])

8: AllZones
$\leftarrow\texttt{AllZones}\cup\\{\texttt{ZoneSequences[route]}\\}$

9:end for

10:CountMatrix $\leftarrow$ ComputeCountMatrix(AllZones, ZoneSequences)

Algorithm 2 Predicting Routes

1:CountMatrix: output of Algorithm 1

2:Zones, Stops: zones and stops characterizing the route

3:station, TravelTimes

4:PredictedStops: predicted sequence of stops

5:ZoneCenters $\leftarrow$ EstimateZoneCenter(Zones, Stops)

6:DistanceMatrix $\leftarrow$ ComputeDistMatrix(Zones, ZoneCenters)

7:CostMatrix $\leftarrow$ ComputeCostMatrix(DistanceMatrix, CountMatrix)

8:PredictedZoneSeq $\leftarrow$ TSP(station, Zones, CostMatrix)

9:PredictedStopSeq $\leftarrow$ [station]

10:for zone in PredictedZoneSeq do

11: PrevStop $\leftarrow$ PredictedStopSeq[-1] $\triangleright$ Last element
of PredictedStopSeq

12: AddStop $\leftarrow$ ClosestToNextZoneCenter(zone, Stops)

13: PredictedStopSeq $\leftarrow$ PredictedStopSeq \+ [OTSP(PrevStop,
Stops[zone], AddStop, TravelTimes)]

14:end for


# Atomic-scale description
of the paramagnetic susceptibility of non-magnetic Ba2YMoO6

R. J. Radwanski  Center of Solid State Physics, Snt Filip 5, 31-150 Krakow,
Poland. email<EMAIL_ADDRESS>
###### Abstract

We succeeded in remarkably nice theoretical description of the temperature
dependence of the paramagnetic susceptibility of Ba2YMoO6 in the whole
temperature range as originating from the atomic-scale discrete electronic
structure of Mo5+ ions in the 4$d^{1}$ configuration resulting from crystal-
field and spin-orbit interactions. A strong violation of the Curie-Weiss law
and a nonmagnetic state of BYMo experimentally-observed down to 50 mK is
caused by the dramatical reduction of the Mo-ion magnetic moment. The
quantitative reproduction of $\chi$(T) proves that almost all Mo ions are in
the Mo5+ state over the whole temperature range.

##  I Introduction

Ba2YMoO6 (BYMo) is presently widely studied oxide with the 4$d$ open-shell ion
1 ; 2 ; 3 ; 4 ; 5 ; 6 ; 7 ; 8 ; 9 ; 10 ; 11 ; 12 ; 13  exhibiting non-trivial
temperature dependence of the paramagnetic susceptibility $\chi$(T), as is
seen in Fig. 1, with an dramatic decrease with the lowering temperature below
50 K of the $\chi^{-1}$(T) plot 2 ; 5 ; 8 . BYMo does not order magnetically
down to 50 mK 2 ; 5 ; 7  though $\chi$(T) in the 200-300 K region points to
strong antiferromagnetic interactions as one would expect from $\theta_{CW}$
of -(160-200) K, see Fig. 1. 1 ; 2 ; 7 ; 8 . Authors of Ref. 2  explained this
persistent low-temperature non-magnetic state by amorphous arrangements of
spin singlets. Quite recently Romhanyi et al. 11  proposed, in contrast to a
chiral spin-liquid model of Ref. 9 , a spin-orbital model with a dimer-singlet
phase, composed of random arrangement of spin-orbit dimers, without any type
of the long-range order. Such explanations sound much scientifically. They
involve sophisticated mathematics but physical meaning of ”spin singlets” or
”dimer-singlet state” do not provide physical description of the anomalous
behaviour of $\chi$(T) and physical mechanism of so large reduction of the
effective magnetic moment and of the local magnetic moment itself.

In this paper we report the microscopic atomic-scale description of the
temperature dependence of the paramagnetic susceptibility of Ba2YMoO6 as
associated with the Mo5+ ions in the crystal. In our description the on-site
crystal-field (CEF) and spin-orbit (s-o) interactions are fundamentally
important realising a general concept ”From the atomic physics to the solid
state physics” 14 .

Figure 1: Temperature dependence of the paramagnetic susceptibility of Ba2YMoO6 shown in the $\chi^{-1}$ plot where off-Curie-Weiss behaviour is more visible than in the $\chi$(T) plot, after Ref. 2  (much similar data are in Refs 5 ; 8 ). High-temperature dependence of the paramagnetic susceptibility of Ba2YMoO6 yields $\theta_{CW}$ of -160 K what would indicate on strong antiferromagnetic interactions, being in sharp contrast with a nonmagnetic state down to 50 mK, and with a quite substantial effective magnetic moment, about 1.5 $\mu_{B}$, at 300 K. The solid line is the result of this paper. See also Fig. 2.  E [K] |  | +2 | -2 | +1 | -1 | 0 | Sz | Lz | Jz | mz | mx
---|---|---|---|---|---|---|---|---|---|---|---
+22250 | 1 | +0.727$\uparrow$ | +0.686$\uparrow$ |  | +0.021$\downarrow$ |  | +0.4996 | +0.1158 | +0.6154 | +1.1161 | -1.0294
37690 | 10 | +0.686$\downarrow$ | +0.727$\downarrow$ | +0.021$\uparrow$ |  |  | -0.4996 | -0.1158 | -0.6154 | -1.1161 | 1.0294
+21650 | 5 |  |  | +0.037$\downarrow$ |  | 0.999$\uparrow$ | +0.4986 | +0.0014 | +0.5000 | +0.9998 | +1.0899
37090 | 6 |  |  |  | +0.037$\uparrow$ | 0.999$\downarrow$ | -0.4986 | -0.0014 | -0.5000 | -0.9998 | -1.0899
-13460 | 3 | +0.317$\downarrow$ | -0.325$\downarrow$ | +0.891$\uparrow$ |  |  | +0.2942 | +0.7843 | +1.0785 | +1.3733 | -0.7844
1980 | 8 | -0.325$\uparrow$ | +0.317$\uparrow$ |  | +0.891$\downarrow$ |  | -0.2942 | -0.7843 | -1.0785 | -1.3733 | +0.7844
-15000 | 4 |  |  | 0.999$\downarrow$ |  | -0.037$\uparrow$ | -0.4986 | +0.9986 | +0.5000 | +0.0002 | +0.0886
440 | 7 |  |  |  | 0.999$\uparrow$ | -0.037$\downarrow$ | +0.4986 | -0.9986 | -0.5000 | -0.0002 | -0.0886
-15440 | 2 | +0.655$\downarrow$ | -0.605$\downarrow$ | -0.453$\uparrow$ |  |  | -0.2946 | +0.3315 | +0.0369 | -0.2584 | -0.2484
0 | 9 | -0.605$\uparrow$ | +0.655$\uparrow$ |  | -0.453$\downarrow$ |  | +0.2946 | -0.3315 | -0.0369 | +0.2584 | +0.2484
Table 1: Realised spin-orbital states, energies and their eigenfunctions
($|LSL_{z}S_{z}>$ space, with $L$=2 and $S$=1/2) of the Mo5+ ion in Ba2YMoO6
calculated with the octahedral crystal-field parameter $B_{4}$= +300 K
(10$Dq$=3.1 eV), an in-plane tetragonal-distortion parameter $B_{4}^{4}$= +50
K (so $B_{4}^{4}$= +1550 K) and the spin-orbit coupling $\lambda_{s-o}$ =
+1100 K. $L_{z}$ values as +2, -2, +1, -1, and 0 contributions are shown,
whereas $S_{z}$ = +1/2 or -1/2 are shown by up and down arrows, respectively.
In the first column the absolute and the relative energies are given. The
second column is not informative - it shows a number of the row of the output
results with eigenvalues and eigenfunctions from our home-made computer
program. In last five columns the resulting magnetic characteristics from
those eigenfunctions are shown. Note, a rather small $S_{z}$ and quite large
values of $L_{z}$ in some cases and $J_{z}$ values much different from 3/2.

##  II Theoretical outline

We have used the crystal-field approach developed already 50 years ago in a
text-book by Abragam and Bleaney for paramagnetic impurities 15 . Although it
is a text book of 50 years old it only recently has started to be popular for
description of oxides containing $d$ open-shell paramagnetic ions with strong
s-o interactions in iridates and ruthenates 16 . In approach, which we call as
the Quantum Atomistic Solid-State Theory (QUASST), the magnetic and electronic
properties of open-shell oxides are largely determined by single-ion
electronic structure determined by crystal-field and spin-orbit interactions
14 . CEF interactions result from charge interactions of the given
paramagnetic ion in a crystal with its charge surroundings. A benchmark of the
CEF+s-o-based QUASST is the integer valency of the involved transition-metal
ion and the existence of the discrete low-energy electronic structure. As
quite early applicability of such description to rare-earth and actinides
compounds we can mention description of ErNi5 and UPd2Al3 17 . For 3$d$-ion
compounds we can mention description of FeBr2 18  and LaCoO3 19 . As ionic
compound with the 3$d^{1}$ configuration we can mention BaVS3, where we
pointed out importance of the s-o interactions even in 3$d$ oxide manifested
by substantial cancellation of the spin moment by the orbital moment 20 .

##  III Results and discusssion

In this Letter we report on a found atomic-scale description of $\chi$(T) of
BYMo in the whole temperature range, without a need for a $\chi_{o}$
correction, from 0.01 K to 300 K using only three physical parameters;

  1. 1.

octahedral crystal-field parameter $B_{4}^{0}$ = +300 K, demanding $B_{4}^{4}$
= 5$\cdot B_{4}^{0}$ = 1500 K and yielding the largest energy scale of 10$Dq$
= 3.1 eV,

  2. 2.

spin-orbit interactions with $\lambda_{s-o}$ = +1100 K yielding energy scale
of 150 meV, and

  3. 3.

a small off-octahedral distortion, visible here in a change of the $B_{4}^{4}$
parameter to 1550 K, yielding a distortion energy of 40 meV.

Ba2YMoO6 crystallises in the double perovskite structure (Fm3̄m; SG 225). In
the double-perovskite structure the Mo5+ (4$d^{1}$) ions are in the oxygen
octahedral surroundings.

The usual single-ion Hamiltonian for the octahedral crystal field with a
tetragonal distortion, completed by the s-o interaction term, takes a form in
the Stevens notation 14

$\hat{H}=B_{2}^{0}\hat{O}_{2}^{0}(L,L_{z})+B_{4}^{0}\hat{O}_{4}^{0}(L,L_{z})+B_{4}^{4}\hat{O}_{4}^{4}(L,L_{z})+\lambda_{s-o}\vec{L}\cdot\vec{S}$
(1)

where $\hat{O}_{2}^{0}$, $\hat{O}_{4}^{0}$ and $\hat{O}_{4}^{4}$ depend only
on $L$ and $L_{z}$ (for the ${}^{2}D$ subterm relevant to the 4$d^{1}$
configuration of the Mo5+ ion $L$=2 and $L_{z}$= 0, $\pm$1, $\pm$2) and where,
for instance, $\hat{O}_{2}^{0}$ = 3$L_{z}^{2}$ \- $L$($L$+1), and $S$=1/2. For
the exactly octahedral local symmetry $B_{4}^{4}$= 5$\cdot B_{4}^{0}$.

The effect of the octahedral CEF is well known - it splits 5 orbital states of
the ${}^{2}D$ subterm ($L$=2, $S$=1/2) into the lower orbital triplet $t_{2g}$
and higher orbital doublet $e_{g}$ separated by 10$Dq$ (= 24$\cdot
B_{4}^{4}$).

Figure 2: Calculated temperature dependence of the paramagnetic susceptibility of Ba2YMoO6 associated with the atomic-scale electronic structure of the Mo5+ (4$d^{1}$) ions (parameters are given in the figure, the realised eigenfunctions are collected in Table I) together with experimental data taken from Ref.2  (much similar data are in Refs 5 ; 8 ). The dashed-point line (red) shows $\chi$(T) for exactly octahedral CEF (10$Dq$ = 3.1 eV with the spin-orbit coupling $\lambda$ = +95 meV) whereas the dashed-two-point line (blue) shows $\chi$(T) without the spin-orbit coupling. The solid and dashed lines show $\chi$(T) with small off-octahedral distortion with $\lambda$= +95 meV for parallel (z) and perpendicular directions.  E [K] |  | +2 | -2 | +1 | -1 | 0 | Sz | Lz | Jz | mz
---|---|---|---|---|---|---|---|---|---|---
+21650 | 1 | +0.728$\uparrow$ | +0.685$\uparrow$ |  | +0.021$\downarrow$ |  | +0.4995 | +0.1196 | +0.6192 | +1.1199
36650 | 10 | +0.685$\downarrow$ | +0.728$\downarrow$ | +0.021$\uparrow$ |  |  | -0.4995 | -0.1196 | -0.6192 | -1.1199
+21650 | 5 |  |  | +0.037$\downarrow$ |  | 0.999$\uparrow$ | +0.4986 | +0.0014 | +0.5000 | +0.9998
36650 | 6 |  |  |  | +0.037$\uparrow$ | 0.999$\downarrow$ | -0.4986 | -0.0014 | -0.5000 | -0.9998
-13300 | 3 | +0.408$\downarrow$ | -0.408$\downarrow$ | +0.816$\uparrow$ |  |  | +0.1671 | +0.6667 | +0.8333 | +1.0000
1700 | 8 | -0.408$\uparrow$ | +0.408$\uparrow$ |  | +0.816$\downarrow$ |  | -0.1671 | -0.6667 | -0.8333 | -1.0000
-15000 | 4 |  |  | 0.999$\downarrow$ |  | -0.037$\uparrow$ | -0.4986 | +0.9986 | +0.5000 | +0.0002
0 | 7 |  |  |  | 0.999$\uparrow$ | -0.037$\downarrow$ | +0.4986 | -0.9986 | -0.5000 | -0.0002
-15000 | 2 | +0.603$\downarrow$ | -0.551$\downarrow$ | -0.577$\uparrow$ |  |  | -0.1671 | +0.4529 | +0.2858 | +0.1183
0 | 9 | -0.551$\uparrow$ | +0.603$\uparrow$ |  | -0.557$\downarrow$ |  | +0.1671 | -0.4529 | -0.2858 | -0.1183
Table 2: Realised spin-orbital states, energies and their eigenfunctions
($|LSL_{z}S_{z}>$ space, with $L$=2 and $S$=1/2) of the Mo5+ ion in Ba2YMoO6
calculated with the octahedral crystal-field parameter $B_{4}$= +300 K
(10$Dq$= 3.1 eV) and the spin-orbit coupling $\lambda_{s-o}$ = +1100 K.
$L_{z}$ values as +2, -2, +1, -1, and 0 contributions are shown, whereas
$S_{z}$ = +1/2 or -1/2 are shown by up and down arrows, respectively. In the
first column the absolute and the relative energies are given. The second
column is not informative - it shows a number of the row of the output results
with eigenvalues and eigenfunctions from our home-made computer program. In
last five columns the resulting magnetic characteristics from those
eigenfunctions are shown.

We think that we have found a nice reproduction of $\chi$(T) of Ba2YMoO6, see
Fig. 1. There $\chi^{-1}$(T) is presented, because an off-Curie-Weiss
dependence is more visible than in the $\chi$(T) plot. By introducing a small
tetragonal distortion, by changing only $B_{4}^{4}$ from 1500 K to 1550 K, we
have got an abrupt lowering of the $\chi^{-1}$(T) (or an abrupt increase of
$\chi$(T)) below T= 50 K exactly as is experimentally observed. Simultaneously
our calculations reproduce large lowering of the effective magnetic moment and
reproduce the observed Curie-Weiss temperature $\theta_{CW}$ of -(160-200 K).
In our calculations large value of $\theta_{CW}$ is an effect of crystal-field
and s-o interactions, not at all related to antiferromagnetic interactions.
The reproduction of the absolute value of $\chi$(T) proves that all Mo atoms
contribute equally to the macroscopic susceptibility and that all Mo atoms
occur in BYMo as the Mo5+ ions in the whole temperature range.

Fig. 2 illustrates influence of the s-o and off-octahedral CEf interactions on
temperature dependence of $\chi$(T) (of course, the change of $\chi$(T) is due
to the changing of the energies of 5 orbitals and their eigenfunctions, see
Table I). In the absence of s-o interactions (then 3 $t_{2g}$ orbitals are
degenerated) $\chi$(T) follows exactly the Curie law with the effective moment
of 2.26 $\mu_{B}$ (dashed-two-point line in Fig. 2). The s-o coupling in the
presence of the octahedral CEF removes 3-fold orbital degeneracy for a lower
doublet and higher by 1700 K orbital singlet, as is shown in Table II, where
also their eigenfunctions are collected. In this case the resulting $\chi$(T)
down to T= 50 K looks like a Pauli susceptibility (the dashed-point line on
the top of Fig. 2) having, however, so huge effective moment as 6.0 $\mu_{B}$.
The inferred distortion $B_{4}^{4}$, which could be produced in the reality,
for instance, by a slight shortening in-plane oxygen distances, causes
anisotropy of the local susceptibility (the experimentally realised $\chi$(T)
direction is along the larger susceptibility shown in Fig. 2 as the solid
line, which well follows the experimental points read from Ref. 2 ) and
removes double orbital degeneracy yielding three Kramers-doublet states at
energies 0, 38 and 170 meV, see Table I.

The Kramers-doublet ground state is characterised by largely reduced both spin
(sz= $\pm$0.295) and magnetic moment (mz= $\pm$0.258 $\mu_{B}$). One sees,
that the spin is reduced by the factor of 2, whereas the magnetic moment much
more, by the factor of 4. Such large reduction of the magnetic moment is
caused by a substantial orbital moment of mz= $\pm$0.33 $\mu_{B}$). This
ground-state magnetic moment largely determines the lowest-temperature
susceptibility yielding in the temperature range up to 25 K the Curie-Weiss
law with $p_{e}ff$= 0.574 $\mu_{B}$. This value is in perfect agreement with
the experimental value of $p_{e}ff$= 0.57 $\mu_{B}$ 8 .

Table I collects the spin-orbital states, energies and their eigenfunctions,
calculated within the $|LSL_{z}S_{z}>$ space, with $L$=2 and $S$=1/2, relevant
to the Mo5+ ion in Ba2YMoO6 under the action of the octahedral crystal field,
$B_{4}^{0}$= +300 K and $B_{4}^{4}$= +1500 K, of the relativistic on-site
spin-orbit coupling $\lambda_{s-o}$ = +1100 K and a small tetragonal
distortion obtained to be realised by extra $B_{4}^{4}$ = +50 K (so, finally
$B_{4}^{4}$ = +1550 K). An occurrence of a local distortion of MoO6 octahedra,
despite of the overall cubic symmetry, has been also infreed fro infrared
transmission spectroscopy experiments 8 . Table II collects the spin-orbital
states and their characteristics for the purely octahedral crystal field.
Table II is presented for noticing the effect of the tetragonal distortion, if
compared with Table I. One can notice a substantial increase of the local
magnetic moment and the spin of the CEF+s-o ground state by the tetragonal
distortion, to mz= $\pm$0.258 $\mu_{B}$ and sz= $\pm$0.295 from mz= $\pm$0.118
$\mu_{B}$ and sz= $\pm$0.167. Thus our calculations prove that distortions
usually help the formation of the magnetic order.

Inspecting Table I one notices that each state has large orbital moment, which
only partly cancels the spin moment (which on the other side can take itself
quite small values). The ground state moment amounts to 0.25 $\mu_{B}$, but in
both directions. Interesting is the diminishing anisotropy of the
susceptibility when temperature approaches zero temperature - it is
responsible for so persistent nonmagnetic state of BYMo or, in other words,
for formation of an amorphous spin state of atomic-scale Mo ions, because the
moments cannot establish one preferred direction. However, due to that the
Kramers degeneracy of the ground state has not been yet removed we expect
enormous low-temperature specific heat resembling heavy-fermion phenomena in
anomalous Ce, U and Yb intermetallics. In both of these type compounds the
intra-atomic relativistic s-o coupling plays the fundamental role for the
final determination of the shape of the ground-state eigenfunctions with
sophisticated values of the spin and orbital moments. The relativistic s-o
interactions has allowed for the revealing of the orbital magnetic moment in
NiO so large as 0.54 $\mu_{B}$ (20% of the total moment) 21  and as 0.80
$\mu_{B}$ associated with the Fe2+ ions in FeBr2 18 .

Although we have got a nice reproduction of $\chi$(T) in the whole temperature
range with revealing the origin of the dramatic reduction of the local Mo
moment with approaching the lowest temperatures we would like to mention the
experimental specific heat of BYMo, Fig. 4 of Ref. 2 , which has revealed a
Schottky-like peak with the heat maximum at about 80 K. Assuming that this
peak is really the Schottky-like one one derives the existence of the local
state at about 200 K. In our case we have got the first excited orbital
Kramers doublet at 440 K. At present, we did not succeed to get local
parameters which would recoincile the $\chi$(T) and c(T) dependences but we
think that within our approach it will be possible by introducing lower
symmetry CEF parameters. If one is at this moment not convinced about our
present contribution to the theoretical description of $\chi$(T) experiments
on BYMo he is welcome to remember that all this discussion is going within 20
meV energy, whereas up to now the electronic structure has been discussed in
the best case in a 0.3 eV-energy scale.

##  IV Conclusions

We succeeded in remarkably nice theoretical description of the temperature
dependence of the paramagnetic susceptibility of Ba2YMoO6 in the whole
temperature range as originating from the atomic-scale discrete electronic
structure of Mo5+ ions in the 4$d^{1}$ configuration resulting from crystal-
field and spin-orbit interactions. A strong violation of the Curie-Weiss law
and a nonmagnetic state experimentally observed down to 50 mK, despite of s=
1/2 spins, results from the on-site CEF + spin-orbit interactions. The
relativistic s-o and crystal-field interactions cause the dramatic decrease of
the spin and magnetic moment of the Mo5+ (4$d^{1}$) ions. Our results indicate
that BYMo should be described as a spin-glass system with incoherent intersite
local distortions in which exchange interactions become too weak to produce a
magnetic order above 50 mK due to the dramatic lowering of the Mo-ion magnetic
moment. We treat our single-valence description, with all Mo ions being in the
Mo5+ state, as much physically superior over a description with many mixed Mo-
ion valences 13 .

email<EMAIL_ADDRESS>www.actaphysica.eu

## References

  * (1) E. J. Cussen, D. R. Lynham, and J. Rogers, Chem. Mater. 18, 2855 (2006).
  * (2) M. A. de Vries, A. C. Mclaughlin, and J. W. G. Bos, Phys. Rev. Lett. 104, 177202 (2010).
  * (3) K. Stitzer, M. Smith, and H. -C. zur Loye, Solid State Sci, 4, 311 (2002).
  * (4) A. S. Erickson, S. Misra, G. J. Miller, R. R. Gupta, Z. Schlesinger, W. A. Harrison, J. M. Kim, and I. R. Fisher, Phys. Rev. Lett. 99, 016404 (2007).
  * (5) T. Aharen, J. E. Greedan, C. A. Bridges, A. A. Aczel, J. Rodriguez, G. MacDougall, G. M. Luke, T. Imai, V. K. Michaelis, S. Kroeker, H. Zhou, C. R. Wiebe, and L. M. D. Cranswick, Phys. Rev. B 81, 224409 (2010).
  * (6) J. P. Carlo, J. P. Clancy, T. Aharen, Z. Yamani, J. P. C. Ruff, J. J. Wagman, G. J. Van Gastel, H. M. L. Noad, G. E. Granroth, J. E. Greedan, H. A. Dabkowska, and B. D. Gaulin, Phys. Rev. B 84, 100404(R) (2011).
  * (7) M. de Vries, J. Piatek, M. Misek, J. Lord, H. Rnnow, and J. -W. G. Bos, New Journal of Physics 15, 043924 (2013).
  * (8) Z. Qu, Y. Zou, S. Zhang, L. Ling, L. Zhang, and Y. Zhang, J. Appl. Phys. 113, 17E137 (2013).
  * (9) W. M. H. Natori, E. C. Andrade, E. Miranda, and R. G. Pereira, Phys. Rev. Lett. 117, 017204 (2016).
  * (10) L. Xu, N. A. Bogdanov, A. Princep, P. Fulde, J. van den Brink, and L. Hozoi, npj Quant Mater. 1, 16029 (2016).
  * (11) J. Romhanyi, L. Balents, and G. Jackeli, Phys. Rev. Lett. 118, 217202 (2017).
  * (12) P. M. M. Thygesen, J. A. M. Paddison, R. Zhang, K. A. Beyer, K. W. Chapman, H. Y. Playford, M. G. Tucker, D. A. Keen, M. A. Hayward, and A. L. Goodwin, Phys. Rev. Lett. 118, 067201 (2017).
  * (13) J. Karp, M. Bramberger, M. Grundner, U. Schollwock, A. J. Millis, and M. Zingl, Phys. Rev. Lett. 125, 166401 (2020).
  * (14) R. J. Radwanski, R. Michalski, and Z. Ropka, Acta Physica Polonica B 31, 3079 (2000).
  * (15) A. Abragam and B. Bleaney, Electron Paramagnetic Resonance of Transition Ions (Clarendon Press, Oxford, 1970); ch. 7.
  * (16) G. Jackeli and G. Khaliullin, Phys. Rev. Lett. 102, 017205 (2009); G. Jackeli and G. Khaliullin, Phys. Rev. Lett. 103, 067205 (2009).
  * (17) R. J. Radwanski, N. H. Kim-Ngan, F. E. Kayzel, J. J. M. Franse, D. Gignoux, D. Schmitt, and F. Y. Zhang, J. Phys. Condens. Matter 4, 8853 (1992).
  * (18) Z. Ropka, R. Michalski, and R. J. Radwanski, Phys. Rev. B 63, 172404 (2001).
  * (19) Z. Ropka and R. J. Radwanski, Phys. Rev. B 67, 172401 (2003).
  * (20) R. J. Radwanski, R. Michalski, and Z. Ropka, Physica B 312-313, 628 (2002).
  * (21) R. J. Radwanski and Z. Ropka, Acta Physica Polonica A 97, 963 (2000).


# The implementation and optimization of Bitonic sort algorithm based on CUDA

Qi Mu. Liqing Cui. Yufei Song

###### ABSTRACT

This paper describes in detail the bitonic sort algorithm,and implements the
bitonic sort algorithm based on cuda architecture.At the same time,we conduct
two effective optimization of implementation details according to the
characteristics of the GPU,which greatly improve the efficiency. Finally,we
survey the optimized Bitonic sort algorithm on the GPU with the speedup of
quick sort algorithm on the CPU.Since Quick Sort is not suitable to be
implemented in parallel,but it is more efficient than other sorting algorithms
on CPU to some extend.Hence,to see the speedup and performance,we compare
bitonic sort on GPU with quick Sort on CPU. For a series of 32-bit random
integer,the experimental results show that the acceleration of our work is
nearly 20 times.When array size is about $2^{16}$,the speedup ratio is even up
to $30$.

KEY WORDS: Parallel; CUDA; Quick Sort; Bitonic sort

##  1 INTRODUCTION

Sorting is a well-studied topic and a fundamental problem in computer science.
Sorting appears as an internal step of many programs and processes. According
to statistics data, it indicates that more than 25 percent of time CPU costs
is sorting.Hence, efficient sorting algorithm and implementations are one of
the basic keys for performance improvement. Currently, there are many sorting
algorithms,Bubble sort,Odd-even sort, Insertion sort,Heap sort,Selection
sort,Radix sorting,Quick sort,Bitonic sort…
There is an increasing need for programming to address parallelism on a
variety of architectual approaches. However,due to some physical limitation of
technology and material,the frequency level of single-core CPU is also
influenced.Today’s graphics cards contain very powerful multi-core processors.
The processors are specialized for compute-intensive,highly parallel
computations. They could be used to assist the CPU in solving problems that
can be efficiently data-paralleled. Meanwhile,more and more and more
scientists,researchers and software developers are using GPU to accelerate
their algorithms and application.
GPU [1] processing power and memory bandwidth has obvious advantages in terms
of cost and power consumption do not need to pay too much with respect to
CPU.Due to the highly parallel graphics rendering, making GPU processing power
can be increased by increasing the memory bandwidth and parallel processing
unit and a memory unit of the control mode. GPU designers put more transistors
used as an execution unit, not like CPU used as a complex control unit and a
cache in order to improve the efficiency of a small number of execution units.
CPU integer calculation, branch, logic and floating-point arithmetic
operations are performed by different units, in addition to a floating point
accelerator. Thus, CPU havs different performance to compute tasks of
different types. While GPU computes tasks of different types by an integer and
floating-point arithmetic unit. so the power of integer GPU computes is
similar to its floating-point capability. At present, the mainstream has
adopted a unified architecture GPU unit. In adition, GPU computing has a huge
advantage, its memory subsystem. The Ultra-high bandwidth of memory not only
makes tremendous floating-point capability keep a stable throughput, but
guarantees the efficient operation of data-intensive tasks. This is the why
GPU is more and more popular in video games, Film industry, industrial design,
medical imaging, space exploration, telecommunications, etc.

##  2 CUDA

###  2.1 CUDA Program Model

CUDA [2] is a parallel programming model and software environment for NVIDIA’s
GPUs. CUDA allows programmer to program kernels executed on the GPU. When one
kernel is executed in parallel, a number of CUDA threads are also created at
the same time. Threads are organized in warps and blocks. A warp is a fixed-
size group of threads (32 threads on GPUs with compute capability 1.x, 2.0),
while a block consists of up to 16 (1.x) or 32 (2.0) warps. The number of
warps per block and the number of blocks (and therefore the total number of
threads) can be specified for each individual kernel launch. When calling
kernal, the data transmitted from the CPU main memory to the GPU memory, and
is sent back to CPU after calculation process.

###  2.2 CUDA Memory Model

A thread that executes a kernel has access to high-performance thread local
registers with the same lifetime as the lifetime of the thread. Each thread
block has a shared memory (visible to all threads of the block) with the
lifetime of the block. The shared memory is as fast as registers, as long as
there are no conflicting accesses by threads of the same block. All threads
have random access to global memory that has the lifetime of the application,
i.e. data reside in global memory for several kernel launches. In general,
global memory is significantly slower than access to registers or shared
memory. When threads of the same half-warp (upper or lower 16 threads of a
warp) access the global memory simultaneously, the accesses can be coalesced
into a single memory transaction, if the accessed memory lies in the same
global memory segment. If not, the simultaneous access causes multiple
sequential memory transfers.
CUDA provides a barrier mechanism for synchronization of threads of the same
block, but does not provide any synchronization mechanisms for threads of
different blocks. In general,the host is not blocked when launching a kernel,
but it can explicitly wait for a kernel to finish(host synchronization). Since
variables reside in global memory for multiple kernel launches, tasks can be
synchronized by distributing the workload to multiple kernel launches using
host synchronization, while keeping information in global memory.
Efficient synchronization and memory access are two of the most important
factors that influence the performance of a GPGPU application. Since kernel
launches cause small startup delay, the application should use CUDA ’s barrier
mechanism instead of host synchronization or even better should avoid explicit
synchronization. The host can only access the global memory of the GPU, thus
global memory access cannot be completely avoided, but the global memory
accesses of a kernel should be reduced to a minimum.

Figure 1:  cuda memory model

##  3 BITONIC SORT ALGORITHM

###  3.1 Algorithm Description

Bitonic sort [3] is a binary merge sort, used in the parallel processor with a
good parallel performance.Bitonic sequence, for example, 1,5,9,10,12,8,7,2 is
a bitonic sequence,the first half of sequence is Monotonically increasing, the
second half of this sequence is Monotonically decreasing.
Similarly,12,8,7,2,1,5,9,10is also a bitonic sequence. Let A be an arbitrary
input sequence to sort and let $n=2^{k}$ be the length of A. The process of
sorting A then consists of k phases. The subsequences of length 2

(A[0],A[1]),(a[2],A[3]),…,A[n-2],A[n-1])

are bitonic sequences by definition. In the first phase these subsequences are
sorted using bitonic merge (as shown above) alternating descending and
ascending, which makes the subsequences of length 4 bitonic sequences:

(A[0],A[1],A[2],A[3]),…,(A[n-4],A[n-3],A[n-2],A[n-1])

In the second phase these subsequences of length 4 are sorted alternating
descending and ascending, resulting in subsequences of length 8 being bitonic
sequences. In the r th phase of bitonic sort the total number of subsequences
being sorted is $2^{k-r}$ and the length of each of these subsequences is
$2^{r}$. Sorting a sequence of length $2^{r}$ using bitonic merge consists of
r steps. After the $[k.1]$th phase sequence A is a bitonic sequence. A is
sorted in the last phase k. Figure 2 shows the bitonic sorting network for
input sequences of size $n=8$. The sorting network consists of $3=log8$
phases, phase p having p steps. Every step consists of $4=n/2$
compare/exchange operations.

Figure 2:  Bitonic sorting network

###  3.2 Algorithm Analysis

To our knowledge,quick sort algorithm often used for sorting on the CPU, and
high efficiency.Although the average time complexity of quick sort is
$O(nlogn)$,it may not be suitable for implementation on the GPU.
Based on the introduction of section 3.1,we know the time complexity of
bitonic sort is $O(n(logn)^{2})$, it is very suitable for implementation on
the GPU for the reasons: the fixed and stable procedure, data independent, and
no data cross in each round. It needs
$\sum_{i=1}^{logn}i=\frac{logn(logn+1)}{2}$ rounds, at the same time, it
completes $\sum_{i=1}^{logn}i*\frac{n}{2}=\frac{nlogn(logn+1)}{4}$ compare-
exchange operations.

###  3.3 Implementation

Suppose an array of length $N=2^{k}$, the sequence needs $k(k+1)/2$ rounds.
And each round calls a kernal function, which completes $N/2$ compare-exchange
operations.All these operations can be done by a lot of threads created
parallelly, after the end, data is synchronized. Successively,then next kernal
is called.

##  4 OPTIMIZATIONS

The method described above has two inevitable shortcomings: too many kernals
are called, access global memory consumes too much time, i.e. a long delay.So
we consider these two aspects of optimization, reducing the number of memory
accesses and the number of kernal launches[4].

###  4.1 Optimization1:using shared memory

In our implementation,every step needs $N/2$ threads,and each thread completes
a compare-exchange operation.If a subsequence od length $2^{s}$ processed in
step $s$ completely fits into shared memory of one block,the block first
transfers the subsequence into shared memory.That way we can process the steps
s,s-1,…,1 using the shared memory and efficient block-synchronization
mechanisms between multisteps instead of host-synchronization in one kernal
function.So time of accessing global memory is decreased.

###  4.2 Optimization2:using the register

Another optimization is using register.According to the Figure 2,in the phase
3, $4$ address of [0,2,4,6] are used in two successive step 3 and step 2 that
we can attach to the register, so that we can reduce two steps to one.
Analogously,using this methods can reduce the number of kernal launches and
the number of access global memory.

##  5 RESULTS

The test result is obtained on a experimental Platform with CPU(Intel Xeon
E5-2620) and GPU(Kepler architecture, K10). Experimental data is 32-bit random
integer from 128K to 256M. As shown in below table1.

[b]  CPU Times(ms) GPU BitonicSort Times(ms) Array size QuickSort BitonicSort
Basic Semi Optimized Ratio $128K$ $-$ $30.00$ $0.76$ $0.46$ $0.36$ – $256K$
$20.00$ $60.00$ $1.21$ $0.87$ $0.66$ 30.2 $521K$ $30.00$ $110.00$ $2.22$
$1.78$ $1.31$ 22.7 $1M$ $80.00$ $250.00$ $4.58$ $3.89$ $2.80$ 28.5 $2M$
$150.00$ $550.00$ $8.90$ $7.95$ $5.87$ 25.5 $4M$ $280.00$ $1230.00$ $18.14$
$16.59$ $12.30$ 22.7 $8M$ $590.00$ $2670.00$ $38.13$ $35.29$ $26.36$ 22.3
$16M$ $1230.00$ $5880.00$ $80.09$ $75.52$ $56.27$ 21.8 $32M$ $2570.00$
$12900.00$ $173.77$ $162.56$ $120.93$ 21.3 $64M$ $5360.00$ $27780.00$ $373.52$
$350.87$ $258.61$ 20.7 $128M$ $11180.00$ $59860.00$ $803.16$ $756.94$ $553.49$
20.1 $256M$ $23260.00$ $128660.00$ $1727.23$ $1631.92$ $1185.02$ 19.6

Table 1: results

  * 1

Basic : no optimized.

  * 2

Semi : optimization1.

  * 3

Optimized : optimization1 and optimization2.

  * 4

Ratio : acceleration ratio = Times(CPU Quick Sort)/Times(GPU Bitonic Sort)

Table 1 shows that quick sort is much faster than bitonic sort on the CPU,this
is due to the differences in time complexity. After our two step
optimization,the efficiency of bitonic sort algorithm on the GPU has increased
significantly.

##  6 CONCLUSION AND FUTURE WORK

We have presented an efficient implementation of bitonic sort based on CUDA.To
achieve this we carefully optimized our implementation with respect to the
number of accesses to global memory and the number of kernal launches. For a
series of 32-bit random integer,our experimental results show that the
acceleration of CPU/GPU is more than 20 times.
Due to our current work is a relatively simple, so we will strengthen two main
directions in the future work.The first is to test different types of
data,such as 64-bit integer,32-bit float,64-bit double.The second is to
further explore and compare the performance of a multicore GPU bitonic sort
implementation.

## References

  * [1] E.Lindholm,J. Nickolls,S. Oberman,and J. Montrym: NVIDIA Tesla:A unified graphs and computing architecture.IEEE Micro,28:39-55,Mar/Apr 2008.
  * [2] J. Nickolls,I. Buck,M. Garland,and K. Skadron: Scalable parallel programming with CUDA,6(2):40-53,Mar/Apr 2008.
  * [3] http://en.wikipedia.org/wiki/Bitonic_sorter.
  * [4] Hagen Peters,Ole Schulz-Hildebrandt and Norbert Luttenberger: Fast in-place,comparison-based sorting with CUDA:a study with bitonic sort.PPAM09:Processings of the International Conference on Parallel Processing and Applied Mathematics,Wroclaw,Poland,September 2009.