Predicting Drivers' Route Trajectories in Last-Mile Delivery Using A
  Pair-wise Attention-based Pointer Neural Network
In last-mile delivery, drivers frequently deviate from planned delivery
routes because of their tacit knowledge of the road and curbside
infrastructure, customer availability, and other characteristics of the
respective service areas. Hence, the actual stop sequences chosen by an
experienced human driver may be potentially preferable to the theoretical
shortest-distance routing under real-life operational conditions. Thus, being
able to predict the actual stop sequence that a human driver would follow can
help to improve route planning in last-mile delivery. This paper proposes a
pair-wise attention-based pointer neural network for this prediction task using
drivers' historical delivery trajectory data. In addition to the commonly used
encoder-decoder architecture for sequence-to-sequence prediction, we propose a
new attention mechanism based on an alternative specific neural network to
capture the local pair-wise information for each pair of stops. To further
capture the global efficiency of the route, we propose a new iterative sequence
generation algorithm that is used after model training to identify the first
stop of a route that yields the lowest operational cost. Results from an
extensive case study on real operational data from Amazon's last-mile delivery
operations in the US show that our proposed method can significantly outperform
traditional optimization-based approaches and other machine learning methods
(such as the Long Short-Term Memory encoder-decoder and the original pointer
network) in finding stop sequences that are closer to high-quality routes
executed by experienced drivers in the field. Compared to benchmark models, the
proposed model can increase the average prediction accuracy of the first four
stops from around 0.2 to 0.312, and reduce the disparity between the predicted
route and the actual route by around 15%.
Baichuan Mo, Qing Yi Wang, Xiaotong Guo, Matthias Winkenbach, Jinhua Zhao
Predicting Drivers' Route Trajectories in Last-Mile Delivery Using A Pair-wise Attention-based Pointer Neural Network

###### Abstract

In last-mile delivery, drivers frequently deviate from planned delivery routes because of their tacit knowledge of the road and curbside infrastructure, customer availability, and other characteristics of the respective service areas. Hence, the actual stop sequences chosen by an experienced human driver may be potentially preferable to the theoretical shortest-distance routing under real-life operational conditions. Thus, being able to predict the actual stop sequence that a human driver would follow can help to improve route planning in last-mile delivery. This paper proposes a pair-wise attention-based pointer neural network for this prediction task using drivers' historical delivery trajectory data. In addition to the commonly used encoder-decoder architecture for sequence-to-sequence prediction, we propose a new attention mechanism based on an alternative specific neural network to capture the local pair-wise information for each pair of stops. To further capture the global efficiency of the route, we propose a new iterative sequence generation algorithm that is used after model training to identify the first stop of a route that yields the lowest operational cost. Results from an extensive case study on real operational data from Amazon's last-mile delivery operations in the US show that our proposed method can significantly outperform traditional optimization-based approaches and other machine learning methods (such as the Long Short-Term Memory encoder-decoder and the original pointer network) in finding stop sequences that are closer to high-quality routes executed by experienced drivers in the field. Compared to benchmark models, the proposed model can increase the average prediction accuracy of the first four stops from around 0.2 to 0.312, and reduce the disparity between the predicted route and the actual route by around 15%.

keywords: Route planning, Trajectory prediction, Sequence-to-sequence model, Last-mile delivery, Pointer network, Attention +
Footnote †: journal:

## 1 Introduction

The optimal planning and efficient execution of last-mile delivery routes is becoming increasingly important for the business operations of many logistics service providers around the globe for a variety of reasons. E-commerce volumes are growing rapidly and make up a constantly growing share of overall retail sales. For instance, in the US, the share of e-commerce sales in total retail sales has grown from around 4% in 2010 to around 13% in 2021. Even by the end of 2019, i.e., before the outbreak of the COVID-19 pandemic,
it had reached 11% (US Census Bureau, 2021). Undoubtedly, the pandemic further accelerated the growth of e-commerce (postnord, 2021; McKinsey and Company, 2021). In the medium to long run, its growth will continue to be fueled by an ongoing trend towards further urbanization, which is particularly pronounced in developing and emerging economies (United Nations Department of Economic and Social Affairs, 2019). The share of the global population living in urban areas is currently projected to rise from around 55% in 2018 to around 68% by 2050. The associated increase in population density in most urban areas will likely lead to growing operational uncertainties for logistics service providers, as increasing congestion levels, less predictable travel times, and scarce curb space make efficient and reliable transport of goods into and out of urban markets increasingly challenging (Rose et al., 2016).

As a result of the continued boom of e-commerce and constantly growing cities, global parcel delivery volumes have been increasing rapidly in recent years and are expected to continue to do so. Across the 13 largest global markets, including the US, Brazil, and China, the volume of parcels delivered more than tripled from 43 billion in 2014 to 131 billion in 2020 (Pitney Bowes, 2020). At the same time, customer expectations towards last-mile logistics services are rising. For instance, there is a growing demand for shorter delivery lead times, including instant delivery services and same-day delivery, as well as customer-defined delivery preferences when it comes to the time and place of delivery (Lim and Winkelbach, 2019; Cortes and Suzuki, 2021; Snoeck and Winkelbach, 2021). The rapid growth and increasing operational complexity of urban parcel delivery operations also amplifies their negative externalities, including their contribution to greenhouse gas and other pollutant emissions, public health safety risks, as well as overall urban congestion and a corresponding decline in overall mobility and accessibility of cities (Jaller et al., 2013; World Economic Forum, 2020).

When applied to realistically sized instances of a last-mile delivery problem, solving the underlying traveling salesman problem (TSP) or vehicle routing problem (VRP) to (near) optimality becomes challenging, as both problem classes are known to be NP-hard. Traditional TSP and VRP formulations aim to minimize the total distance or duration of the route(s) required to serve a given set of delivery stops. The operations research literature has covered the TSP, VRP, and their many variants extensively, and in recent years important advances have been made with regards to solution quality and computational cost. However, in practice, many drivers, with their own tacit knowledge of delivery routes and service areas, divert from seemingly optimal routes for reasons that are difficult to encode in an optimization model directly. For example, experienced drivers may have a better understanding of which roads are hard to navigate, at which times traffic is likely to be bad, when and where they can easily find parking, and which stops can be conveniently served together. Therefore, compared to the theoretically optimal (i.e., distance or time minimizing) route, the _deviated actual route sequence_ chosen by an experienced human driver is potentially preferable under real-life operational conditions.

An important challenge in today's last-mile delivery route planning is therefore to leverage historical route execution data to propose planned route sequences that are close to the actual trajectories that would be executed by drivers, given the delivery requests and their characteristics. Note that, while distance and time-based route efficiency is still an important factor for planning route sequences, it is not the sole objective, as tacit driver knowledge is also incorporated in the proposed route sequences. Unlike a typical VRP in which the number of vehicles and their respective route sequences need to be determined simultaneously, in this study, we focus on solving a problem that is similar to a TSP at the individual vehicle level. That is, we aim to solve a stop sequence to serve a given set of delivery requests, and expect that the proposed stop sequence is as close to the actual trajectories that would be executed by drivers as possible.
To this end, we propose a pair-wise attention-based pointer neural network to predict the actual route sequence taken by delivery drivers using drivers' historical delivery trajectory data. The proposed model follows a typical encoder-decoder architecture for the sequence-to-sequence prediction. However, unlike previous studies, we propose a new attention mechanism based on an alternative specific neural network (ASNN) to capture the local pair-wise information for each stop pair. To further capture the global efficiency of the route (i.e., its operational cost in terms of total distance or duration), after model training, we propose a new sequence generation algorithm that iterates over different first stops and selects the route with the lowest operational cost.

The main contribution of this paper is three-fold: First, we propose a new ASNN-based attention mechanism to capture the local information between pairs of stops (e.g., travel time, geographical relation), which can be well adapted to the original pointer network framework for sequence prediction. Second, we propose a new sequence generation algorithm that iterates over different first stops in the predicted route sequences and selects the lowest operational cost route. The intuition is that the stop-to-stop relationship (referred to as the _local view_) is easier to learn from data than the stop sequence of the route as a whole (referred to as the _global view_). Lastly, we apply our proposed method to a large set of routes executed by Amazon delivery drivers in the US. The results show that our proposed model can outperform traditional optimization-based approaches and other machine learning methods in finding stop sequences that are closer to high-quality routes executed by experienced drivers in the field.

