Original Article Journal Section [2]Equally contributing authors. Juan Pablo
Mesa López, School of Applied Sciences and Engineering, Universidad EAFIT,
Medellín, Antioquia, 05001, Colombia<EMAIL_ADDRESS>
# A Bi-Objective Approach to Last-Mile Delivery Routing Considering Driver
Preferences

Juan Pablo Mesa  Alejandro Montoya  Raul Ramos-Pollán  Mauricio Toro

###### Abstract

The Multi-Objective Vehicle Routing Problem (MOVRP) is a complex optimization
problem in the transportation and logistics industry. This paper proposes a
novel approach to the MOVRP that aims to create routes that consider drivers’
and operators’ decisions and preferences. We evaluate two approaches to
address this objective: visually attractive route planning and data mining of
historical driver behavior to plan similar routes. Using a real-world dataset
provided by Amazon, we demonstrate that data mining of historical patterns is
more effective than visual attractiveness metrics found in the literature.
Furthermore, we propose a bi-objective problem to balance the similarity of
routes to historical routes and minimize routing costs. We propose a two-stage
GRASP algorithm with heuristic box splitting to solve this problem. The
proposed algorithm aims to approximate the Pareto front and to present routes
that cover a wide range of the objective function space. The results
demonstrate that our approach can generate a small number of non-dominated
solutions per instance, which can help decision-makers to identify trade-offs
between routing costs and drivers’ preferences. Our approach has the potential
to enhance the last-mile delivery operations of logistics companies by
balancing these conflicting objectives.

###### keywords:

Last-mile delivery operations, Historical patterns, Visual attractiveness
metrics, Data mining, Bi-objective problem, Pareto front

##  1 Introduction

The Vehicle Routing Problem (VRP) is a well-known problem in the area of
combinatorial optimization with applications in many different fields, such as
logistics and transportation. It has been studied extensively and has many
variations, and extensions [1]. The VRP is usually formulated as a single-
objective problem, where the goal is to route a fleet of vehicles to satisfy
the demand of a set of customers while minimizing operations costs. However,
in many real-world applications, multiple objectives need to be taken into
account. This leads to a variation of the VRP called the multi-objective VRP
(MOVRP).

Traditionally, the objective has been to minimize costs. However, in recent
years, a new objective in the industry has emerged that focuses on creating
routes that take into account the decisions and preferences of operators and
drivers. This objective has been explored from two approaches. Firstly, there
has been research on route planning that aims to design visually attractive
routes. This line of research is evident in the multiple works [2, 3, 4, 5].
Secondly, another approach involves utilizing data or pattern mining
techniques to analyze the historical behavior of drivers and leverage this
information for planning similar routes. This approach is evident in the works
[6, 7, 8, 9, 10, 11, 12].

The first approach consists of planning vehicle routes that consider visual
attractiveness factors. Even if visual attractiveness is subjective and
difficult to express in mathematical terms, in a review by Rossit et al. [4],
the authors identify three features that can be computed, which attractive
vehicle routes should have: Compactness or proximity of the customers visited.
Not overlapping or crossing of routes with each other. Low complexity of each
single route. In this work, we further explore this strategy within a last-
mile delivery context.

In 2021, Amazon and MIT’s Center for Transportation & Logistics initiated the
Amazon Last-Mile Routing Research Challenge [13]. The goal was to incorporate
the know-how of Amazon’s experienced drivers into route planning. Participants
aimed to predict delivery routes just as these drivers would execute them.
According to the results of different participant teams [14, 15, 16], the
objective can be attained by using information extracted from routes
previously executed by drivers and integrating this information within
traditional routing heuristics. We consider that this approach works well when
the executed routes come from planned routes where a route optimization
software was involved at some point during the planning. On the contrary, in
cases where the real operation routes of a company are planned without any
routing optimization step, the route sequences predicted using the previous
methods will probably vary considerably with respect to the executed routes.

These objective functions that consider the preferences of drivers and
decision-makers can be highly beneficial to them. However, this may conflict
with the interests of the company, whose primary concern is minimizing
transportation costs. As a result, a bi-objective problem naturally arises,
which involves minimizing costs while simultaneously taking driver preferences
into account. There is already a study in the literature that addresses this
issue, namely the work of Rocha et al. (2022) [5]. In their research, the
authors propose the utilization of clustering as an indicator of visual
attractiveness in vehicle routing problems. They present a bi-objective
vehicle routing problem, examine the impact of visual attractiveness on
routing costs, and evaluate the effectiveness of their proposed algorithms in
generating VRP solutions that are both cost-effective and visually appealing.

Following the idea of Rocha et al. we propose a bi-objective vehicle routing
problem that instead of focusing on the visual attractiveness of the routes,
it leverages information from historical routes to plan routes that consider
drivers’ preferences, while also aiming to minimize routing costs. This work
is divided into two parts to demonstrate the potential of our method.

We evaluate and compare both approaches in the first part of this work. We use
the model proposed by Mesa et al. [12] for the approach of leveraging
historical data. Additionally, we modify some of the components of this same
model to compute visual attractiveness features and use them in the objective
function of the other approach. We evaluate both strategies using the dataset
and evaluation metrics provided by Amazon [17, 13] to determine which method
works better for predicting the order in which drivers will execute route
sequences. We demonstrate that the method that leverages data from historical
routes is more effective than using visual attractiveness metrics.

In the second part of this work, we argue that minimizing routing cost and
favoring the routing order of historical routes executions can be conflicting
in last-mile delivery. Drivers might be more inclined to execute routes in a
familiar or known way. In contrast, the company might prefer that the routes
are executed with the minimum possible cost. This leads to a bi-objective
problem, with planned routes that are as similar as possible to previous
routes and routes that aim to minimize routing costs. Therefore, we formulate
the bi-objective problem, and then we propose a simple two-stage GRASP
algorithm with heuristic box splitting to solve it. The proposed algorithm
aims to approximate the Pareto front [18] and present routes that cover a wide
range of the objective function space to help decision-makers.

The contributions of this work are summarized as follows:

  * •

We explore using visual attractiveness metrics to predict route sequences in
the way drivers would execute them.

  * •

We compare the use of visual attractiveness metrics versus data mining
strategies for predicting last-mile delivery routes.

  * •

We propose a relevant last-mile delivery bi-objective problem for decision-
makers.

  * •

We propose a two-stage GRASP algorithm with heuristic box splitting to
approximate the Pareto front of this problem.

The remainder of this paper is organized as follows. Section 2 discusses a
brief literature review on multi-objective last-mile delivery problems and the
prediction of delivery route sequences. Section 3 explains the motivation for
this research. Section 4 provides a description of the routing problems that
this work studies. Section 5 describes the three approaches studied in this
work and their respective algorithms. Section 6 presents the computational
experiments and discusses the results. Finally, Section 7 presents the
conclusions and future work directions.

##  2 Related works

This section of the research article will provide an overview of existing
methods and studies related to considering decision-makers and drivers’
preferences and incorporating cost considerations along with these preferences
when planning routes. The review is divided into two parts, each addressing a
specific aspect of the problem. The first part explores methods that focus on
incorporating drivers’ preferences into route planning, while the second part
reviews the literature on methods that consider both travel costs and
decision-makers’ preferences. By reviewing relevant research studies, we aim
to identify gaps, challenges, and opportunities for further research in the
context of multi-objective last-mile delivery problems and the prediction of
delivery route sequences.

In recent years, incorporating the preferences of drivers or decision-makers
has gained attention within the vehicle routing community. Several authors
have explored an approach that relies on the opinion that practitioners may
drive their logistics decisions based on aesthetical considerations [2, 4, 5].
Tang and Miller-Hooks [2] proposed planning visually attractive routes that
are compact and do not cross one another. This was achieved by considering
some shape measures that would improve the visual layout of routes while
maintaining satisfactory results in terms of “conventional" VRP measures. A
comprehensive review of different visual attractiveness metrics used in VRP
can be found in Rossit et al. [4].

Another approach has been recently focused on the last-mile delivery industry,
where routes are planned considering data related to the drivers and decision-
makers’ experience, knowledge, and preferences [7, 19, 11]. Multiple works
have proposed strategies to use data from past routes to learn the routing
preferences of drivers [6, 8, 9, 20, 10]. These works rely on GPS data and
Markov models to predict future route decisions. Krumm [6] predicts the next
route segment that a driver will take in the current route, while the works of
Wang et al. [8] and Neto et al. [20] go a step further and aim to predict the
remaining part of a route and the destination.

Similarly, Canoy and Guns [10] also used weighted Markov models to plan
multiple routes in VRP considering driver’s preferences, while avoiding the
need to specify drivers’ preferences explicitly. Additionally, Canoy et al.
[16] extended the use of Markov models and transition probabilities together
with distance information to solve the Traveling Salesman Problem (TSP) in the
Amazon Routing Challenge [13]. The approach of Canoy et al. focuses on vehicle
route planning with constraints and how to learn the implicit constraints that
can be extracted from past routes. Other solutions that were presented in the
Amazon Challenge also leveraged the use of constraints extracted from past
routes. Some notable ones were the solutions of the three winning teams of
Cook et al. [15], Guo et al. [21], and Arslan and Abay [22]. All these
solutions encoded patterns found in the data either as clustering constraints,
precedence, or antecedence constraints and coupled them with powerful local
search heuristics. Among the other competitors, Mesa et al. [14] also used the
strategy of constructing a probability transition matrix based on historical
trips taken by a driver to predict a route. However, instead of using the
extracted information as constraints for route planning, the information was
considered within the objective function as penalizations. Instead of relying
on a Markov model to calculate probabilities, the authors used vector fields
to represent historical data. They used the directions and magnitudes of the
vectors to calculate the transition probabilities. Additionally, they
developed a two-stage approach for the TSP where one first calculates the
transition probabilities between clusters based on historical data and then
uses these probabilities, along with other factors, within a TSP metaheuristic
to predict route sequences. This approach was later improved by the same
authors in [12].

The problem becomes more complex when multiple objectives need to be
considered, such as minimizing transportation costs, reducing vehicle
emissions, improving customer satisfaction, or including the drivers’
preferences. These objectives may conflict with each other, and finding an
optimal solution that satisfies all of them is challenging [23]. Multi-
objective optimization techniques aim to find multiple solutions, called
Pareto optimal solutions, that represent trade-offs between conflicting
objectives. This approach provides decision-makers with a range of
alternatives from which they can choose the most suitable solution based on
their preferences.

Some companies aim to optimize vehicle routes to reduce costs, while also
considering consistency in visiting each customer at similar times of the day
with familiar drivers to increase customer satisfaction and lifetime value.
Lian et al. [24] proposed a multi-objective strategy to reduce the maximum
number of drivers visiting a customer to improve driver consistency and
minimize the maximum difference between the earliest and latest arrival time
to achieve arrival time consistency. The authors included consistency
considerations such as driver and time consistency as objectives in their
multi-objective variant of the Consistent Vehicle Routing Problem. By
including these considerations as objectives rather than constraints, their
approach can help decision-makers balance competing objectives effectively,
allowing them to choose between a minimal-cost routing plan with poor
consistency, a high-cost plan with perfect consistency, or any plan in
between.

Rocha et al. [5] explore delivery cases where incorporating visual
attractiveness into route planning can lead to more favorable routes for
drivers and present a bi-objective vehicle routing model to support their
proposal. They proposed the use of customer clustering as a proxy for visual
attractiveness in vehicle routing problems, assessing the impact of visual
attractiveness on routing cost and evaluating the effectiveness of NSGA-II
heuristics in producing low-cost and visually attractive VRP solutions.

In conclusion, there are two research gaps that this work aims to address.
Firstly, while previous works have separately explored methods based on data
from driver experience or visual attractiveness metrics for route planning,
this paper stands out as the only one that directly compares these two
approaches using routes from a real-life last-mile delivery operation. By
conducting this comparison, we can determine which method yields better
results in terms of route prediction based on driver preferences. Secondly,
this work introduces a novel perspective by formulating the problem as a bi-
objective optimization task that considers both the total cost of the route
and the objective of predicting the route based on historical data and machine
learning models. To the best of our knowledge, this is the first study that
tackles a bi-objective problem in this context, offering decision-makers
valuable insights into optimizing their route planning process.

##  3 Motivation

In 2021, Amazon and MIT’s Center for Transportation & Logistics launched the
2021 Amazon Last-Mile Routing Research Challenge [13]. The challenge focused
on improving last-mile delivery by addressing the issue of drivers deviating
from planned routes. Traditionally, route planning does not consider driver
preferences and their tacit knowledge of the areas they traverse. To bridge
this gap, the challenge encouraged the use of machine learning to devise or
predict route sequences that mirror the choices of seasoned drivers.

The routes in the dataset provided in the challenge are characterized for
being a TSP solution that ends at the same station from which it departed,
with asymmetric travel times, given that the TSP instances come from real-
world routing operations, and each stop has an assigned zone ID. In previous
studies from [22, 25, 12], different authors determined that the drivers
usually respect the zone ordering and visit the stops with the same zone ID
contiguously before visiting other stops with a different zone ID. This
routing with zones is considered in the literature as a Clustered TSP.
Therefore, in this work, we are going to address the routing problem as a
Clustered Asymmetric TSP.

The research challenge pushed participants to propose techniques that aimed to
capture the decision-making and tacit knowledge of seasoned delivery drivers
and generate route sequences that a driver would consider convenient,
efficient, and safe. However, we consider that even if the challenge affronted
a reality in last-mile delivery routing, it still lacked considering the
minimization of the total cost of the route, which is one of the main
objectives that the industry considers. For example, using data from the
Amazon Last-Mile Routing Research Challenge and a VRP solver, we obtain the
routes observed in Figure 1.

In Figure 1(a), we show the solution with the minimal traveling time cost that
the solver was able to find for this problem, with a traveling time of 172
minutes. A colored circle represents each stop, and each color represents a
different zone ID. However, an Amazon driver performed the route as it is
shown in Figure 1(b). The traveling time cost of this route is 177 minutes.
The drivers or the route planners might give different arguments, such as this
route is safer or easier to execute, which can justify the extra travel-time
cost. A third possible route is shown in Figure 1(c), with a traveling time
cost of 174.5 minutes, which is longer than the first route, but shorter than
the way in which an Amazon driver performed the route. One can argue that this
route is more similar to the Amazon route than the route in Figure 1(a) in
terms of a similarity score metric that will be explained later in this work.

(a) Planned Clustered ATSP route with minimized traveling time of 172 minutes.

(b) Clustered ATSP route as planned by Amazon with traveling time of 177
minutes.

(c) Planned Clustered ATSP route similar to Amazon’s route with traveling time
of 174.5 minutes.

Figure 1: Different routing solutions for the same instance of the Amazon
dataset.

It can be observed that predicting or planning the routes that consider the
preferences of the drivers does not necessarily translate into minimizing the
cost of the routes. Therefore, the problem that we aim to solve should
consider both objectives simultaneously. However, it is important to determine
how to measure the amount of “preference" that a driver has toward a route. In
the literature, we found two approaches. One approach is the use of visual
attractiveness metrics [4] to evaluate and create routes based on the
“attractiveness” of the routes to drivers and decision-makers. The other
approach consists of using data mining strategies to extract historical
patterns about the order in which stops were visited along the routes; the
patterns are then used as constraints or penalizations to guide the route
planning.

##  4 Problem description

The Clustered Asymmetric Traveling Salesman Problem Problem (CATSP) can be
formulated using graph notation as follows. Let’s consider a directed graph
$G=(V,A)$, where $V$ is the set of vertices representing stops and $A$ is the
set of arcs representing directed edges between stops. There is one depot
station $s$, from which a route must start and end, $n$ stops, and $m$
clusters.

Let $Z$ be the set of clusters, and $V_{g}$ be the set of vertices in cluster
$g$ for $g=1,2,\ldots,m$, such that $\bigcup_{k=1}^{m}V_{g}=V$ and $V_{i}\cap
V_{j}=\emptyset$ for $i\neq j$. Let $a_{ij}$ be the arc from vertex $i$ to
vertex $j$ for $i,j\in V$ and $i\neq j$. The arc set $A$ is defined as
$A=\\{(i,j):i,j\in V,i\neq j\\}$, and the cost of traversing an arc $(i,j)$ is
$c_{ij}$. Let $x_{ij}$ be a binary decision variable indicating whether arc
$(i,j)$ is included in the tour. $x_{ij}=1$ if arc $(i,j)$ is selected, and
$x_{ij}=0$ otherwise.

The “traditional" objective of the CATSP is to minimize the total cost of the
tour, which can be represented as:

$\min f_{o}=\sum_{i\in V}\sum_{j\in V}\,c_{ij}x_{ij}$ (1)

Subject to the following constraints:

1\. Each stop should be visited exactly once:

$\sum_{j\in V,j\neq i}x_{ij}=1\quad\forall i\in V$

2\. Each stop should be left exactly once:

$\sum_{i\in V,i\neq j}x_{ij}=1\quad\forall j\in V$

3\. Elimination of subtours:

$\sum_{i\in S}\sum_{j\in S,j\neq i}x_{ij}\leq|S|-1\quad\forall S\subset
V,|S|\geq 2$

4\. Stops within each cluster should be visited consecutively:

$\sum_{i\in V_{g}}\sum_{j\in V_{g}}x_{ij}=|V_{g}|-1\quad\forall V_{g}\subset
V,|V_{g}|\geq 1,g=1,2,\dots,m$

5\. Binary variable constraints:

$x_{ij}\in\\{0,1\\}\quad\forall i\in V,\forall j\in V,i\neq j$

The objective is to find the values of $x_{ij}$ that minimize the objective
function while satisfying the given constraints.

The 2021 Amazon Last Mile Routing Research Challenge provided a dataset of
real-life operations with over 6000 last-mile delivery route instances. Each
instance consists of a collection of stops that must be visited by a single
driver. Additionally, each instance contains information about a vehicle
driver’s sequence of stops. Each stop has attributes such as latitude and
longitude, delivery time window, and zone ID (cluster to which the stop should
belong). The asymmetric travel times between stops in a route are also
provided on the research challenge dataset introduced by Merchan et al. [17].

In this work, for each instance, we are going to aim to minimize the total
travel time of the route as a first objective. As the second objective, we aim
to predict routes that a driver would prefer to execute. For this second
objective, we are going to evaluate the two approaches that were previously
mentioned. One approach is the use of visual attractiveness metrics [4] to
evaluate and create routes based on the “attractiveness" of the routes to
drivers and decision-makers. The other approach consists in using data mining
strategies to extract historical patterns about the order in which stops were
visited along the routes. The patterns are then used as constraints or
penalizations to guide the route planning. Later in this work, we will
demonstrate with computational experiments which of these two approaches works
better for predicting routes in the order that drivers would prefer to execute
them. In some cases, finding vehicle routes that optimize both objectives
simultaneously will prove to be impossible. Therefore, the problem we tackle
is to approximate the Pareto front [18] for each instance.

##  5 Solution method