The remainder of this paper is structured as follows. In Section 2 we define the problem setting under investigation in a more formal way. Section 3 then reviews previous studies in the literature related to this paper. Section 4 presents our methodology and elaborates on the detailed architecture of the proposed pair-wise attention-based pointer neural network. Section 5 presents the experimental setup and numerical results of our case study, applying our proposed method to real-world data made available by the Amazon Last-Mile Routing Research Challenge (Merchan et al., 2022; Winkenbach et al., 2021). Section 6 concludes this paper and discusses future research directions.

## 2 Problem Setting

In the last-mile delivery routing problem considered here, a set of stops \(\mathcal{S}=\{s_{1},...,s_{n}\}\) to be served by a given delivery vehicle is given to the route planner. The planner's objective is to find the optimal stop sequence that has the minimal operational cost. In this case, we consider total cost as total travel time. The planner is given the expected operational cost (i.e., travel times) between all pairs of stops \((s_{i},s_{j})\). The theoretically optimal stop sequence, denoted by \((s_{(1)}^{\mathsf{T}},...,s_{(n)}^{\mathsf{T}})\), can be found by solving a TSP formulation. This stop sequence is referred to as the planned stop sequence. However, as discussed in Section 1, minimizing the theoretical operational cost (i.e., total travel time) of the route may not capture drivers' tacit knowledge about the road network, infrastructure, and recipients. Therefore, the actual driver executed stop sequence \((s_{(1)},...,s_{(n)})\) can be different from the planned route sequence. Note that here, \(s_{(i)}\in\mathcal{S}\) denotes the \(i\)-th stop that is actually visited by the driver.

The objective of the model presented in this study is to predict the actual driver executed sequence \((s_{(1)},...,s_{(n)})\) given a set of stops \(\mathcal{S}\) and the corresponding delivery requests and characteristics \(X^{\mathcal{S}}\) (such as the number of packages, estimated service time for each package, geographical information for each stop, travel time between each stop pairs, etc.). All drivers are assumed to start their routes from a known depot \(D^{\mathcal{S}}\) and return back to \(D^{\mathcal{S}}\). Therefore, the complete trajectory should be a tour \((D^{\mathcal{S}},s_{(1)},...,s_{(n)},D^{\mathcal{S}})\). For the convenience of model description, we ignore the depot station in the sequence.
Figure 1 provides a simple example for illustration. In this example, we are given four stops \(\mathcal{S}=\{s_{1},s_{2},s_{3},s_{4}\}\) and a depot \(D^{\mathcal{S}}\). The planned stop sequence for the driver is \((s_{4},s_{1},s_{2},s_{3})\), while the actual stop sequence executed by the driver is \((s_{4},s_{2},s_{1},s_{3})\). The proposed model aims to predict the actual sequence \((s_{4},s_{2},s_{1},s_{3})\) given the depot location \(D^{\mathcal{S}}\), the set of stops to be visited \(\mathcal{S}\), and characteristics of the stops \(X^{\mathcal{S}}\). This problem setup is inspired by the Amazon Last-Mile Routing Research Challenge (cf., Winkelbach et al., 2021). Note that this study only focuses on the stop sequence prediction. The routing between stops is not considered. It is assumed that the drivers always take the optimal route between stops, which is reflected by the travel time matrix between stops in our problem setup.

## 3 Literature Review

The problem setting defined in Section 2 involves both solving a cost-minimizing routing problem (i.e., the TSP) and capturing tacit driver knowledge to learn systematic deviation of drivers from the planned and theoretically optimal stop sequences. Therefore, we will first review the extant literature on the TSPs and its most relevant variants. We will then go through various machine learning approaches that have been proposed by the extant literature to generate sequences, with a section on methods specifically for solving the TSP. Note that although these machine learning approaches are used to solve the TSP instead of the actual routes taken by drivers, their architectures may be helpful to learn the actual route as well.

### Travelling salesman problems

First, given the travel times between stops, a solution to the TSP, which finds the route with the minimum cost or distance (i.e., the planned route), can be a close approximation of the actual route. Since the drivers are paid for the number of packages delivered, all drivers' goal is to deliver the packages in the minimum amount of time. Most of the drivers do follow large parts of the planned routes.

The TSP is a well-known NP-hard problem that has been studied extensively over the last century, with a lot of books and review papers published on its history, formulations, solution approaches, and applications (Applegate et al., 2006; Matai et al., 2010; Davendra and Bialic-Davendra, 2020). An overview of the relevant TSP variants and solution approaches are presented below.

The basic setup of TSP has one traveler and requires the traveler to return to the starting point after visiting each node exactly once, and that the traveling cost matrix (represented by distance and/or time) is

Figure 1: Illustrative example of the problem setting.
symmetric (cost between \(i\) and \(j\) is the same with that between \(j\) and \(i\)). In most real-world applications, the basic setup needs to be modified. For example, the cost matrix, if represented by travel times, is likely asymmetric. This variant of TSP is thus named asymmetric TSP (ATSP) (Jonker and Volgenant, 1983). In some applications, the vehicle does not need to return to the original depot (Traub et al., 2021), or it can charge/refuel and potentially load additional delivery items at intermediate stops (Kucukoglu et al., 2019). In many last-mile delivery applications, some packages are time-sensitive, and therefore time window constraints to their delivery need to be considered in a so-called TSP with time windows (TSPTW) (da Silva and Urrutia, 2010; Mladenovic et al., 2012). In large systems, there might be more than one salesman serving a set of stops, resulting in multiple traveling salesmen problems (MTSPs) (Cheikhrouhou and Khoufi, 2021).

Different variants of TSP further impose different constraints on the solution. While some problems can be reduced to the basic setup in the formulation stage, others require more versatile solution algorithms. In general, the solution approaches to the TSP can be divided into exact approaches and approximate approaches. Exact approaches include branch-and-cut (Yuan et al., 2020) and branch-and-bound (Salman et al., 2020). Since the TSP is a well-known NP-hard problem, exact approaches can only be applied on problems of smaller scale, or aid in heuristics to cut the solution space. Among approximate approaches, there are heuristics designed for the TSP specifically, as well as meta-heuristics that are generic and treat the problem like a blackbox. The most commonly used heuristics and meta-heuristics include nearest neighbor searches, local searches, simulated annealing, and genetic algorithms. A more comprehensive review of existing solution approaches can be found in Halim and Ismail (2017); Purkayastha et al. (2020). Despite the TSP being NP-hard, modern mixed-integer optimization solvers (e.g., Gurobi, CPLEX, or GLPK) can solve it efficiently for real-world instances by combining exact approaches with heuristics.

### Sequence-to-sequence prediction using deep learning

The TSP and its variants are a viable option for sequence generation only when the objective is clearly defined. They fall short when the sequence generation problem does not have a well-defined cost-minimization objective. In a lot of applications, the rule of sequence generation cannot be simply defined and optimized.

A standard example for a sequence learning problem is machine translation, where a sequence of words in one language needs to be translated to another language. Another type of sequence learning is time series modeling, where a sequence of historical observations is given to predict future states of the system. In both cases, the primary modeling task is to learn the sequence generation rules. In recent years, deep learning has successfully achieved great performance in various settings of sequence generation. These models are often referred to as sequence-to-sequence (seq2seq) models.

seq2seq models often consist of an encoder and a decoder, where the encoder encodes the input sequence into a fixed-length vector representation, and the decoder generates a sequence based on the generated vector representation. Most encoder-decoder architectures adopt recurrent neural network (RNN) layers and its variants such as Long Short-Term Memory (LSTM) (Hochreiter and Schmidhuber, 1997) and gated recurrent layers (GRU) (Cho et al., 2014) to learn long-range dependencies. Early works using LSTM alone were able to generate plausible texts (Graves, 2013) and translate between English and French (Sutskever et al., 2014) with long-range dependencies. Chung et al. (2014) demonstrate the superiority of GRU compared to LSTMs in music and speech signal modeling.

Attention-based mechanisms, first introduced by Bahdanau et al. (2015), have been shown to be a great addition since it allows the decoder to selectively attend to parts of the input sequence and relieves the encoder of the task of encoding all the information into a fixed-length vector representation. Most sequence generation
problems benefit from keeping track of long-range dependencies and global context while decoding. To address that, multi-level attention was proposed to capture the local and global dependency, and has shown to be effective in speech recognition (Chorowski et al., 2015), text generation (Liu et al., 2018), and machine translation tasks (Luong et al., 2015).

The encoder-decoder architecture combined with attention is very versatile, and it can be combined with other deep learning architectures to perform sequence learning in addition to language tasks. The LSTM and attention architecture is applied to semantic trajectory prediction (Karatzoglou et al., 2018), text summarization (Liang et al., 2020), demand modelling (Ren et al., 2020), and wind power forecasting (Zhang et al., 2020). When the goal is set to recover the original sequence, unsupervised learning of molecule embedding can be obtained for downstream classification tasks (Xu et al., 2017). When the spatial dimension is added, a convolutional neural network (CNN) layer can be added, and the dimension of the sequence generated can be expanded. For example, Wang et al. (2020) predict a city's crowd flow patterns, and Wu et al. (2020) generate 3D shapes via sequentially assembling different parts of the shapes.

While RNN-based architectures are still a widely adopted choice for seq2seq modeling, attention can also be used as a standalone mechanism for seq2seq translations independent of RNNs. The idea was proposed by Vaswani et al. (2017) in an architecture named transformer. Without recurrence, the network allows for significantly more parallelization, and is shown to achieve superior performance in experiments, and powered the popularity of transformer-based architectures in various sequence generation tasks (Huang et al., 2018; Lu et al., 2021). A separate line of work by Zhang et al. (2019) also demonstrated that a hierarchical CNN model with attention outperforms the traditional RNN-based models.

### Using deep learning to generate TSP solutions

The above seq2seq translation mechanisms work well when the input data is naturally organized as a sequence, and the output sequence corresponds to the input sequence, such as in music and language. However, in our paper, the input is an unordered sequence, and the output has the same but re-ordered elements of the same input sequence. In this case, the concept of attention is helpful and has been successfully used to produce solutions to the TSP. The pointer network, proposed by Vinyals et al. (2015) and further developed in Vinyals et al. (2016), uses attention to select a member of the input sequence at each decoder step. While it is not required that the input sequence is ordered, an informative ordering could improve the performance (Vinyals et al., 2016).

While the original pointer network was solved as a classification problem and cross-entropy loss was used, it is not necessarily the most efficient choice. The cross-entropy loss only distinguishes between a correct prediction and an incorrect prediction. But in instances like routing, the distances between the predicted position and the correct position, as well as the ordering of subsequences, could incur different costs in practice. Further developments in solving TSP with machine learning methods involve reinforcement learning (RL), which enables the optimization of custom evaluation metrics (Bello et al., 2019; Kool et al., 2019; Ma et al., 2019; Liu et al., 2020). Joshi et al. (2019) compared the performance of RL and supervised learning (SL) on TSP solutions and found that SL and RL models achieve similar performance when the graphs are of similar sizes in training and testing, whereas RL models have better generalizability over variable graph sizes. However, RL models require significantly more data points and computational resources, which is not always feasible.

Although this seq2seq and attention framework has only been used to reproduce TSP solutions, it provides an opportunity to learn and incorporate additional information beyond the given travel times and potentially
learn individual differences when more information is given to the neural network. In this paper, we combine the ideas of seq2seq modeling and attention to predict the actual route executed by a driver.

## 4 Methodology

This section details the methodology proposed to address the problem. First, the high-level seq2seq modelling framework is introduced, followed by the explanation of the novel pair-wise attention and sequence generation and selection mechanism used within the modelling framework.

### Sequence-to-sequence modeling framework

Let the input sequence be an arbitrarily-ordered sequence \((s_{1},...,s_{n})\). Denote the output sequence as \((\hat{s}_{(1)},...,\hat{s}_{(n)})\). Let \(c_{i}\) indicate the "position index" of stop \(\hat{s}_{(i)}\) with respect to the input sequence (where \(c_{i}\in\{1,...,n\}\)). For example, for input sequence \((B,A,C)\) and output sequence \((A,B,C)\), we have \(c_{1}=2\), \(c_{2}=1\), \(c_{3}=3\), which means the first output stop \(A\) is in the second position of the input sequence \((B,A,C)\) and so on.

The seq2seq model computes the conditional probability \(\mathbb{P}(c_{1},...,c_{n}\mid\mathcal{S};\theta)\) using a parametric neural network (e.g., recurrent neural network) with parameter \(\theta\), i.e.,

\[\mathbb{P}(c_{1},...,c_{n}\mid\mathcal{S},X^{\mathcal{S}};\theta)=\mathbb{P} (c_{1}\mid\mathcal{S},X^{\mathcal{S}};\theta)\cdot\prod_{i=2}^{n}\mathbb{P}(c _{i}\mid c_{1},...,c_{i-1},\mathcal{S},X^{\mathcal{S}};\theta) \tag{1}\]

The parameters of the model are learnt by empirical risk minimization (maximizing the conditional probabilities on the training set), i.e.,

\[\theta^{*}=\operatorname*{arg\,max}_{\theta}\sum_{\mathcal{S}}\mathbb{P}(c_{1 },...,c_{n}\mid\mathcal{S},X^{\mathcal{S}};\theta) \tag{2}\]

where the summation of \(\mathcal{S}\) is over all training routes. In the following section, we will elaborate how \(\mathbb{P}(c_{i}\mid c_{1},...,c_{i-1},\mathcal{S},X^{\mathcal{S}};\theta)\) is calculated using the pair-wise attention-based pointer neural network.

### Pair-wise attention-based pointer neural network

Figure 2 uses a four-stop example to illustrate the architecture of the proposed model. The whole model is based on the LSTM encoder and decoder structure. In particular, we use one LSTM (i.e., encoder) to read the input sequence, one time step at a time, to obtain a large fixed dimensional vector representation, and then to use another LSTM (i.e., decoder) to extract the output sequence. However, different from the typical seq2seq model, we borrow the idea of the pointer network (Vinyals et al., 2015) to add a pair-wise attention mechanism to predict the output sequence based on the attention mask over the input sequence. The pair-wise attention is calculated based on an ASNN which was previously used for travel mode prediction (Wang et al., 2020). Model details will be shown in the following sections.

Intuitively, the LSTM encoder and decoder aim to capture the global view of the input information (i.e., overall sequence pattern) by embedding the input sequence to hidden vector representation. While the ASNN-based pair-wise attention aims to capture the local view (i.e., the relationship between two stops). Our experiments in Section 5 demonstrate the importance of both global and local views in the sequence prediction.
#### 4.2.1 LSTM encoder.

Given an arbitrary stop sequence \((s_{1},...,s_{n})\) as the input, let \(x_{i}\in\mathbb{R}^{K}\) be the features of stop \(s_{i}\), where \(x_{i}\) may include the package information, the customer information, and the geographical information of the stop \(s_{i}\). \(K\) is the number of features. The encoder computes a sequence of encoder output vectors \((e_{1},...,e_{n})\) by iterating the following:

\[h_{i}^{\text{E}},e_{i}=\text{LSTM}(x_{i},h_{i-1}^{\text{E}};\ \theta^{\text{E}}) \quad\forall i=1,...,n \tag{3}\]

where \(h_{i}^{\text{E}}\in\mathbb{R}^{K_{h}^{\text{E}}}\) is the encoder hidden vector with \(h_{0}^{\text{E}}:=0\). \(e_{i}\in\mathbb{R}^{K_{e}}\) is the encoder output vector. \(K_{h}^{\text{E}}\) and \(K_{e}\) are corresponding vector dimensions. \(\theta^{\text{E}}\) is the learnable parameters in an encoder LSTM cell. The calculation details of an LSTM cell can be found in A. The encoding process transforms a sequence of features \((x_{1},...,x_{n})\) into a sequence of embedded representation \((e_{1},...,e_{n})\). And the hidden vector of the last time step (\(h_{n}^{\text{E}}\)) includes the global information of the whole sequence, which will be used for the LSTM decoder.

#### 4.2.2 LSTM decoder.

The role of a decoder in the traditional seq2seq model (Figure 4) is to predict a new sequence one time step at a time. However, in the pointer network structure with attention, the role of the decoder becomes producing