In this section, we propose a solution method for the problem involving the
CATSP and the previously mentioned objectives. The solution method is
subdivided into four sections, each addressing a specific aspect of the
problem. Firstly, we present a two-stage algorithm as a general framework for
solving the CATSP, which aims to minimize the total travel time. This
algorithm serves as the foundation for subsequent adaptations. The following
subsection explains how the method is adapted to incorporate visual
attractiveness metrics. We discuss the specific visual attractiveness metrics
used and the composition of the objective function. Additionally, we explore
the adaptation of the method to incorporate patterns found in historical
routes based on ML methods. Finally, we address the bi-objective nature of the
problem. Regardless of whether the second objective is based on visual
attractiveness or ML-based, we aim to find solutions that conform to the
Pareto frontier. Thus, we showcase how the general framework is modified to
handle the bi-objective problem effectively.

###  5.1 Two-stage GRASP algorithm for the CATSP

We start by introducing a multi-start-based metaheuristic as a solver for the
CATSP with the objective of minimizing total route time. The solver is a Two-
Stage Greedy Randomized Adaptive Search Procedure (GRASP) based on the method
introduced by Mesa et al. [12].

The CATSP solver is divided into two main steps to generate route sequences:
ordering the clusters and ordering the stops within the cluster. First, we
compute the average travel times between clusters (zones) and use these times
to solve a TSP of clusters. Then, after the order of the clusters has been
fixed, the stops are randomly assigned to each cluster, and a TSP of all the
stops is solved with heuristics that only allow intra-clusters movements. The
pseudo-code of the CATSP solver is presented in Algorithm 1.

1

Input: $V,T,L,f_{o}$

2

3$cost_{best}\leftarrow\infty$

4$V$ $\leftarrow$ Assign missing zone IDs to stops ($V$)

5$M=\emptyset\leftarrow$Matrix of average cluster-to-cluster (zone-to-zone)
travel times

6for  _$g=1,2,\dots,m$_ do

7

8 for  _$h=1,2,\dots,m\quad h\neq g$_ do

9

10 for  _$i\in V_{g}$_ do

11

12 for  _$j\in V_{h}$_ do

13

14 $M_{gh}=M_{gh}+T_{ij}$

15

16 $M_{gh}=\frac{M_{gh}}{|V_{g}|+|V_{h}|}$

17

18

19for  _$iteration\leq L$_ do

20 Zones solution $z\leftarrow$ Random_order($Z$)

21 $z\leftarrow$ relocate_zones($z,M$)

22 Solution $x\leftarrow$ Randomize ($z,N$)

23 $x\leftarrow$ VND($x,neighborhoods,T$)

24 if  _$f_{o}(x) <cost_{best}$_ then

25 $cost_{best}\leftarrow f_{o}(x)$

26 $x_{best}\leftarrow x$

27

28

29return _$x_{best}$_

Algorithm 1 CATSP Solver. Adapted from [12]

Algorithm 1 describes the proposed Two-Stage GRASP. “The input parameters $V$
and $T$ correspond to the set of stops and the travel times between stops,
respectively. The input parameter $L$ indicates the maximum number of
iterations performed by the GRASP. Line 1 initializes the cost of the best
incumbent solution to infinite. Then in Line 2, when we Assign missing zone
IDs to stops, to each stop that does not have a zone ID, we assign it the zone
ID of the nearest stop (that does have a zone ID) in terms of travel time.
After that, in Lines 3-9, we use the travel time information between the stops
to calculate the estimated travel time $M_{gh}$ between each cluster $g$ to
cluster $h$, which we call zone-to-zone (or cluster-to-cluster) travel time.
Lines 10-17 execute the main loop of the GRASP for $L$ iterations. Line 11 of
the algorithm creates an initial zones solution route $z$, composed of each
cluster (zone) and the depot station $s$. The depot is assigned at the
beginning and the end of route $z$, while all the clusters are randomly
ordered. Next, line 12 applies the classical relocate neighborhood with a
first-improvement strategy, using the cluster-to-cluster travel time
information from $M$ to calculate the cost of the solution. When no further
relocate movements that would improve the cost of $z$ can be applied, the
algorithm proceeds to line 13. The algorithm constructs an initial solution
$x$ of stops $i\in V$, assigning all stops to their corresponding cluster in a
randomized order while keeping the order of the clusters from solution $z$.
After constructing the initial solution $x$, the algorithm applies the local
search phase with a Variable Neighborhood Descent (VND) [26] in line 14. In
lines 15 to 17, the overall best solution is updated when a better solution
(lower objective function value) is found. Finally, line 18 returns the best
solution after all iterations are performed. The VND in this GRASP uses the
following neighborhoods with a best-improvement configuration. Relocate stop
relocates the position of one stop at a time within the same zone. Swap stops
exchanges the position of two stops from the same zone."[12]

###  5.2 Adapting the CATSP solver with visual attractiveness metrics

We can modify the CATSP solver, whose objective function is the minimization
of total travel time, to consider visual attractiveness metrics in an attempt
to plan routes preferred by drivers or decision-makers. In this section, we
explain the visual attractiveness metrics used, the changes to the objective
function, and how the CATSP solver is adapted to incorporate these metrics.

The visual attractiveness metrics are usually formulated in the literature for
VRP variants. Rossit et al. [4] identify three common characteristics that
many authors consider that attractive vehicle routes should have: compactness,
non-overlap, and no-complexity. Different authors have proposed multiple
metrics for each characteristic. However, given that in this work we are
tackling a CATSP, we selected four metrics that can be applied to individual
routes. The visual attractiveness metrics $\upsilon_{k}$ that we use are as
follows.

  * •

Average Distance Compactness introduced by Matis [3]:

$\upsilon_{adc}=\frac{AvgDist}{AvgMaxDist}$ (2)

$AvgDist$ is the average distance (or travel time) between two consecutive
stops in the route, and $AvgMaxDist$ is the average distance of the 20%
longest distances between two consecutive stops in the route. A larger
$\upsilon_{adc}$ indicates that the route is more compact.

  * •

Center Distance Compactness introduced by Kant et al. [27] :

$\upsilon_{cdc}=\sum_{i\in V}{distance(i,c)}$ (3)

$c$ is the stop located at the center (intermediate position) of the route.
This work represents $distance()$ by travel time between stops. A smaller
$\upsilon_{cdc}$ indicates that the route is more compact.

  * •

Number of crossings introduced by Poot et al. [28]:

Indicates the total number of crossings ($\upsilon_{nc}$) between arcs of the
route. This value is an approximation of the real number of crossings, as the
coordinates of the stops are used to determine the crossings geometrically
instead of using a real network of roads. A smaller number of crossings
indicates more visual attractiveness.

  * •

Bending energy introduced by Gretton and Kilby [29]:

$\upsilon_{be}=\frac{\sum_{i=2}^{|V|}{(alpha_{i-2,i-1,i})}}{|V|}$ (4)

We take the vectors formed by the stops $i-2$, $i-1$, and $i$, and compute the
smallest angle $alpha_{i-2,i-1,i}$, in radians, between these vectors. The
smaller the sum of these angles, the more visually attractive.

Each one of these visual attractiveness metrics $\upsilon_{k}$ can be
considered within the objective function when solving the CATSP. We can modify
the objective function from Equation 1 by adding the metric that we want to
consider, as indicated in Equation 5. If the metric that we want to consider
is the Average Distance Compactness $\upsilon_{adc}$, we should subtract it
instead of adding it because, contrary to the other metrics that we want to
minimize, that metric indicates more visual attractiveness when it is
maximized.

$\min f_{o}=\upsilon_{k}+\sum_{i\in V}\sum_{j\in V}\,c_{ij}x_{ij}$ (5)

With the new formulation of the objective function and the means to calculate
the visual attractiveness of a route, we can modify Algorithm 1. We insert the
visual attractiveness metric $\upsilon_{k}$ computation immediately after Line
14. Additionally, we modify lines 15 and 16 to include $\upsilon_{k}$ within
the computation of the objective function $f_{o}(x)$, which now is
$f_{o}(x,\upsilon_{k})$. The pseudo-code of the CATSP solver with visual
attractiveness metrics is presented in the Appendix in Algorithm 5.

###  5.3 Adapting the CATSP solver with data mining from historical routes

Another way in which we can modify the CATSP solver to plan routes preferred
by drivers or decision-makers is to consider information that can be extracted
from historical routes. This section explains the method for extracting
historical patterns from routes, the changes applied to the objective
function, and how the CATSP solver is adapted to incorporate these metrics.

There are multiple approaches that use data mining techniques to find patterns
in delivery routes [6, 9, 10, 11, 12]. For this work, we selected the strategy
proposed by Mesa et al. [12]. The strategy consists in extracting information
from historic routes using a probability estimation method that learns
transition probabilities between clusters and then incorporates these
transition probabilities into route planning algorithms as penalizations on
the objective function. This strategy of data mining (DM) information from
routes works as follows.

We first define a zone centroid for each cluster/zone by taking the average
latitude and longitude coordinates of the stops within that zone. A vector
between each zone centroid and the next zone centroid in a route can be traced
to represent the route. If we do this for every route in a city, we will
construct a vector field across the city using historical routes. However,
each of these vectors might have different magnitudes. Therefore, we divide
them into step vectors, each with a length determined by the parameter
$\beta$. This field of step vectors enables us to estimate a preferred heading
for any arbitrary point on the 2D plane by averaging the observed directions
of vectors in its vicinity.

We then construct a vector field that comprises all step vectors from the
observed sequences of zones in the available historical routes. Let $Q$ be the
set of historic routes. From $Q$, we derive the set of step vectors $V_{j}$
associated with the sequence of zones in each route $j\in Q$. Consequently, we
define the vector field $F$ as the union of all $V_{j}$ sets, i.e.,
$F=\bigcup_{j\in Q}V_{j}$. Each step vector $\mathbf{v}$ in $F$ is
characterized by an origin $o_{v}$, a final position $f_{v}$, and a length of
$\beta$ meters.