Figure 3: Illustration of LSTM encoder

Figure 2: Overall architecture of the pair-wise attention-based pointer neural network (adapted from Vinyals et al. (2015))
a vector to modulate the pair-wise attention over inputs. Denote the output sequence as \((\hat{s}_{(1)},...,\hat{s}_{(n)})\). Let \(x_{(i)}\) be the feature of stop \(\hat{s}_{(i)}\).

At decoder step \(i\), we have

\[h^{\text{D}}_{(i+1)},d_{(i)}=\text{LSTM}\left(\begin{bmatrix}x_{(i)}\\ w_{(i)}\end{bmatrix},h^{\text{D}}_{(i)};\ \theta^{\text{D}}\right)\quad\forall i=0,1,...,n \tag{4}\]

where \(h^{\text{D}}_{(i)}\in\mathbb{R}^{K^{\text{D}}_{h}}\) is the decoder hidden vector with \(h^{\text{D}}_{(0)}=h^{\text{E}}_{n}\), \(d_{(i)}\in\mathbb{R}^{K_{d}}\) is the decoder output vector, \(K^{\text{D}}_{h}\) and \(K_{d}\) are corresponding vector dimensions, and \(\theta^{\text{D}}\) are learnable parameters of the decoder LSTM cell. Note that we set \(x_{(0)}=x_{D}\) and \(d_{(0)}=d_{D}\), representing the features and the decoder output of the depot, respectively. \(w_{(i)}\) is the context vector calculated from the attention component, which will be explained in the next section.

#### 4.2.3 ASNN-based pair-wise attention.

The pair-wise attention aims to aggregate the global and local information to predict the next stop. Specifically, at each decoder time step \(i\in\{0,...,n\}\), we know that the last predicted stop is \(\hat{s}_{(i)}\). To predict \(\hat{s}_{(i+1)}\), we consider all candidate stops \(s_{j}\in\mathcal{S}\), which is the set of all stops not yet visited. We want to evaluate how possible that \(s_{j}\) will be the next stop of \(\hat{s}_{(i)}\). The information of the stop pair \(\hat{s}_{(i)}\) and \(s_{j}\) can be represented by the following concatenated vector:

\[v^{j}_{(i)}=\text{concat}(z^{j}_{(i)},\ \phi(x_{(i)},x_{j}),\ d_{(i)},\ e_{j}) \tag{5}\]

where \(z^{j}_{(i)}\) is a vector of features associated with the stop pair (such as travel time from \(\hat{s}_{(i)}\) to \(s_{j}\)), and \(\phi(x_{(i)},x_{j})\) represents a feature processing function to extract the pair-wise information from \(x_{(i)}\) and \(x_{j}\). For example, \(\phi(\cdot)\) may return geographical relationship between stops \(\hat{s}_{(i)}\) and \(s_{j}\), and it may also drop features not useful for the attention calculation. Intuitively, \(z^{j}_{(i)}\) and \(\phi(x_{(i)},x_{j})\) contains only local information of the stop pair, while \(d_{(i)}\) and \(e_{j}\) contain the global information of the whole stop set and previously visited stops.

Figure 4: Illustration of LSTM decoder
Given the pair-wise information vector \(v^{j}_{(i)}\), we can calculate the attention of stop \(\hat{s}_{(i)}\) to stop \(s_{j}\) as:

\[u^{j}_{(i)} =\text{ASNN}(v^{j}_{(i)};\ \theta^{\text{A}})\quad\forall i,j=1,...,n \tag{6}\] \[a^{j}_{(i)} =\frac{\exp(u^{j}_{(i)})}{\sum_{j^{\prime}=1}^{n}\exp(u^{j^{\prime }}_{(i)})}\quad\forall i,j=1,...,n \tag{7}\]

where \(a^{j}_{(i)}\in\mathbb{R}\) is attention of stop \(\hat{s}_{(i)}\) to stop \(s_{j}\). \(\text{ASNN}(\cdot;\ \theta^{\text{A}}))\) is a multilayer perception (MLP) with the output dimension of one (i.e., \(u^{j}_{(i)}\in\mathbb{R}\)). \(\theta^{\text{A}}\) are the learnable parameters of the ASNN. The name "alternative specific" is because the same parametric network will be applied on all alternative stops \(s_{j}\in\mathcal{S}\) separately (Wang et al., 2020). Finally, we calculate the conditional probability to make the prediction:

\[\mathbb{P}(c_{i+1}=j\ |\ c_{1},...,c_{i},\mathcal{S},X^{ \mathcal{S}};\theta)=a^{j}_{(i)}\quad\forall i=0,1,...,n,\ j=1,...,n \tag{8}\] \[\hat{s}_{(i+1)}=\operatorname*{arg\,max}_{s_{j}\in\mathcal{S} \setminus\mathcal{S}^{\text{V}}_{(i)}}a^{j}_{(i)}\quad\forall i=0,1,...,n \tag{9}\]

where \(\mathcal{S}^{\text{V}}_{(i)}=\{\hat{s}_{(1)},...,\hat{s}_{(i)}\}\) is the set of stops that have been predicted (i.e., previously visited) until decoder step \(i\). Eqs. 8 and 9 indicate that the predicted next stop at step \(i\) is the one with highest attention among all stops that have not been visited.

The pair-wise attention framework also leverages the attention information as the input for the next step. This was achieved by introducing the context vector (Bahdanau et al., 2015):

\[w_{(i)}=\sum_{j=1}^{n}a^{j}_{(i)}\cdot e_{j} \tag{10}\]

The context vector is a weighted sum of all the encoder output vectors with attention as the weights. As the attention provides the emphasis for stop prediction, \(w_{(i)}\) helps to incorporate the encoded representation of the last predicted stop for the next stop prediction. The inputs for the next LSTM cell thus will be the

Figure 5: Illustration of ASNN-based pair-wise attention
concatenation of the stop features and \(w_{(i)}\), i.e., \(\begin{bmatrix}x_{(i)}\\ w_{(i)}\end{bmatrix}\).

It is worth noting that, the specific architecture of \(\text{ASNN}(\cdot;\ \theta^{\text{A}}))\) can be flexible depending on the input pair-wise information. For example, if the information includes images or networks, convolutional neural network or graph convolutional networks can be used for better extract features. In this study, we use the MLP for simplification as it already outperforms benchmark models. The key idea is of the ASNN is to share the same trainable parameter \(\theta^{\text{A}}\) for all stop pairs so as to better capture various pair-wise information in the training process.

### Sequence generation and selection

During inference, given a stop set \(\mathcal{S}\), the trained model with learned parameters \(\theta^{*}\) are used to generate the sequence. Typically, in the seq2seq modeling framework, the final output sequence is selected as the one with the highest probability, i.e.,

\[(s_{j_{1}^{*}},...,s_{j_{n}^{*}}),\text{where }j_{1}^{*},...,j_{n}^{*}= \operatorname*{arg\,max}_{j_{1},...,j_{n}\in\mathcal{C}^{\mathcal{S}}}\mathbb{ P}(c_{1}=j_{1},...,c_{n}=j_{n}\mid\mathcal{S},X^{\mathcal{S}};\theta^{*}) \tag{11}\]

where \(\mathcal{C}^{\mathcal{S}}=\{\text{All permutations of }\{1,...,n\}\}\)

Finding this optimal sequence is computationally impractical because of the combinatorial number of possible output sequences. And so it is usually done with the greedy algorithm (i.e., always select the most possible next stop) or the beam search procedure (i.e., find the best possible sequence among a set of generated sequences given a beam size). However, in this study, we observe that the first predicted stop \(\hat{s}_{(1)}\) is critical for the quality of the generated sequence. The reason may be that the local relationship between a stop pair (i.e., given the last stop to predict the next one) is easier to learn than the global relationship (i.e., predict the whole sequence). Hence, in this study, we first generate sequences using the greedy algorithm with different initial stops, and select the one with the lowest operational cost. The intuition behind this process is that, once the first stop is given, the model can follow the learned pair-wise relationship to generate the sequence with relatively high accuracy. For all the generated sequences with different first stops, the one with the lowest operation cost captures the global view of the sequence's quality. Therefore, the final sequence generation and selection algorithm is as follows:

```
0: Trained model, \(\mathcal{S}\)
0: Predicted stop sequence
1:for\(s\) in \(\mathcal{S}\)do
2: Let the first predicted stop be \(\hat{s}_{(1)}=s\)
3: Predict the following stop sequence \((\hat{s}_{(2)},...,\hat{s}_{(n)})\) using the greedy algorithm. Denote the predicted sequence as \(P_{s}\).
4: Calculate the total operation cost of the whole sequence (including depot), denoted as \(OC_{s}\). return\(P_{s^{*}}\) where \(s^{*}=\operatorname*{arg\,min}_{s\in\mathcal{S}}OC_{s}\)
```

**Algorithm 1** Sequence generation

## 5 Case Study

### Dataset

The data used in our case study was made available as part of the Amazon Last Mile Routing Research Challenge (Merchan et al., 2022). The dataset contains a total of 6,112 actual Amazon driver trajectories
for the last-mile delivery from 5 major cities in the US: Austin, Boston, Chicago, Los Angeles, and Seattle. Each route consists of a sequence of stops. Each stop represents the actual parking location of the driver, and the package information (package numbers, package size, and planned service time) associated with each stop is given. The stops are characterized by their latitudes and longitudes, and expected travel time between stops are known.

Figure 6 shows the distribution of the number of stops per route and an example route. Most routes have around 120 to 180 stops, and the maximum observed number of stops is around 250. Figure 5(b) shows an example of an actual driver trajectory in Boston. Since the depot is far from the delivery stops, we attach the complete route (with the depot indicated by a red dot) at the bottom left of the figure, while the main plot only shows the delivery stops.

In this data set, each stop is associated with a zone ID (indicated by different colors in Figure 5(b)). When Amazon generates planned routes for drivers, they usually expect drivers to finish the delivery for one zone first, then go to another zone. And the actual driver trajectories also follow this pattern as shown in Figure 5(b) (but the actual zone sequence may be different from the planned one). Therefore, in this study, we focus on the problem of _zone sequence prediction_. That is, \(s_{i}\) in the case study section now represents a specific zone, \(\mathcal{S}\) represents the set of zones, and \(X^{\mathcal{S}}\) represents zone features. This transformation does not affect the model structure proposed in Section 4. The only difference is that the new problem has a relatively smaller scale compared to the _stop sequence prediction_ because the number of zones in a route is smaller than that of stops. The zone-to-zone travel time is calculated as the average travel time of all stop pairs between the two zones. Figure 7 presents an illustrative example of the relationship between zone and stop sequences. As the dataset does not contain the original planned sequence, we assume the planned zone sequence is the one with the lowest total travel time (generated by a TSP solver, \((s_{1}^{\text{T}},...,s_{n}^{\text{T}})\)). After generating the zone sequence, we can restore the whole stop sequence by assuming that drivers within a specific zone follow an optimal TSP tour. Details of the zone sequence to stop sequence generation can be found in B.

### Experimental setup

We randomly select 4,889 routes for model training and cross-validation, and the remaining 1,223 routes are used to evaluate/test model performance.

We consider a one-layer LSTM for both the encoder and decoder with the hidden unit sizes of 32 (i.e., \(K_{h}^{\text{D}}=K_{e}=K_{h}^{\text{E}}=K_{d}=32\)). And the ASNN is set with 2 hidden layers with 128 hidden units in each layer. We train the model using Adam optimizer with a default learning rate of 0.001 and 30 training epochs. To utilize the planned route information, the input zone sequence for the LSTM encoder is set as the TSP

Figure 7: Relationship between stop sequence and zone sequence.
Figure 6: Description of dataset
result (i.e., lowest travel time). That is, the input sequence \((s_{1},...,s_{n})=(s_{1}^{\intercal},...,s_{n}^{\intercal})\).

In the case study, \(x_{i}\) represents zone features, including the latitude and longitude of the zone center, number of stops in the zone, number of intersections in the zone, number of packages in the zone, total service time in the zone, total package size in the zone, and the travel time from this zone to all other zones. The zone pair features \(z_{(i)}^{j}\) includes the travel time from \(\hat{s}_{(i)}\) to \(s_{j}\) and zone ID relationship characteristics. For example, the zone IDs "B-6.2C" and "B-6.3A" signal that they belong to the higher-level cluster "B-6". As we assume all pair-wise features are captured by \(z_{(i)}^{j}\), \(\phi(x_{(i)},x_{j})\) is not specified in this case study.

Consistent with the Amazon Last Mile Routing Research Challenge, we evaluate the quality of the predicted stop sequences using a "disparity score" defined as follows:

\[R(A,B)=\frac{SD(A,B)\cdot ERP_{\text{norm}}(A,B)}{ERP_{\text{e}}(A,B)} \tag{12}\]

where \(R(A,B)\) is the disparity score for the actual sequence \(A\) and predicted sequence \(B\), and \(SD(A,B)\) is the sequence deviation defined as

\[SD(A,B)=\frac{2}{n(n-1)}\sum_{i=2}^{n}\left(|c_{[B_{i}]}-c_{[B_{i-1}]}|-1\right) \tag{13}\]

where \(n\) is the total number of stops, \(B_{i}\) is the \(i\)-th stop of sequence \(B\), \(c_{[B_{i}]}\) is the index of stop \(B_{i}\) in the actual sequence \(A\) (i.e., its position in sequence \(A\)). In the case of \(A=B\) (i.e., perfectly predicted), we have \(c_{[B_{i}]}-c_{[B_{i-1}]}=1\) for all \(i=2,...,n\), and \(SD(A,B)=0\).

\(ERP_{\text{norm}}(A,B)\) is the Edited Distance with Real Penalty (ERP) defined by the following recursive formula:

\[ERP_{\text{norm}}(A,B)=ERP_{\text{norm}}(A_{2:|A|},B_{2:|B|})+\text{Time}_{ \text{norm}}(A_{1},B_{1}) \tag{14}\]

where \(\text{Time}_{\text{norm}}(s_{i},s_{j})=\frac{\text{Time}(s_{i},s_{j})}{\sum_{ j^{\prime}\in\{1,...,n\}}\text{Time}(s_{i},s_{j^{\prime}})}\) is the normalized travel time from stop \(s_{i}\) to stop \(s_{j}\). \(ERP_{\text{e}}(A,B)\) is the number of edit operations (insertions, substitutions, or deletions) required to transform sequence \(A\) to sequence \(B\) as when executing the recursive \(ERP_{\text{norm}}\) formulation. Hence, the ratio \(\frac{ERP_{\text{norm}}(A,B)}{ERP_{\text{e}}(A,B)}\) represents the average normalized travel time between the two stops involved in each ERP edit operation. In the case of \(A=B\), we have \(\frac{ERP_{\text{norm}}(A,B)}{ERP_{\text{e}}(A,B)}=0\).

The disparity score \(R(A,B)\) describes how well the model-estimated sequence matches the known actual sequence. Lower score indicates better model performance. A score of zero means perfect prediction. The final model performance is evaluated by the mean score over all routes in the test set.

In addition to the disparity score, we also evaluate the prediction accuracy of the first four zones in each route. We choose the first four because the minimum number of zones in a route is four.

### Benchmark models

The following optimization and machine learning models are used as benchmarks to compare with the proposed approach.

**Conventional TSP**. The first benchmark model is the zone sequence generated by conventional TSP, which we treat as the planned route with the lowest travel time.

**ASNN model**. The ASNN component can be trained to predict the next zone given the current zone, and the prediction sequence can be constructed in a greedy way starting from the given depot. The training zone
pairs (including from depot to the first zone) are extracted from all sequences in the training routes. And the input features are the same as the ASNN component in the proposed model except for \((d_{(i)},e_{j})\) (i.e., output vectors from LSTM decoder and encoder, respectively). All hyper-parameters of the ASNN model are the same as the attention component.

Inspired by the importance of the first zone, we also implement another sequence generation method similar to Section 4.3. That is, we go through all zones in a route and assume it is the first zone, then use the trained ASNN to predict the remaining sequence. The final sequence is selected as the one with the lowest travel time.

**LSTM-encoder-decoder**. The LSTM-encoder-decoder (LSTM-E-D) architecture is a typical seq2seq model proposed by Sutskever et al. (2014). The model structure is shown in Figure 8. In the decoder stage, the model outputs the predicted zone based on last predicted zone's information. The model formulation can be written as