Once we have a vector field constructed with the data from all available
historical routes, we can use the information to construct matrices of
transition probability matrices and deviation probability for each new route
that we want to plan. For a given route $x$ to be planned, a transition
probability matrix $H_{ij}$ represents the probability that zone $i\in Z$ is
visited before zone $j\in Z$. The matrix of deviation probability $A_{ij}$
represents the probability that the direction taken by heading from zone $i$
to zone $j$ deviates from the historical direction taken after visiting zone
$i$. These matrices are calculated as follows.

First, we introduce the concept of heading vectors. These vectors point in the
direction where the route sequence should go for each zone centroid of zone
$p\in Z$ after visiting all the stops in the current zone. The estimation of a
heading vector $\mathbf{h}(p)$ for zone $p$ is obtained using Equation (6),
where $\alpha$ is a weighting factor parameter and $D(p,o_{v})$ represents the
distance from zone centroid $p$ to the origin $o_{v}$ of a step vector
$\mathbf{v}$.

$\mathbf{h}(p)=\sum\limits_{\mathbf{v}\in F}\mathbf{v}\cdot
e^{-D(p,o_{v})\cdot\alpha}$ (6)

To efficiently compute the heading vectors, we take advantage of the fact that
the contribution of a step vector $\mathbf{v}$ diminishes as it moves further
away from the zone centroid $p$. This allows us to focus on step vectors
located within a maximum radius Euclidean distance from $p$. We employ the
KDTree algorithm [30], based on a binary tree data structure, to store the
origins of the step vectors. By dividing the vector space in half at each tree
level, the algorithm enables us to find all step vectors $\mathbf{v}$ within
the maximum distance $max_{d}$ (determined such that $e^{max_{d}}<10^{-3}$)
from each zone centroid $p\in Z$ efficiently using a single query.

Once we have computed the heading vectors of a new route, we can estimate the
transition probability matrix $H_{ij}$ using Equation (7). “This equation
takes into account the normalized Euclidean distance $\delta_{ij}$ between
zone centroids $i$ and $j$, as well as the normalized cosine distance
$\epsilon_{ij}$ between the estimated heading vector $\mathbf{h}(i)$ at zone
$i$ and the direction vector from zone $i$ to zone $j$. Both distances are
normalized based on the maximum distances between any pair of zone centroids
within the route. It is important to note that when there is no historical
information available regarding the order between zones $i$ and $j$, we assign
a value of $0.5$ to $H_{ij}$ to avoid any bias in their sequencing."[12]

$H_{ij}=\frac{1}{2}\delta_{ij}+(1-\delta_{ij})\epsilon_{ij}$ (7)

Using the heading vectors, we can also estimate the deviation probability
matrix $A_{ij}$. First, we estimate a direction vector $u_{ij}$ by subtracting
the zone centroid of $i\in Z$ from the zone centroid of $j\in Z$ for each zone
in $Z$. Using the direction vector, we can calculate the deviation angle
$a_{ij}$ (in degrees) between $u_{ij}$ and the heading vector originating from
the zone centroid of $i$. This deviation angle is used in Equation (8) to
construct the deviation probability matrix. Notably, when $A_{ij}=1$, this
indicates that the direction the route takes after visiting zone $i$ is
probably deviating from the historical route order. Conversely, when
$A_{ij}=0$, the route follows the estimated historical direction.

$A_{ij}=\frac{a_{ij}}{180}$ (8)

We can also leverage data from the zone centroids that we construct for the
new route. We can estimate the travel time $r_{ij}$ between the zone centroids
of $i\in Z$ and $j\in Z$ using Equation (9). In this equation, $N_{i}$ and
$N_{j}$ represent the sets of stops associated with zone $i$ and $j$
respectively, and $t_{kl}$ denotes the travel time between stop $k\in N_{i}$
and stop $l\in N_{j}$. The zone-to-zone travel time $M_{ij}$ is calculated as
the sum of all travel times between stops in $N_{i}$ and $N_{j}$, divided by
the total number of stops in both sets. This calculation can easily be
extended to include the estimated travel time from each zone to the delivery
station $s$, assuming $s$ to represent a zone with only one stop.

$M_{ij}=\frac{\sum\limits_{k\in N_{i}}\sum\limits_{l\in
N_{j}}t_{kl}}{|N_{i}|+|N_{j}|}$ (9)

We now have the matrices $M_{ij}$, $H_{ij}$, and $A_{ij}$ that contain
information about route $x$ and from the historical routes that were executed
nearby route $x$. We incorporate this information to plan a route as follows.

Using the route $x$, we can estimate the binary variable $m_{ij}$, which is
equal to 1 if zone $i\in Z$ immediately precedes zone $j\in Z$, and 0
otherwise. We can estimate the deviation from the historical angular direction
$\phi$ using Equation (11). This equation computes $\phi$ by summing the
product of $m_{ij}$ and $A_{ij}$ for all pairs $i,j\in N$, and normalizing it
by the total number of zones $|Z|$. The resulting value of $\phi$ provides an
estimation of the extent to which the route’s actual order deviates from the
expected historical direction. We can also calculate the value of the
deviation from historical heading vectors $\eta$ as defined in Equation (10).
This equation computes $\eta$ by summing the product of $w_{ij}$ and
$(1-H_{ij})$ for all pairs $i,j\in Z$, and normalizing it by the total number
of zones $|N|$. Finally, we can estimate the average zone-to-zone travel time
of the route $\lambda$. Equation (12) provides an estimation of $\lambda$ by
summing the product of $m_{ij}$ and $M_{ij}$ for all pairs $i,j\in Z$, and
dividing it by the maximum travel time between any two zones. The resulting
value of $\lambda$ captures the overall impact of travel times on the route,
considering both inter-zone travel and travel between zones and the delivery
station.

$\eta=\frac{\sum\limits_{i,j\in Z}m_{ij}\cdot(1-H_{ij})}{|N|}$ (10)

$\phi=\frac{\sum\limits_{i,j\in N}m_{ij}\cdot A_{ij}}{|Z|}$ (11)

$\lambda=\frac{\sum\limits_{i,j\in Z}m_{ij}\cdot M_{ij}}{max(M_{ij})}$ (12)

We can include this information within a routing optimization algorithm by
modifying the objective function from Equation 1 as follows. First, we
normalize the total travel time and add it as the component normalized total
time ($\tau$). Then, we add the components: historic headings deviation
($\eta$), historic angles deviation ($\phi$), and zones-to-zones travel time
($\lambda$). Each component has different magnitudes and is normalized to
ensure comparable values. We assign a weighting factor $\theta_{i}$ to each
component, reflecting its importance in the objective function. The objective
function that considers travel time and information from historical routes is
shown in Equation 13.

$f_{o}=\theta_{1}\tau+\theta_{2}\eta+\theta_{3}\phi+\theta_{4}\lambda$ (13)

We modify the single-objective CATSP solver1 following the Two-stage data-
driven GRASP metaheuristic introduced by Mesa et al. [12]. The modified CATSP
solver uses two phases (routines) to generate route sequences: the historical
information extraction phase and the route planning phase. First, it extracts
information from historical routes in the form of transition probability
matrices. Then, it proceeds with the CATSP solver, using traditional routing
information (travel time) and adding the transition probability matrices to
plan the routes. The pseudo-code of the modified CATSP solver with historical
routing information is presented in the Appendix in Algorithm 6. For a
complete explanation of the method, we refer the reader to the work of Mesa et
al. [12].

1

Input: $V,T,Q$

2

3$F$ $\leftarrow$ Create vector field ($Q$)

4$H_{ij}$, $A_{ij}$ $\leftarrow$ Calculate transition probability matrices
($F,V$)

5$M_{ij}$ $\leftarrow$ Calculate zone-to-zone travel-time ($V$)

Algorithm 2 Extract historical information phase. Adapted from [12]

The first stage of the approach is the extract historical information phase.
The Create vector field step constructs the vector field $F$ using historical
route information from the set of historic routes $Q$. The Calculate
transition probability matrices step leverages $F$ and its set of zones $Z$ to
calculate the transition probability matrices $H_{ij}$ and $A_{ij}$ for a new
route $x$ that needs planning. We condense lines 4 to 9 from Algorithm1 into
the Calculate zone-to-zone travel-time method. Moving on to the second stage,
the route prediction phase, the routing process of Algorithm1 is executed,
with the addition of $H_{ij}$, $A_{ij}$, and $M_{ij}$ into the objective
function $f_{o}$ to generate the new route sequence $x$.

In summary, we can construct a vector field by analyzing historical routes
using the step vectors between consecutive zone centroids. This vector field
captures the preferred directions of routes across the city, allowing us to
estimate the heading or order of zones for new routes based on their proximity
to the observed routes in the vector field. Finally, this information, along
with stop-to-stop travel times and zone-to-zone estimated travel times, can be
fed into a weighted objective function that can be used within a CASTP solver.

###  5.4 Adapting the CATSP solver for solving the bi-objective routing
problem

Independently of whether the method that uses visual attractiveness metrics or
the one that leverages data mining information from historical routes works
better, we are going to tackle the problem as a bi-objective problem. To this
end, we aim to find solutions that make up the Pareto front. In this section,
we show how we adapt the general framework of the CATSP solver in order to
solve a bi-objective vehicle routing problem.

Let’s assume we have two objective functions $f_{1}$ and $f_{2}$ for the
CATSP, $f_{1}$ considers traveling time cost, while $f_{2}$ considers
preferences from the drivers or decision makers. To solve this bi-objective
problem, we integrate the single-objective CATSP solver from Algorithm 1
within a heuristic box splitting algorithm [31] to enable the ability to
generate a set of Pareto-optimal solutions (just like any multi-objective
algorithm).