\[h_{i}^{\text{E}},e_{i} =\text{LSTM}(x_{i},h_{i-1}^{\text{E}};\ \theta^{\text{E}})\quad \forall i=1,...,n \tag{15}\] \[h_{(i+1)}^{\text{D}},d_{(i)} =\text{LSTM}(x_{(i)},h_{(i)}^{\text{D}};\ \theta^{\text{D}})\quad \forall i=0,1,...,n \tag{16}\]

The decoder output vector \(d_{(i)}\) are, then feed into a fully-connected (FC) layer to calculate probability of the next stop:

\[g_{(i)}=\text{FC}(d_{(i)};\ \theta^{\text{FC}})\quad\forall i=1,...,n \tag{17}\] \[\mathbb{P}(c_{i+1}\ |\ c_{1},...,c_{i},\mathcal{S},X^{\mathcal{S}} ;\theta)=\text{Softmax}(g_{(i)})\quad\forall i=1,...,n \tag{18}\]

where \(g_{(i)}\in\mathbb{R}^{K_{z}}\), \(K_{z}\) is the maximum number of zones in the dataset. And the next predicted stop is selected by maximizing \(\mathbb{P}(c_{i+1}=j\ |\ c_{1},...,c_{i},\mathcal{S},X^{\mathcal{S}};\theta)\) for all \(s_{j}\in\mathcal{S}\setminus\mathcal{S}_{(i)}^{\text{V}}\) (i.e., the zones that are not in the route and that have been visited are excluded).

**Original Pointer Network**. Another benchmark model is the original pointer network (Pnt Net) proposed by (Vinyals et al., 2015). The overall architecture of the pointer network is similar to the proposed model

Figure 8: Model architecture of the LSTM-E-D seq2seq prediction model.
except for the attention component. Specifically, the pointer network calculates attention as:

\[u_{(i)}^{j} =W_{1}^{T}\text{tanh}(W_{2}e_{j}+W_{3}d_{(i)})\quad\forall i,j=1,...,n \tag{19}\] \[a_{(i)}^{j} =\frac{\exp(u_{(i)}^{j})}{\sum_{j^{\prime}=1}^{n}\exp(u_{(i)}^{j^ {\prime}})}\quad\forall i,j=1,...,n \tag{20}\]

The original pointer network does not include the pair-wise local information \((z_{(i)}^{j},\ \phi(x_{(i)},x_{j}))\), and the attention calculation is only quantified from three learnable parameters \(W_{1},W_{2}\), and \(W_{3}\), which may limit its capacity in prediction. We observe that the original pointer network without local information performs extremely badly. For a fair comparison, we add the local information with the similar format in Eq. 19 as:

\[u_{(i)}^{j}=W_{1}^{T}\text{tanh}(W_{2}e_{j}+W_{3}d_{(i)})+W_{4}\begin{bmatrix} z_{(i)}^{j}\\ \phi(x_{(i)},x_{j})\end{bmatrix}\quad\forall i,j=1,...,n \tag{21}\]

After training the model, we generate the final sequence with the greedy algorithm and Algorithm 1, respectively.

### Results

#### 5.4.1 Model comparison.

Table 1 presents the performance of different models. Note that for all approaches except for the TSP, we generate sequences based on two different methods (greedy and Algorithm 1) for comparison. The standard deviation of disparity scores is taken over all testing routes. Results show that sequence generation with Algorithm 1 (i.e., iterating different first zones) can consistently reduce the disparity score for all machine learning methods. It implies that the first zone prediction and the global view (i.e., shortest path) are important for estimating the driver's trajectory.

The proposed method outperforms all other models, both in disparity scores and prediction accuracy. This means the proposed pair-wise ASNN-based attention (Eq. 6) has better performance than the original content-based attention (Eq. 21). The comparison between LSTM-E-D and Pnt Net models demonstrates the effectiveness of the attention mechanism. All machine learning models except for LSTM-E-D can outperform the baseline TSP sequence with Algorithm 1 sequence generation method, suggesting that the hidden trajectory patterns can be learned from the data.

Another observation is that, the prediction accuracy and disparity score do not always move in the same direction. For example, the LSTM-E-D model with Algorithm 1 sequence generation, though has lower accuracy, shows a better disparity score. This is because the accuracy metric does not differentiate "how wrong an erroneous prediction is". By the definition of disparity score, if a stop is \(s_{i}\) but the prediction is \(s_{j}\), and \(s_{j}\) and \(s_{i}\) are geographically close to each other, the score does not worsen too much. This suggests a future research direction in using disparity score as the loss function (e.g., training by RL) instead of cross-entropy loss.

Figure 9 shows the distribution of disparity scores for our proposed method with Algorithm 1 sequence generation (i.e., the best model). We observe that the prediction performance varies a lot across different routes. There is a huge proportion of routes with very small disparity scores (less than 0.01). The mean score is impacted by outlier routes. The median score is 0.0340, which is smaller than the mean value.
#### 5.4.2 Factors on trajectory predictability.

As our proposed model exhibits various levels of predictability across different routes, we aim to investigate which attributes of a route cause high (or low) predictability. This can be done by running a regression model with the disparity score as the dependent variable and route attributes (e.g., locations, departure time, package numbers) as independent variables. The variables used are defined as follows:

* Total planned service time: The estimated time to deliver all packages in the route (service time only, excluding travel time).
* Earliest time window constraint: The earliest due time to deliver packages with time window constraint minus the vehicle departure time. The smaller the value, the tighter the time limit.
* Avg. # traffic signals: Average number of traffic signals in each zone of the route (obtained from OpenStreetMap data).

\begin{table}
\begin{tabular}{c c c c c c c} \hline \hline \multirow{2}{*}{Sequence generation} & \multirow{2}{*}{Model} & \multicolumn{2}{c}{Disparity score} & \multicolumn{3}{c}{Prediction accuracy} \\ \cline{3-8}  & & Mean & Std. Dev & 1st zone & 2nd zone & 3rd zone & 4th zone \\ \hline \multirow{2}{*}{-} & TSP & 0.0443 & 0.0289 & 0.207 & 0.185 & 0.163 & 0.168 \\ \hline \multirow{4}{*}{Greedy} & ASNN & 0.0470 & 0.0289 & 0.150 & 0.141 & 0.119 & 0.123 \\  & LSTM-E-D & 0.0503 & 0.0313 & 0.207 & 0.183 & 0.161 & 0.166 \\  & Pnt Net & 0.0460 & 0.0309 & 0.224 & 0.204 & 0.186 & 0.165 \\  & Ours & **0.0417** & 0.0306 & **0.241** & **0.231** & **0.224** & **0.221** \\ \hline \multirow{4}{*}{Algorithm 1} & ASNN & 0.0429 & 0.0299 & 0.221 & 0.213 & 0.203 & 0.195 \\  & LSTM-E-D & 0.0501 & 0.0305 & 0.182 & 0.156 & 0.142 & 0.149 \\  & Pnt Net & 0.0382 & 0.0301 & 0.286 & 0.273 & 0.262 & 0.274 \\  & Ours & **0.0369** & 0.0301 & **0.320** & **0.310** & **0.303** & **0.314** \\ \hline \hline \end{tabular}
\end{table}
Table 1: Model performance

Figure 9: Disparity score distribution of the best model
* If high-quality route: A dummy variable indicating whether the route is labeled as "high quality" by Amazon or not (Yes = 1). High quality means the actual travel time of the route is similar to or better than Amazon's expectation.
* If in Location: A dummy variable indicating whether the route is in a specific city or not (Yes = 1).
* If departure Time: A dummy variable indicating the (local) departure time (e.g., before 7AM, after 10AM).

Table 2 shows the results of the regression. Since the dependent variable is _disparity scores_, a negative sign indicates a positive impact on the predictability. We observe that routes with tighter time window constraints and more stops are easier to predict. This may be due to the fact that these routes are usually harder to deliver. Hence, to avoid the risk of violating time constraints or delay, drivers tend to follow the planned routes and thus the route sequences are easier predict. We also find that routes associated with larger vans (i.e., larger vehicle capacity) are more predictable. The reason may be that larger vans are less flexible in choosing different routes, thus drivers are more likely to follow the navigation. Another important factor for better predictability is high-quality routes. This may be because high-quality routes are closer to the TSP sequence which we use as inputs. Finally, routes in LA are more predictable than in other areas such as Chicago and Boston.

#### 5.4.3 Impact of input sequence

All machine learning models in Table 1 (except for ASNN) have the LSTM encoder component, which requires the specification of input zone sequence. As mentioned in Section 5.2, we currently use the TSP sequence as input. It is worth exploring the model performance if we use a random zone sequence instead, which corresponds to the scenario without planned route information. Table 3 shows the model performance without the TSP sequence information. Since the ASNN result does not rely on TSP information, it is not listed in the table. Results show that the LSTM-E-D model becomes much worse with a random sequence as inputs, while the performance of Pnt Net and our method is only slightly affected. Even without the planned route information, the proposed model can still provide a reasonable estimation of driver trajectories.

\begin{table}
\begin{tabular}{l l l l} \hline \hline Variables & Coefficients (\(\times 10^{-3}\)) & Variables & Coefficients (\(\times 10^{-3}\)) \\ \hline Intercept & 91.07 ** & If high quality route & -1.66\(\times 10^{-14}\) ** \\ Total \# of packages & 0.059 & If in LA & -4.998 * \\ Total planned service time & -0.476 & If in Chicago & 0.783 \\ Earliest time window constraint & -3.047 ** & If in Boston & -3.354 \\ Avg. \# traffic signals & -3.255 & If on weekends & 1.775 \\ Total \# of stops & -0.142 ** & If departure before 7AM & 0.582 \\ Vehicle capacity (\(m^{3}\)) & -6.041 * & If departure after 10AM & -2.704 \\ \hline Number of routes: 1,002. & & & \\ \(R^{2}\): 0.065; & & & \\ \({}^{**}\): \(p\)-value \(<0.01\); \({}^{*}\): \(p\)-value \(<0.05\). & & \\ \hline \hline \end{tabular}
\end{table}
Table 2: Factors on trajectory predictability
### Summary

Our numerical results show that our proposed model outperforms its benchmarks in terms of disparity scores and prediction accuracy, meaning that it can better predict the actual route trajectories taken by drivers. The comparison with benchmark models shows that our proposed ASNN-based pair-wise attention mechanism and our sequence generation algorithm (Algorithm 1) are both helpful for the prediction. Moreover, we can observe that the predictive performance varies across different routes. Factors such as route quality, delivery time windows, and the total number of stops of a route affect predictability. Finally, the proposed model is insensitive to the input sequence. The prediction performance only slightly decreases when the input sequence is changed from the TSP solution to a random stop sequence. This property implies that we only need the set of stops to implement the model and obtain high-quality solution, while information on the planned route sequence is not strictly required.

## 6 Conclusion and Future Research

In this paper, we propose a pair-wise attention-based pointer neural network that predicts actual driver trajectories on last-mile delivery routes for given sets of delivery stops. Compared to previously proposed pointer networks, this study leverages a new alternative specific neural network-based attention mechanism to incorporate pair-wise local information (such as relative distances and locations of stops) for the attention calculation. To better capture the global efficiency of a route in terms of operational cost (i.e., total travel time), we further propose a new sequence generation algorithm that finds the lowest-cost route sequence by iterating through different first stops.

We apply our proposed method to a large set of real operational route data provided by the Amazon Last-Mile Routing Research Challenge in 2021. The results show that our proposed method can outperform a wide range of benchmark models in terms of both the disparity score and prediction accuracy, meaning that the predicted route sequence is closer to the actual sequence executed by drivers. Compared to the best benchmark model (original pointer network), our method reduces the disparity score from 0.0382 to 0.0369, and increases the average prediction accuracy of the first four zones from 0.229 to 0.312. Moreover, our proposed sequence generation method can consistently improve the prediction performance for all models. The disparity scores are reduced by 10-20% across different models. Lastly, we show that the proposed methodology is robust against changes in the input sequence pattern. Compared to an optimal TSP solution as the input sequence, a random input sequence only slightly increases the disparity score from 0.0369 to 0.0376.

\begin{table}
\begin{tabular}{c c c c c c c c} \hline \hline \multirow{2}{*}{Sequence generation} & \multirow{2}{*}{Model} & \multicolumn{2}{c}{Disparity score} & \multicolumn{3}{c}{Prediction accuracy} \\ \cline{3-8}  & & Mean & Std. Dev & 1st zone & 2nd zone & 3rd zone & 4th zone \\ \hline \multirow{3}{*}{Greedy} & LSTM-E-D & 0.1176 & 0.0498 & 0.045 & 0.047 & 0.041 & 0.050 \\  & Pnt Net & 0.0512 & 0.0323 & 0.090 & 0.096 & 0.097 & 0.096 \\  & Ours & **0.0426** & 0.0311 & **0.204** & **0.192** & **0.195** & **0.196** \\ \hline \multirow{3}{*}{Algorithm 1} & LSTM-E-D & 0.1054 & 0.0463 & 0.103 & 0.061 & 0.049 & 0.052 \\  & Pnt Net & 0.0398 & 0.0311 & 0.298 & 0.284 & 0.273 & 0.273 \\ \cline{1-1}  & Ours & **0.0376** & 0.0307 & **0.316** & **0.298** & **0.302** & **0.298** \\ \hline \hline \end{tabular}
\end{table}
Table 3: Model performance without TSP information
The data-driven route planning method proposed in this paper has several highly relevant practical implications. First, our proposed model performs well at predicting stop sequences that would be preferable to delivery drivers in a real operational environment, even if it is not provided with a theoretically optimal (i.e., minimal route duration) planned TSP sequence as an input. Therefore, the model can be used to generate a predicted actual stop sequence that a driver would likely be taking for a given set of delivery stops. The prediction can serve as a new 'empirical' planned route that is informed by historical driver behavior and thus more consistent with the driver's experience and preferences. Second, by comparing the stop sequence predicted by our model with the traditional, TSP-based planned stop sequence, a route planner may infer potential reasons for the drivers' deviations and adjust the company's planning procedures and/or driver incentives if necessary. Third, as stop sequence generation using machine learning models is computationally more efficient than traditional optimization-based approaches, a trained machine learning model can be applied in real-time to quickly re-optimize routes when drivers are unexpectedly forced to deviate from their original stop sequence (e.g., due to road closures) and need updated routing strategies.

Based on the work presented in this paper, a number of fruitful future research avenues arise. First, instead of focusing on stop sequence prediction, future work may improve the interpretability of such prediction models and develop machine learning approaches that better explain which factors cause drivers to deviate from a planned stop sequence and how they affect their actual route trajectories. Second, future work should attempt to combine the strengths of optimization-based route planning approaches and machine learning by incorporating tacit driver knowledge learned via machine learning models into route optimization algorithms.

## Appendix A Mathematical Formulation of a LSTM Cell

The details of an LSTM cell, \(h_{t},e_{t}=\text{LSTM}(x_{t},h_{t-1};\ \theta)\), is shown below:

\[f_{t} =\sigma_{g}(W_{f}x_{t}+U_{f}h_{t-1}+b_{f})\] (A.1) \[i_{t} =\sigma_{g}(W_{i}x_{t}+U_{i}h_{t-1}+b_{i})\] (A.2) \[o_{t} =\sigma_{g}(W_{o}x_{t}+U_{o}h_{t-1}+b_{o})\] (A.3) \[\tilde{c}_{t} =\sigma_{c}(W_{c}x_{t}+U_{c}h_{t-1}+b_{c})\] (A.4) \[c_{t} =f_{t}\circ c_{t-1}+i_{t}\circ\tilde{c}_{t}\] (A.5) \[h_{t} =o_{t}\circ\sigma_{h}(c_{t})\] (A.6) \[e_{t} =h_{t}\text{ (if this is a single layer one-directional LSTM)}\] (A.7)