The heuristic box splitting (HBS) algorithm proposed by Matl et al. [31] is
based on representing the limits of the objective space with an initial
rectangle and iteratively solving and splitting the rectangle space into
smaller rectangles while respecting the limits imposed on the objective
function by the current rectangle. HBS can be applied to a route optimization
algorithm. Therefore, we take the GRASP algorithm in Algorithm 1, which
corresponds to lines 10 to 17, and use it as the route optimization heuristic
with some modifications. The modifications consist of using the objective
functions $f_{1}$ and $f_{2}$, by using $f_{1}$ as the objective function to
minimize subject to $f_{2}\leq c$ as a constraint. We present the Route
Optimization Heuristic (ROH) in Algorithm 3.

1

Input: $N,Z,L,c$

2$cost_{best}\leftarrow\infty$

3for  _$i\leq L$_ do

4 Zones solution $z\leftarrow$ Random_order($Z$)

5 $z\leftarrow$ relocate_zones($z,c$)

6 Solution $x\leftarrow$ Random_order ($z,N$)

7 $x\leftarrow$ VND($x,neighborhoods,c$)

8 if  _$f_{1}(x) <cost_{best}$ and $f_{2}(x)\leq c$_ then

9 $cost_{best}\leftarrow f_{1}(x)$

10 $x_{best}\leftarrow x$

11

12

13return _$x_{best}$_

Algorithm 3 Route optimization heuristic (ROH).