where \([W_{f},W_{i},W_{o},W_{c},U_{f},U_{i},U_{o},U_{c},b_{f},b_{i},b_{o},b_{c}]=\theta\) is the vector of learnable parameters. \(x_{t}\) is the input vector to the LSTM unit. \(f_{t}\) is the forget gate's activation vector. \(i_{t}\) is the input/update gate's activation vector. \(o_{t}\) is the output gate's activation vector. \(h_{t}\) is the hidden state vector. \(e_{t}\) is the output vector of the LSTM. Note that for a multi-layer or bidirectional LSTM, \(e_{t}\) may not equal to \(h_{t}\). In this study, we use a single layer one-directional LSTM and thus have \(e_{t}=h_{t}\). More details on the output vector can be found in Pytorch (2021). \(\tilde{c}_{t}\) is the cell input activation vector. \(c_{t}\) is the cell state vector. "\(\circ\)" indicates the component-wise multiplication.
## Appendix B From Zone Sequence to Stop Sequence

The complete stop sequence is generated based on the given zone sequence. The detailed generation process is shown in Algorithm 2.

```
1:functionCompleteSeqGeneration(\((\hat{z}_{(1)},..,\hat{z}_{(n)}),\{\mathcal{S}_{(i)},i=1,,,n\}\))
2:\(s_{prev}\gets D^{\mathcal{S}}\)
3:\(s_{complete}^{*}\leftarrow(s_{prev})\)\(\triangleright\) Initialize the complete stop sequence with depot
4:for\(i\in\{1,...,n-1\}\)do
5:\(\mathcal{S}_{first}\leftarrow\) Set of three stops in \(\mathcal{S}_{(i)}\) that are closest to \(s_{prev}\)
6:\(\mathcal{S}_{last}\leftarrow\) Set of three stops in \(\mathcal{S}_{(i)}\) that are closest to all stops in \(\mathcal{S}_{i+1}\) on average
7:\(\mathcal{P}_{(i)}\leftarrow\emptyset\)\(\triangleright\) Initialize the set of optimal paths in zone \(\hat{z}_{(i)}\)
8:for\(s_{first}\in\mathcal{S}_{first}\)do
9:for\(s_{last}\in\mathcal{S}_{last}\)do
10:if\(s_{first}=s_{last}\)then
11:\(\hat{p}_{\text{temp}},t_{\text{temp}}=\texttt{T}\texttt{P}\texttt{P}( \mathcal{S}_{(i)})\)\(\triangleright\) Solve the optimal tour and travel time for zone \(\hat{z}_{(i)}\)
12: Delete the last edge back to \(s_{first}\) in the tour \(\hat{p}^{\text{temp}}\). Let the new path and travel time be \(\hat{p}^{\prime}_{\text{temp}}\) and \(t^{\prime}_{\text{temp}}\)
13: Add \(\hat{p}^{\prime}_{\text{temp}}\) and \(t^{\prime}_{\text{temp}}\) to \(\mathcal{P}_{(i)}\)
14:else
15:\(\hat{p}_{\text{temp}},t_{\text{temp}}=\texttt{P}\texttt{P}\texttt{P}\texttt{P} (\mathcal{S}_{(i)},s_{first},s_{last})\)\(\triangleright\) Solve the optimal path and travel time for zone \(i\)
16: Add \(\hat{p}_{\text{temp}}\) and \(t_{\text{temp}}\) to \(\mathcal{P}_{(i)}\)
17:\(\hat{p}_{(i)}\leftarrow\) Path in \(\mathcal{P}_{(i)}\) with the minimum travel time
18:\(s_{complete}^{*}\leftarrow(s_{complete}^{*},\ \hat{p}_{(i)})\)\(\triangleright\) Concatenate two sequence
19:\(s_{prev}\leftarrow\) Last stop of path \(\hat{p}_{(i)}\)
20:\(s_{complete}^{*}\leftarrow(s_{complete}^{*},\ D^{\mathcal{S}})\)\(\triangleright\) Concatenate the last stop as the depot
21:return\(s_{complete}^{*}\)
```

**Algorithm 2** Complete sequence generation. Input: zone sequence \((\hat{z}_{(1)},..,\hat{z}_{(n)})\), depot \(D^{\mathcal{S}}\), set of stops in each zone \(\mathcal{S}_{(i)}\), \(i=1,...,n\). \(\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P} \texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P} \texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P} \texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P} \texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P} \texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P} \texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P} \texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P} \texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P} \texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P} \texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P} \texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P} \texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P} \texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P} \texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P} \texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P} \texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P} \texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P} \texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P} \texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P} \texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P} \texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P} \texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P} \texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P} \texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P} \texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P} \texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P} \texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P} \texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P} \texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P} \texttt{P}\texttt{P}\texttt{P}\texttt{P}{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P} {P}\texttt{P}\texttt{P}{P}\texttt{P}\texttt{P}{P}\texttt{P}\texttt{P}{P}\texttt{P} \texttt{P}{P}\texttt{P}{P}\texttt{P}{P}\texttt{P}\texttt{P}{P}\texttt{P}{P} \texttt{P}{P}\texttt{P}{P}\texttt{P}{P}\texttt{P}{P}\texttt{P}{P}\texttt{P} {P}\texttt{P}{P}\texttt{P}{P}\texttt{P}{P}\texttt{P}{P}\texttt{P}{P}\texttt{P} {P}\texttt{P}{P}\texttt{P}{P}\texttt{P}{P}\texttt{P}{P}\texttt{P}{P}\texttt{P} {P}\texttt{P}{P}\texttt{P}{P}\texttt{P}{P}\texttt{P}{P}\texttt{P}{P}\texttt{P} {P}\texttt{P}{P}\texttt{P}{P}\texttt{P}{P}\texttt{P}{P}\texttt{P}{P}\texttt{P} {P}\texttt{P}{P}\texttt{P}{P}\texttt{P}{P}\texttt{P}{P}\texttt{P}{P}\texttt{ {P}}\texttt{P}{P}\texttt{P}{P}\texttt{P}{P}\texttt{P}\texttt{P}{P}\texttt{ {}}\texttt{P}{P}\texttt{P}{P}\texttt{P}{P}\texttt{P}{P}\texttt{P} {}\texttt{P}{P}\texttt{P}{P}\texttt{P}{P}\texttt{}{P}\texttt{P}{P}\texttt{ {}}\texttt{P}{P}\texttt{P}{P}\texttt{P}{P}\texttt{P}{P}\texttt{P} {P}\texttt{P}{P}\texttt{P}\texttt{P}{P}\texttt{}\texttt{P}{P}\texttt{P} {}\texttt{P}{P}\texttt{P}\texttt{P}{P}\texttt{}\texttt{P}{P}\texttt{P} {P}\texttt{}\texttt{P}{P}\texttt{P}\texttt{}}{P\texttt{}\texttt{P}{P}\texttt{} {P}\texttt{P}\texttt{P}\texttt{}P{}\texttt{}\texttt{P}}{P}\texttt{}\texttt{ {P}\texttt{P}\texttt{P}}{P}\texttt{}\texttt{P}\texttt{}\texttt{P}{P}\texttt{} \texttt{P}\texttt{}\texttt{P}\texttt{}\texttt{P}\texttt{}\texttt{P}\texttt{} {P}\texttt{P}\texttt{}\texttt{P}\texttt{}\texttt{P}\texttt{}\texttt{P}{P}\texttt{} {P}\texttt{P}\texttt{P}\texttt{}\texttt{P}\texttt{}\texttt{P}\texttt{P}\texttt{} {P}\texttt{P}\texttt{P}\texttt{}\texttt{P}\texttt{P}\texttt{}\texttt{P}\texttt{} {P}\texttt{P}\texttt{P}\texttt{}\texttt{P}\texttt{P}\texttt{}\texttt{P}\texttt{P} \texttt{P}\texttt{P}\texttt{}\texttt{P}\texttt{P}\texttt{P}\texttt{}\texttt{P}\texttt{P} \texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P}\texttt{P} \texttt{P
nine small-scale path TSP problems will be solved since there might be overlapping between the first and the last stops. If the first and the last stops are identical, we solve a tour TSP problem and output the path by deleting the last edge which traverses back to the first stop in the tour.

After having all potential inner zone paths and total path travel time between any first and last stop pair, we keep the path with the minimum travel time as the inner zone sequence, \(\hat{p}_{(i)}\). The key assumption we make here about drivers is that they will deliver packages within a zone following a path that minimizes their total travel time. With the optimal inner zone stop sequence of the current zone, we then move to the next visited zone in the optimal zone sequence and repeat the same procedure until we generate the complete stop sequence.