Now, we take ROH and the CATSP solver, and integrate them into the heuristic
box algorithm. The pseudocode of the proposed bi-objective CATSP solver is
presented in Algorithm 4. The proposed algorithm uses rectangles or boxes to
represent areas in the objective space that may contain non-dominated
solutions. The algorithm maintains an archive of these rectangles $R_{areas}$
and a set of non-dominated solutions $A_{sol}$ that have been found.
Initially, $f_{2}^{min}$ is a valid lower bound on objective $f_{2}$ (We use
$f_{2}^{min}=0$), $f_{1}^{max}$ is any upper bound on objective $f_{1}$ (we
use $f_{1}^{max}=cost$(random nearest neighbor solution)), and $\delta_{min}$
is the minimal difference in objective $f_{2}$ when considering two feasible
solutions (we use $\delta_{min}=0.01$. The other inputs of the algorithm
($V,T,L$) are the same inputs of Algorithm 1.

The algorithm starts by Assigning missing zone IDs, where the missing zone IDs
of the stops in a given route are assigned. If a stop does not have a zone ID,
the algorithm assigns it the zone ID of the nearest stop (that already has a
zone ID) based on the estimated travel time between them. Then, in Calculate
zone-to-zone travel-time, the algorithm utilizes the travel time information
between the stops in a given route to calculate the estimated travel time
between zones $M_{ij}$. After calculating the zone centroids travel-time, the
algorithm proceeds with the phase of Routes prediction phase: Heuristic Box
Splitting Two-Stage GRASP. First, we execute the ROH of Algorithm 3 to attempt
to minimize objective $f_{1}$ without assigning an upper bound to objective
$f_{2}$ (which is the same as subjecting it to the constraint
$f_{2}\leq\infty$). This yields the “upper-left” extreme point ($f_{1}^{min}$,
$f_{2}^{max}$). The “initial lower-right” extreme point corresponds to the two
initial bounds ($f_{1}^{max}$, $f_{2}^{min}$). Together, these two points
delimit the initial box (or rectangle). The “true lower-right" point is
progressively approximated by splitting unexamined regions of the objective
space and discarding redundant regions. The proposed algorithm iteratively
splits the largest remaining rectangle in half by setting the $\epsilon-
constraint$ and then solves the corresponding subproblem using ROH. The
boundaries of the selected rectangle determine the value of $c$, which
indicates the upper bound of ($f_{2}\leq c$) when solving ROH. The subproblem
may be infeasible, so ROH returns NULL instead of a feasible solution. In such
cases, the archive $A_{sol}$ remains unchanged, and the explored half of the
rectangle is discarded, resulting in no new rectangle. If a solution $x$ is
found, this resulting solution lies within the boundaries of the selected
rectangle, and it will dominate a region of the rectangle. This region is
removed, $R_{areas}$ is updated with the new rectangles, and $A_{sol}$ is
updated with the new non-dominated solution while solutions that are dominated
by $x$ are removed from it. The process is repeated until $R_{areas}$ is empty
or after $n_{max}$ iterations are performed. Finally, the algorithm’s output
is the archive $A_{sol}$ with all the non-dominated solutions found during the
execution.

1

Input: $V,T,L,n_{max}$, $f_{1}^{max}$, $f_{2}^{min}$, $\delta^{min}$

Output: Archive of non-dominated solutions $A_{sol}$

2$V$ $\leftarrow$ Assign missing zone IDs ($V$)

3$M_{ij}$ $\leftarrow$ Calculate zone-to-zone travel-time ($V$)

4$x\leftarrow ROH(f_{1},\infty)$

5$A_{sol}\leftarrow\\{x\\}$ , $z^{1}\leftarrow(x_{1},x_{2})$ ,
$z^{2}\leftarrow(f_{1}^{max},f_{2}^{min})$,
$R_{areas}\leftarrow\\{$R$(z^{1},z^{2})\\}$

6for  _$i\leq n_{max}$_ do

7 Find R$(y^{1},y^{2})\in R_{areas}$ with maximal area

8 $c\leftarrow\frac{1}{2}(y_{2}^{1}+y_{2}^{2})$

9 $x\leftarrow ROH(f_{1},c)$

10 if  _$x=NULL$ or $x$ is dominated_ then

11

12 $y^{2}\leftarrow(y_{1}^{2},c)$

13 else

14 Update archive $A_{sol}$ with $x$

15 if  _there exists{R’$(z^{1},z^{2})\in R_{areas}$ $|z_{1}^{1}<x_{1}\leq
z_{1}^{2}$}_ then

16

17 $z*\leftarrow(x_{1},max\\{z_{2}^{2},c\\})$

18 if  _ $z_{2}^{1}-z_{2}^{*}>\delta_{min}$_ then

19

20 Add R$(z^{1},z^{*})$ to $R_{areas}$

21

22

23 if  _there exists{R”$(z^{1},z^{2})\in R_{areas}$ $|z_{2}^{1}\geq x_{2}\geq
z_{2}^{2}$}_ then

24

25 $z^{*}\leftarrow(max\\{x_{1},z_{1}^{1}\\},x_{2}$)

26 if  _$z_{2}^{*}-z_{2}^{2} >\delta_{min}$_ then

27

28 Add $R(z*,z2)$ to $R_{areas}$

29

30

31 Delete R’ and R” if they exist

32 forall  _R $(z^{1},z^{2})\in R_{areas}$ with $z^{1}$ dominated by $x$_ do

33

34 Delete R from $R_{areas}$

35

36 If $R_{areas}=\emptyset$

37 Stop

38

39 return _$A_{sol}$_

Algorithm 4 Bi-objective Two-Stage GRASP with heuristic box splitting. Based
on [31, 12].

$f_{1}=\theta_{1}\tau$ (14)

$f_{2}=\theta_{2}\eta+\theta_{3}\phi+\theta_{4}\lambda$ (15)

Later, in Section 6, we explain in detail the experiments and the results of
applying these three approaches (Visual attractiveness metrics, Data mining
from historical routes, and Bi-objective routing) to last-mile delivery
planning. In summary, the single objective historical pattern mining strategy
works better than the single objective visual attractiveness metrics.
Therefore, we take the method for extracting historical patterns from routes
and integrate it within the bi-objective algorithm, as shown in the Appendix
in Algorithm 7. Additionally, we construct the following objective functions
that formulate the bi-objective mathematical problem. Equation (14) seeks to
optimize the total routing cost. Equation (15) seeks to minimize the weighted
differences in the routing sequence with respect to historical routes.

##  6 Computational experiments

In this section, we present and analyze the results of two sets of experiments
carried out in our work. First, we compare the use of visual attractiveness
metrics vs. leveraging historical information from routes to predict last-mile
delivery routes. Then, taking the method with better performance in the first
experiment, we enhance it to enable the solution of bi-objective problems and
evaluate the algorithm’s ability to approximate the Pareto front when solving
the proposed bi-objective model. For both experiments, we use the set of CATSP
instances provided by Amazon and MIT in the 2021 Amazon Last Mile Routing
Research Challenge [13, 17].

All the routing algorithms were coded in C++ and compiled using g++ 11.2.0 and
Boost 1.78 library [32]. As previously mentioned, to estimate the vector
fields, the heading vectors, and the different probability matrices, we used
the code provided by Mesa et al. [12] 111Available at
https://github.com/mesax1/Data-driven-last-mile-routing, which is written in
Python 3.8. All the experiments were run on an Intel Xeon Gold 6130 CPU,
running at 2.10 GHz, with 64 GB of RAM, on a Rocky Linux 8.5 Operating System.
All experiments were run on a single CPU core. The code repository for this
work is available at https://github.com/mesax1/biobjective_lastmile_routing.

###  6.1 Experiments of visual attractiveness metrics versus historical
patterns mining

We evaluate and compare two approaches to create routes as similar as possible
to the way in which routes have been previously executed. This is with the
objective of minimizing the differences between previously executed routes and
newly planned routes. The two approaches are: using visual attractiveness
metrics to assess the routes and data mining techniques to find patterns in
the routes.

To evaluate which of the two approaches works best, we ran a set of
experiments that determines which strategy generates routes that resemble the
most previously executed routes. To this end, we used Amazon’s dataset of
last-mile delivery routes [17] and a score similarity metric that compares two
routes and outputs a scalar value that indicates their similarity.

The score similarity metric, in Equation 16, calculates the similarity between
two route sequences, $x$ and $B$, based on three factors: the sequence
deviation $SD(x,B)$, the normalized edit distance $ERP_{norm}(x,B)$, and the
number of edit operations $ERP_{e}(x,B)$ required to transform $x$ into $B$.
If $x$ and $B$ are the same sequence, then $score(x,B)$ is 0. Additional
information on how the score similarity metric is computed is given in the
Appendix 19.

$score(x,B)=\frac{SD(x,B)\cdot ERP_{norm}(x,B)}{ERP_{e}(x,B)}$ (16)

The sequence deviation $SD(x,B)$, shown in Equation 17, measures the deviation
of the position of stops in the predicted route sequence $x$ with respect to
the actual route sequence $B$. Assuming both routes visit the same $n+1$ stops
($n$ stops and one depot station). First, we identify the position $g$ of each
stop in sequence $x$ within sequence $B$. We then form a vector
$(g_{0},g_{1},...,g_{n})$. Next, we create vector $(g_{0},g_{1},...,g_{n})$.
$SD(x,B)$ is the sum of the absolute differences between consecutive positions
in this vector, divided by the total number of unique stop combinations.

$SD(x,B)=\frac{2}{n(n-1)}\sum_{i=1}^{n}(|g_{i}-g_{i-1}|-1)$ (17)

“The number of edit operations $ERP_{e}(x,B)$ required to transform sequence
$x$ into sequence $B$ is calculated using the Levenshtein distance, which is a
measure of the difference between two sequences introduced by Levenshtein et
al. [33]. Let $E=E_{0},E_{1},...,E_{n}$ be the edit operations to transform
sequence $x$ into sequence $B$. Additionally, let $C=c_{0},c_{1},...,c_{n}$ be
the normalized travel time costs associated with the operations performed.
Each operation of insertion, deletion, and substitution has a unit cost."[12].
Then, the normalized edit distance is calculated as shown in Equation 18.

$ERP_{norm}(x,B)=\sum_{i=0}^{n}c_{i}E_{i}$ (18)

A lower score indicates that $x$ is more similar to $B$. A score of zero means
they are exactly the same sequence. This metric was utilized in the 2021
Amazon Last Mile Routing Research Challenge [13]222Further details of the
scoring algorithm and related implementations can be found at
https://github.com/MIT-CAVE/rc-cli/tree/main/scoring.

To evaluate the performance of routing guided by historical patterns found in
previous routes, we modified Algorithm 1 following the procedure described in
Section 5.3, using the same values for the parameters and the objective
function weights (Equation 13) as proposed by Mesa et al. [12]. The method
consists of four weighting factors $(\theta_{i})$ and two parameters that need
to be calibrated based on the available data. The four components of the
objective function 13 are assigned a weighting factor $\theta_{i}$. In
addition, there are two parameters to consider: the length of step vectors
($\beta$) in meters for estimating heading vectors and the weighting factor
($\alpha$) for step vectors. Table 1 outlines the values for each parameter.

Parameter Description Value $L$ Number of GRASP iterations 50 $\alpha$
Weighting factor for step vectors 0.01 $\beta$ Length of the step vectors in
meters 200 $\theta_{1}$ Weighting factor of total time $\tau$ 3 $\theta_{2}$
Weighting factor of historic zone headings $\eta$ 1 $\theta_{3}$ Weighting
factor of historic zone angles $\phi$ 1 $\theta_{4}$ Weighting factor of zones
time $\lambda$ 5

Table 1:  Parameters of the data mining from historical routes method. Adapted
from [12].

For evaluating the performance of visual attractiveness metrics, and ensuring
a fair comparison, we use the same code provided by Mesa et al.[12] and modify
it following the procedure described in Section 5.2. We incorporate the four
visual attractiveness metrics: Average Distance Compactness ($adc$), Center
Distance Compactness ($cdc$), Number of Crossings ($nc$), and Bending Energy
($be$), as penalization components of the objective function. We keep the
parameter $L$ with the same value. Table 2 presents the values for each
parameter.

Parameter Description Value $L$ Number of GRASP iterations 50 $\theta_{1}$
Weighting factor of total time $\tau$ 3 $\theta_{4}$ Weighting factor of zones
time $\lambda$ 5 $\theta_{6}$ Weighting factor of average distance compactness
$\upsilon_{adc}$ 1 $\theta_{7}$ Weighting factor of center distance
compactness $\upsilon_{cdc}$ 1 $\theta_{8}$ Weighting factor of number of
crossings $\upsilon_{nc}$ 1 $\theta_{9}$ Weighting factor of bending energy
$\upsilon_{be}$ 1

Table 2:  Parameters of the visual attractiveness metrics method.

We use the 3052 routes of the evaluation dataset to evaluate both methods. For
the historical patterns data mining method, we execute Algorithm 6 once per
route, with the parameter configuration presented in Table 1, considering all
the components of the objective function. We then calculate the score of each
predicted route, by comparing it against how the real route was executed,
using the scoring metric function of Equation 16. For the visual
attractiveness metrics method, we execute Algorithm 5, with the parameters
configuration of Table 2. For each visual attractiveness metric, we execute
the algorithm once per route, keeping the components of total time $\tau$
($\theta_{1}=3$) and zones time $\lambda$ ($\theta_{4}=5$), and including only
the single attractiveness metric that we want to evaluate with a weight of 1,
e.g. the number of crossings $\upsilon_{nc}$ ($\theta_{8}=1$). Then, we
calculate the score for each route using the same procedure as for the
historical patterns method.

We summarize the score results obtained with these experiments in Figure 2 and
Table 3. The results for the historical patterns data mining method are
represented with the name “DM". The results for each visual attractiveness
metric are represented by the name of the parameter, e.g. $\upsilon_{adc}$.
Additionally, we include the results of only using the components of total
time ($\theta_{1}=3$) and zones time ($\theta_{5}=5$) in the objective
function as a baseline, represented with the name “BASE". In Table 3, for each
column, we report the average route score “Mean", the standard deviation
“Std", the best score of a route “Min", the “Median" score, and the worst
score of a route “Max".

Figure 2: Distributions of the scores obtained with different routing methods
on the evaluation dataset.

BASE DM $\upsilon_{adc}$ $\upsilon_{cdc}$ $\upsilon_{nc}$ $\upsilon_{be}$ Mean
0.0484 0.0413 0.0822 0.0817 0.0848 0.1076 Std 0.0335 0.0355 0.0544 0.0540
0.0561 0.0531 Min 0.0025 0.0008 0.0017 0.0020 0.0027 0.0117 Median 0.0428
0.0337 0.0696 0.0688 0.0716 0.0989 Max 0.3543 0.3249 0.4091 0.4036 0.4041
0.5942

Table 3:  Results for scores obtained with different routing methods on the
evaluation dataset.

Visually, we can notice the performance difference between the methods. The
best performance is achieved using the historical patterns data mining method,
followed by the “classical" baseline, which aims to minimize travel times,
both between stops and between zones. When any single visual attractiveness
metric is added to the baseline objective function, the performance worsens.
The average score for the baseline is 0.0484. The “DM" method improves the
score to 0.0411, while the best average score achieved by a visual
attractiveness metric is just of 0.0817 for the $\upsilon_{cdc}$, and the
worst average score of 0.1076 is obtained with the $\upsilon_{be}$ metric. For
this particular situation of predicting single-vehicle last-mile delivery
routes, the advantage of using the historical patterns data mining method is
clear. Therefore, we select this method and use it to address the bi-objective
vehicle routing problem proposed in this work.

###  6.2 Experiments of bi-objective vehicle routing problem

After selecting the method for solving the single-objective VRP problem, we
modified the method following the heuristic box-splitting strategy proposed by
Matl et al. [31]. Now, instead of having the single objective function of
Equation 13, we have the bi-objective problem with Equation 14 and Equation
15. We followed the procedure described in Section 5.4, and turned Algorithm 6
into Algorithm 7. Table 4 presents the parameters of the bi-objective CATSP
solver and their corresponding values.

Parameter Description Value $L$ Number of GRASP iterations 50 $\alpha$
Weighting factor for step vectors 0.01 $\beta$ Length of the step vectors in
meters 200 $\theta_{1}$ Weighting factor of total time $\tau$ 1 $\theta_{2}$
Weighting factor of historic zone headings $\eta$ 1 $\theta_{3}$ Weighting
factor of historic zone angles $\phi$ 1 $\theta_{4}$ Weighting factor of zones
time $\lambda$ 5 $n_{max}$ Maximum number of HBS iterations 50 $f_{1}^{max}$
Initial upper bound of $f_{1}$ (24 hours in seconds) 86400 $f_{2}^{min}$
Initial lower bound of $f_{2}$ 0

Table 4:  Parameters of the bi-objective CATSP solver.

We run Algorithm 7 on the 3052 instances of the evaluation dataset. Due to the
randomness involved in the algorithm, ten runs per instance were performed.
The average running time for an instance is 45 seconds. Table 5 presents the
results on the average number of non-dominated solutions obtained with our bi-
objective CATSP solver. The first column represents the number of non-
dominated solutions found. The second column shows the average number of
instances where the non-dominated solutions were found.

Number of non-dominated solutions found Number of instances 1 604 2 1020 3 769
4 397 5 or more 262

Table 5:  Number of non-dominated solutions found.

Description Value Average maximum travel time 11812 s Average minimum travel
time 11251 s Average difference between maximum and minimum travel time 602 s
Number of routes with more than 10 minutes maximum time difference 1202 Number
of routes with more than 15 minutes maximum time difference 569 Number of
routes with more than 30 minutes maximum time difference 47 Number of routes
with more than 45 minutes maximum time difference 5 Number of routes with more
than 60 minutes maximum time difference 1

Table 6:  Travel time differences in the Pareto front.

From Table 5, it can be seen that the algorithm obtains at least 2 non-
dominated solutions for 80.21% (2448 instances out of 3052) of the instances.
Only for 8.58% of the instances the algorithm found 5 or more non-dominated
solutions. On average, there are at least 2.60 non-dominated solutions for
each instance. Three factors can explain the latter. The first one is related
to how the algorithm explores the objective space: Even if the algorithm tries
to minimize the objective functions $f_{1}$ and $f_{2}$, the local search of
the algorithm is guided by the travel time cost, both between zones (clusters)
and from customer to customer, where solutions with lower travel time cost are
mainly explored. The second factor relates to how the historical routes are
planned. Due to tight operations margins, last-mile delivery routes are
typically planned in a way that minimizes operations costs. Therefore,
information from historical routes will also trend toward minimizing travel
time. Finally, as indicated by Mesa et al. [12], almost 37% of the routes are
located in areas where no route in the training dataset was executed, so there
is not enough information available for the objective function $f_{2}$ to
consider.

The approximated Pareto fronts for multiple instances are plotted in Figures 3
through 6. The approximated Pareto fronts obtained by the proposed bi-
objective CATSP solver are plotted in blue. Additionally, in Table 6 we report
information around the average travel time differences for the 2448 instances
where at least two non-dominated solutions were found. When comparing the
travel time cost for an instance, the average difference between the possible
route found with the minimum travel time and the most similar route to the
historical routes is 10 minutes. This relatively small difference can be
explained by the way Algorithm 7 focuses on searching routes and improvements
based on reducing the traveling time cost, therefore most solutions will tend
to minimize travel time.

Figure 3: Approximated Pareto front for instance 9. Figure 4: Approximated
Pareto front for instance 31. Figure 5: Approximated Pareto front for instance
184. Figure 6: Approximated Pareto front for instance 218.

The obtained results show that the proposed algorithm is able to find non-
dominated solutions. As a result, a Pareto front can be approximated. This is
important because having a set of non-dominated solutions instead of a unique
solution gives support to the decision-maker. For example, a decision-maker
may want to maximize historical route-similarity for new drivers to promote
learning familiar routes, but bearing in mind the total travel time to avoid
long routes that exceed working time limits.

##  7 Conclusions

In this work, we introduced a bi-objective vehicle routing problem with
simultaneous minimization of traveling costs and differences with historical
patterns. Additionally, we compared the strategies of using visual
attractiveness metrics or data mining of patterns, to predict routes that
resemble how routes are executed by drivers. These situations were studied
using data from a real-world last-mile delivery operation. The main objective
of this work is to help route planners and decision-makers to find the trade-
off between travel time cost and learned patterns. We introduce a bi-objective
metaheuristic that uses a simple GRASP and a heuristic box-splitting procedure
to approximate a Pareto front for each instance.

Our results when comparing visual attractiveness metrics versus data mining of
historical patterns, suggest that the latter strategy is better to predict
routes in the way that drivers execute them, compared to the visual
attractiveness metrics found in the literature. The results of the
computational experiments indicate that the algorithm is able to find only a
small amount of non-dominated solutions per instance. The small number of non-
dominated solutions generated per instance can be beneficial to decision-
makers. Using the set of solutions provided by the Pareto front, a decision-
maker can choose the routes that the drivers should follow. Selecting a route
amongst two or three different options can be simpler and more time-efficient
compared to selecting a route amongst 20 or 30 options. This becomes more
important when the operation requires hundreds or thousands of routes that
need to be dispatched daily, and the time spent deciding which route to
execute for each instance also has to be minimized.

Our proposed method uses a simple GRASP as the main route optimization
procedure. As a future research direction, this GRASP can be replaced by a
state-of-the-art routing algorithm that can consider penalizations within its
objective function. This could lead to a broader and better exploration of the
solution space, which can lead to finding more solutions on the Pareto front.
Another idea for further research is to consider other methods for extracting
information from historical routes. For example, some other ideas presented in
[11] could be considered to guide the local search heuristic via constraints,
as Cook et al. [25] further explain in their work.

##  8 Acknowledgments

The authors acknowledge supercomputing resources made available by the
Universidad EAFIT scientific computing center (APOLO) to conduct the research
reported in this work and for its support for the computational experiments.
Finally, the first author is grateful to Universidad EAFIT for the Ph.D. Grant
from project 974-000005.

##  9 Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the
publication of this paper.

## References

  * Toth and Vigo [2014] Toth P, Vigo D.  Vehicle routing: problems, methods, and applications.  Vigo D, Toth P, editors, Philadelphia, PA: Society for Industrial and Applied Mathematics; 2014.  https://epubs.siam.org/doi/abs/10.1137/1.9781611973594.
  * Tang and Miller-Hooks [2006] Tang H, Miller-Hooks E.  Interactive Heuristic for Practical Vehicle Routing Problem with Solution Shape Constraints.  Transportation Research Record 2006;1964(1):9–18.  https://doi.org/10.1177/0361198106196400102.
  * Matis [2008] Matis P.  Decision support system for solving the street routing problem.  Transport 2008;23(3):230–235.
  * Rossit et al. [2019] Rossit DG, Vigo D, Tohmé F, Frutos M.  Visual attractiveness in routing problems: A review.  Computers & Operations Research 2019;103:13–34.
  * Rocha et al. [2022] Rocha D, Aloise D, Aloise DJ, Contardo C.  Visual attractiveness in vehicle routing via bi-objective optimization.  Computers & Operations Research 2022;137:105507.
  * Krumm [2008] Krumm J.  A markov model for driver turn prediction.  In: Society of Automotive Engineers (SAE) 2008 World Congress, April 2008; 2008. .
  * Toledo et al. [2013] Toledo T, Sun Y, Rosa K, Ben-Akiva M, Flanagan K, Sanchez R, et al.  Decision-making process and factors affecting truck routing.  In: Freight Transport Modelling Emerald Group Publishing Limited; 2013.
  * Wang et al. [2015] Wang X, Ma Y, Di J, Murphey YL, Qiu S, Kristinsson J, et al.  Building efficient probability transition matrix using machine learning from big data for personalized route prediction.  Procedia Computer Science 2015;53:284–291.
  * Delling et al. [2015] Delling D, Goldberg AV, Goldszmidt M, Krumm J, Talwar K, Werneck RF.  Navigation made personal: Inferring driving preferences from gps traces.  In: Proceedings of the 23rd SIGSPATIAL international conference on advances in geographic information systems; 2015. p. 1–9.
  * Canoy and Guns [2019] Canoy R, Guns T.  Vehicle routing by learning from historical solutions.  In: International Conference on Principles and Practice of Constraint Programming Springer; 2019. p. 54–70.
  * Winkenbach et al. [2021] Winkenbach M, Parks S, Noszek J.  Technical Proceedings of the Amazon Last Mile Routing Research Challenge.  MIT Libraries; 2021. https://hdl.handle.net/1721.1/131235.
  * Mesa et al. [2023] Mesa JP, Montoya A, Ramos-Pollán R, Toro M.  A two-stage data-driven metaheuristic to predict last-mile delivery route sequences.  Engineering Applications of Artificial Intelligence 2023;125:106653.  https://www.sciencedirect.com/science/article/pii/S0952197623008370.
  * ama [2021] Amazon Last Mile Routing Research Challenge; 2021.  https://routingchallenge.mit.edu/.
  * Mesa et al. [2021] Mesa JP, Montoya A, Ramos-Pollán R.  A Data-Driven Metaheuristic Approach to Predict Delivery Route Sequences.  In: Winkenbach M, Parks S, Noszek J, editors. Technical Proceedings of the Amazon Last Mile Routing Research Challenge MIT Libraries; 2021. p. VIII.1–VIII.12.
  * Cook et al. [2021] Cook W, Held S, Helsgaun K.  Local Search with Learned Constraints for Last Mile Routing.  In: Winkenbach M, Parks S, Noszek J, editors. Technical Proceedings of the Amazon Last Mile Routing Research Challenge MIT Libraries; 2021. p. XXI.1–XXI.12.
  * Canoy et al. [2021] Canoy R, Mandi J, Bucarey V.  TSP with Learned Zone Preferences for Last-Mile Vehicle Dispatching.  In: Winkenbach M, Parks S, Noszek J, editors. Technical Proceedings of the Amazon Last Mile Routing Research Challenge MIT Libraries; 2021. p. XXX.1–XXX.9.
  * Merchan et al. [2022] Merchan D, Arora J, Pachon J, Konduri K, Winkenbach M, Parks S, et al.  2021 Amazon last mile routing research challenge: Data set.  Transportation Science 2022;https://www.amazon.science/publications/2021-amazon-last-mile-routing-research-challenge-data-set.
  * Pareto [1964] Pareto V.  Cours d’économie politique, vol. 1.  Librairie Droz; 1964.
  * Li and Phillips [2018] Li Y, Phillips W, Learning from route plan deviation in last-mile delivery; 2018.  https://ctl.mit.edu/pub/thesis/learning-route-plan-deviation-last-mile-delivery.
  * Neto et al. [2018] Neto FDN, de Souza Baptista C, Campelo CE.  Combining Markov model and prediction by partial matching compression technique for route and destination prediction.  Knowledge-Based Systems 2018;154:81–92.
  * Guo et al. [2021] Guo X, Mo B, Wang Q.  Last-Mile Delivery Trajectory Prediction Using Hierarchical TSP with Customized Cost Matrix.  In: Winkenbach M, Parks S, Noszek J, editors. Technical Proceedings of the Amazon Last Mile Routing Research Challenge MIT Libraries; 2021. p. II.1–II.12.
  * Arslan and Abay [2021] Arslan O, Abay R.  Data-Driven Vehicle Routing in Last-Mile Delivery.  In: Winkenbach M, Parks S, Noszek J, editors. Technical Proceedings of the Amazon Last Mile Routing Research Challenge MIT Libraries; 2021. p. VI.1–VI.14.
  * Jozefowiez et al. [2008] Jozefowiez N, Semet F, Talbi EG.  Multi-objective vehicle routing problems.  European journal of operational research 2008;189(2):293–309.
  * Lian et al. [2016] Lian K, Milburn AB, Rardin RL.  An improved multi-directional local search algorithm for the multi-objective consistent vehicle routing problem.  IIE Transactions 2016;48(10):975–992.
  * Cook et al. [2022] Cook W, Held S, Helsgaun K.  Constrained local search for last-mile routing.  Transportation Science 2022;.
  * Mladenović and Hansen [1997] Mladenović N, Hansen P.  Variable neighborhood search.  Computers & operations research 1997;24(11):1097–1100.
  * Kant et al. [2008] Kant G, Jacks M, Aantjes C.  Coca-Cola Enterprises optimizes vehicle routes for efficient product delivery.  Interfaces 2008;38(1):40–50.
  * Poot et al. [2002] Poot A, Kant G, Wagelmans APM.  A savings based method for real-life vehicle routing problems.  Journal of the Operational Research Society 2002;53(1):57–68.
  * Gretton and Kilby [2013] Gretton C, Kilby P.  A study of shape penalties in vehicle routing.  In: VIII Triennial Symposium on Transportation Analysis. San Pedro de Atacama, Chile; 2013. .
  * Bentley [1975] Bentley JL.  Multidimensional Binary Search Trees Used for Associative Searching.  Commun ACM 1975 Sep;18(9):509–517.
  * Matl et al. [2019] Matl P, Hartl RF, Vidal T.  Leveraging single-objective heuristics to solve bi-objective problems: heuristic box splitting and its application to vehicle routing.  Networks 2019;73(4):382–400.
  * Boo [2022] Boost C++ Libraries; 2022.  http://www.boost.org/.
  * Levenshtein et al. [1966] Levenshtein VI, et al.  Binary codes capable of correcting deletions, insertions, and reversals.  In: Soviet physics doklady, vol. 10 Soviet Union; 1966. p. 707–710.  https://nymity.ch/sybilhunting/pdf/Levenshtein1966a.pdf.

##  Appendix A CATSP Solver with Visual Attractiveness Metrics

1

Input: $V,T,L,f_{o}$

2

3$cost_{best}\leftarrow\infty$

4$V$ $\leftarrow$ Assign missing zone IDs to stops ($V$)

5$M=\emptyset\leftarrow$Matrix of average cluster-to-cluster (zone-to-zone)
travel times

6for  _$g=1,2,\dots,m$_ do

7

8 for  _$h=1,2,\dots,m\quad h\neq g$_ do

9

10 for  _$i\in V_{g}$_ do

11

12 for  _$j\in V_{h}$_ do

13

14 $M_{gh}=M_{gh}+T_{ij}$

15

16 $M_{gh}=\frac{M_{gh}}{|V_{g}|+|V_{h}|}$

17

18

19for  _$iteration\leq L$_ do

20 Zones solution $z\leftarrow$ Random_order($Z$)

21 $z\leftarrow$ relocate_zones($z,M$)

22 Solution $x\leftarrow$ Randomize ($z,N$)

23 $x\leftarrow$ VND($x,neighborhoods,T$)

24 $\upsilon_{k}\leftarrow$ compute_visual_attractiveness_metric($x$)

25 if  _$f_{o}(x,\upsilon_{k}) <cost_{best}$_ then

26 $cost_{best}\leftarrow f_{o}(x,\upsilon_{k})$

27 $x_{best}\leftarrow x$

28

29

30return _$x_{best}$_

Algorithm 5 CATSP solver with visual attractiveness metrics. Adapted from [12]

##  Appendix B CATSP Solver with Historical Routes Information

1

Input: $V,T,Q,L,f_{o}$

2$V$ $\leftarrow$ Assign missing zone IDs to stops ($V$)

Extract historical information phase

3$F$ $\leftarrow$ Create vector field ($Q$)

4$H_{ij}$, $A_{ij}$ $\leftarrow$ Calculate transition probability matrices
($F,V$)

5$M_{ij}$ $\leftarrow$ Calculate zone-to-zone travel-time ($V$)

Route prediction phase: Two-Stage GRASP

6$cost_{best}\leftarrow\infty$

7for  _$i\leq L$_ do

8 Zones solution $z\leftarrow$ Random_order($Z$)

9 $z\leftarrow$ relocate_zones($z,M,H,A$)

10 Solution $x\leftarrow$ Random_order ($z,N$)

11 $x\leftarrow$ VND($x,neighborhoods$)

12 if  _$f_{o}(x,M,H,A) <cost_{best}$_ then

13 $cost_{best}\leftarrow f_{o}(x,M,H,A)$

14 $x_{best}\leftarrow x$

15

16

17return _$x_{best}$_

Algorithm 6 CATSP solver with historical routing information. Adapted from
[12]

##  Appendix C Bi-objective CATSP Solver with Historical Routes Information

1

Input: $Q,V,L,n_{max}$, $f_{1}^{max}$, $f_{2}^{min}$, $\delta^{min}$

Output: Archive of non-dominated solutions $A_{sol}$

2

Extract historical information phase

3$F$ $\leftarrow$ Create vector field ($Q$)

4$H_{ij}$, $A_{ij}$ $\leftarrow$ Calculate transition probability matrices
($V$)

Routes prediction phase: Heuristic Box Splitting Two-Stage GRASP

5$V$ $\leftarrow$ Assign missing zone IDs ($V$)

6$M_{ij}$ $\leftarrow$ Calculate zone-to-zone travel-time ($V$)

7$x\leftarrow ROH(f_{1},\infty,H,A,M)$

8$A_{sol}\leftarrow\\{x\\}$ , $z^{1}\leftarrow(x_{1},x_{2})$ ,
$z^{2}\leftarrow(f_{1}^{max},f_{2}^{min})$,
$R_{areas}\leftarrow\\{$R$(z^{1},z^{2})\\}$

9for  _$i\leq n_{max}$_ do

10 Find R$(y^{1},y^{2})\in R_{areas}$ with maximal area

11 $c\leftarrow\frac{1}{2}(y_{2}^{1}+y_{2}^{2})$

12 $x\leftarrow ROH(f_{1},c,H,A,M)$

13 if  _$x=NULL$ or $x$ is dominated_ then

14

15 $y^{2}\leftarrow(y_{1}^{2},c)$

16 else

17 Update archive $A_{sol}$ with $x$

18 if  _there exists{R’$(z^{1},z^{2})\in R_{areas}$ $|z_{1}^{1}<x_{1}\leq
z_{1}^{2}$}_ then

19

20 $z*\leftarrow(x_{1},max\\{z_{2}^{2},c\\})$

21 if  _ $z_{2}^{1}-z_{2}^{*}>\delta_{min}$_ then

22

23 Add R$(z^{1},z^{*})$ to $R_{areas}$

24

25

26 if  _there exists{R”$(z^{1},z^{2})\in R_{areas}$ $|z_{2}^{1}\geq x_{2}\geq
z_{2}^{2}$}_ then

27

28 $z^{*}\leftarrow(max\\{x_{1},z_{1}^{1}\\},x_{2}$)

29 if  _$z_{2}^{*}-z_{2}^{2} >\delta_{min}$_ then

30

31 Add $R(z*,z2)$ to $R_{areas}$

32

33

34 Delete R’ and R” if they exist

35 forall  _R $(z^{1},z^{2})\in R_{areas}$ with $z^{1}$ dominated by $x$_ do

36

37 Delete R from $R_{areas}$

38

39 if  _$R_{areas}=\emptyset$_ then

40

41 Stop

42

43 return _$A_{sol}$_

Algorithm 7 Bi-objective Two-Stage GRASP with HBS and DM. Based on [31, 12].

##  Appendix D Scoring metric

The scoring metric used on the 2021 Amazon Last-Mile Routing Challenge is
shown in Equation 19.

$score(x,B)=\frac{SD(x,B)\cdot ERP_{norm}(x,B)}{ERP_{e}(x,B)}$ (19)

$B$ is the actual executed route sequence, and $x$ is the predicted route
sequence. $D(x,B)$ measures the sequence deviation of $x$ with respect to $B$.
$E_{norm}(x,B)$ denotes the edit distance weighted with a real penalty applied
to sequences $x$ and $B$ using normalized travel times. $E_{op}(x,B)$ counts
the number of edits required to transform sequence $x$ into sequence $B$. If
both sequences are identical, $x=B$, then $score(x,B)$ is replaced with a
score of 0.

The sequence deviation $D(x,B)$ is calculated as follows. Assuming both routes
visit the same $n+1$ stops ($n$ stops and one depot station). First, we take
each stop from $x$ and trace its position $g$ in $B$. Next, we create vector
$(g_{0},g_{1},...,g_{n})$. Finally, we compute the sequence deviation using
Equation 20.

$D(x,B)=\frac{2}{n(n-1)}\sum_{i=1}^{n}(|g_{i}-g_{i-1}|-1)$ (20)

The number of edit operations $E_{op}(x,B)$ is calculated using the
Levenshtein distance introduced by Levenshtein et al. [33]. Let
$E=E_{0},E_{1},...,E_{n}$ be the minimum number of edit operations to
transform sequence $x$ into sequence $B$. Additionally, let
$C=c_{0},c_{1},...,c_{n}$ be the normalized travel time costs associated with
the operations performed. Each operation of insertion, deletion, and
substitution has a unit cost. Then, the normalized edit distance similarity
metric is calculated as shown in Equation 21.

$E_{norm}(x,B)=\sum_{i=0}^{n}c_{i}E_{i}$ (21)

The more a predicted route sequence $x$ coincides with the actual route $B$,
the lower the score. A score of zero means they are exactly the same sequence
333Details of the scoring algorithm and related implementations can be found
at https://github.com/MIT-CAVE/rc-cli/tree/main/scoring.