[ad_1]

Chen Wu, Yin Song, Baichuan Sun, Eden Duthie

@AWS

In this blog post, we discuss our open source solution to solving the district optimisation problem for last mile planning. This is a collaboration between AWS and Aramex Australia. Aramex Australia operates a courier service under a franchise model. Its logistics network includes 29 regional franchises and over 900 franchise partners across the country. Each regional franchise operates a delivery station over a territory. Currently couriers are assigned packages based on predefined geo-regions (e.g. postcode areas).

Aramex Australis is exploring alternative allocation methods to improve their driver experience and delivery efficiency. AWS had developed new innovations in last-mile routing and sequencing which Aramex Australia was interested in. Aramex Australia thus worked with AWS to solve this problem.

Our method used public data from the Australian government, supervised learning and the graph partitioning algorithm (METIS). The result suggests that our method could improve driver workload balance and shorten travel time. Moreover, it provides a quantitative tool for last mile planners to tinker with the tradeoff between delivery efficiency and driver workload.

Package allocation can be formulated as a districting problem — we are tasked to assemble Basic Geographic Units (BGU), such as streets, post code areas, suburbs (of residents), into larger geographic districts, in order to achieve some optimisation objectives or planning criteria.

In the case of package allocation, these artificially constructed geographic districts are known as “polygons”, and packages are automatically allocated to different polygons based on the spatial matching between their addresses and the boundaries (geo-fence) of the established polygons. More specifically, the optimisation objectives may involve reducing travel time or cost, keep the delivery workload balanced, etc, and the planning criteria for package allocation often boil down to improving delivery efficiency.

In this work, we used Statistical Areas (defined in Section 3.1) as the BGUs. Each geo-unit contains a collection of residential addresses that represent the demand for package deliveries. Multiple geo-units (of the same colour as shown in Figure 0) are grouped into a District, which becomes a *Polygon* to be managed by a single courier.

The key here is to determine which BGUs belong to the same district, i.e., the colour — green, purple, or grey — of each BGU in Figure 0. Once a district is formed, subjects (delivery addresses) from all BGUs in that district are combined as a whole assigned to a given courier for delivery. Since it is not always possible for Aramex to directly involve couriers (e.g. third-party logistics providers) operations such as route sequence and optimisation, the only decision variable that Aramex Australia can optimise is the “colour” of those BGUs, namely, how to form larger district clusters out of smaller BGUs in order to achieve delivery planning objectives and criteria.

## 3.1 Demand Sampling and Basic Geo-Units

We obtained the publicly available Geoscape Geocoded National Address File released by the Australia Government. From those 4.7 million *property* and *building *addresses, we randomly sampled a set ** S** of 4,968

*addresses covered by the Aramex Sydney south west delivery station. No Aramex customer addresses were used in this analysis. Instead, we used geocodes associated with public addresses to simulate the distribution of customer demand in an Aramex delivery region.*

We used Statistical Areas (SA), defined in the Australian Statistical Geography Standard Edition 3 (ASGS) as the basic geo-units** **for polygon optimisation. The main motivation for adopting this standard is that ASGS boundaries are

- well in line with existing land parcel and property boundaries and road networks.
- updated regularly as per changes to dwellings and population to ensure their accuracies

Specifically, we obtained the digital boundary files — Statistical Areas Level 2 (SA2) — released by the Australian Bureau of Statistics. There are in total 2,473 SA2 across Australia, and we only selected a small set ** P** of 62 SA2s such that each SA2 in

**covers or intersects with at least one public address in**

*P***. Figure 1a shows the addresses (green dots) in**

*S***overlaid on the SA2s in**

*S***. It is clear that the address density varies across different SA2s. This is expected since the populations in some regions are relatively sparse. While an SA2 generally has an average of 10,000 people, it also considers many other criteria when determining its boundary such as functional areas, region cohesiveness, road networks, community facilities, boundary stability, etc. For this reason, a single SA2 may consist of several smaller suburbs and vice versa, a spatially large suburb could be decomposed into several SA2s.**

*P*Aramex Australia proposed that the new polygon schemes be similar to the current polygon scheme. This is to minimise the disruption to current driver affinity program and package allocation scheme while still benefiting from reduced cost such as travel time. This requirement has motivated us to use SA2 as the Basic Geo-Units. Figure 1b overlays the current Aramex polygons onto SA2s. It can be seen that many Aramex polygons share the same boundaries with SA2s, and the size of the Aramex polygons are either comparable to or slightly smaller than an SA2. This suggests that merging SA2 to form new polygons involves less dramatic changes than solutions based on, for example, SA1. Figure 2 shows that the average size of an SA1 is much smaller than Aramex polygons. When merging many SA1 to form the new polygon schemes, the degree of freedom — the number of ways these SA1s can be combined in order to achieve the same level of optimal objective (e.g. travel time) — is very high. Consequently, it would be much more challenging for the optimisation process to meet the “similarity” constraint without explicit modelling efforts if the SA1 areas were used.

## 3.2 Problem reformulation

The goal of districting is to reduce the total travel time required to fulfil the delivery. Assuming the number *k* of drivers and the total delivery quantity is constant, in order to reduce the travel time, we can cast the districting problem as a graph partition problem. Given an undirected graph *G* = (*P*, *E*), where each edge *r* ∈ *E* represents the actual travelling route between a pair of SA2s in *P*, we define travel cost *w*** **as the non-negative weights associated with

*r*,

*w*:

*E*→ R. The optimisation goal is to partition all SA2s ∈

*P*into an optimal set

*Π*of polygons {

*B1*,

*B2*, …

*Bk*}, where

*B1*U · · · U

*Bk*=

*P*,

*Bi*∩

*Bj*= ∅ ∀

*i*

**≠**

*j*, such that the sum of weights associated with all

*cut*edges that span between polygons is maximised.

Π* =argmaxΠΣijw(Eij:= {{u,v} ∈E| u ∈Bi,v∈Bj, i≠j} — Eq.1

where *Ei,j* is the **set** of edges that straddle between polygon *Bi *and *Bj*. This formulation suggests that if a pair of SA2s has a large weight (e.g. a high travel cost), they will be placed in different polygons so that their edge can be included in some *Ei,j *in order to increase *w*(*Eij*), thereby maximising the weight sum.* *Likewise, more importantly, SA2s that are close to one another (measured as a low travel cost) will be placed in the same polygon *Bk. *In practice, this means packages whose addresses fall inside any SA2s ∈ *Bk* will be allocated to the same courier for delivery.

Here maximising the total weights of all inter-polygon edges is somewhat equivalent to minimising the accumulated weights of sequential edges inside polygons, which is shown as minimising the total travel time in Figure 0. This is because the total weights of all edges in *E*, which is a sum over all inter- and intra-polygon weights, is a constant. So maximising the former leads to minimising the latter.

Aramex Australia currently uses postal codes to determine the current scheme *H* of polygons {*BH1*, *BH2*, … *BHk*}, where each *BHi* roughly represents a post code area. Section 4 will compare the delivery efficiency under the current polygon scheme *H* and the proposed scheme *Π *respectively.

## 3.3 Distance Matrix Generation

This step produces the input matrix for graph *G* defined in Section 3.2. To start with, we used docker-compose to configure and launch a local Docker container that runs the Valhalla routing service API on an Amazon SageMaker notebook instance. For the data layer, we used the latest Australia OpenStreetMap geographic data file publicly hosted on GEOFABRIK.

Next, we concatenated a list *L* of addresses (namely, their locations expressed as latitude and longitude) selected from ** S** as a batch (e.g. with a size of

*N*= 50, 100 or 200), and used

*L*as both “sources” and “targets” parameters to query the Valhalla Matrix API. The query response is a

*N*x

*N*matrix, and each entry of this matrix contains the both the travel distance and travel time between a pair of

*source-to-target*addresses. Let

*M ( = 4968)*be the length of

*S*, and since

*N*<<

*M*, it takes 600 to 10,000 API invocations to fill up the entire large

*M*by

*M*distance matrix

**.**

*D*Since most optimisation solvers solve a cost minimisation problem (e.g. minimum cut). In order to use these solvers to solve Eq. 1, which represents a cost maximisation formulation, we inverted the cost (time and distance) returned from the matrix API as *new_cost* = *max_cost *— *raw_cost**, *where *max_cost* is the maximum value of ** D**. The

*new_cost*can be interpreted as a measure of

**or adjacency.**

*closeness*Once the distance matrix is generated and saved, it can be re-used repeatedly for the following optimisation steps. In practice, we recommend that the matrix be re-calculated periodically depending on the region and seasonality. There will be multiple matrices generated for all regions of interest.

## 3.4 Statistical Area Graph Generation

This step creates the graph *G* defined in Eq 1. in Section 3.2 where each node of the graph is an SA2. Two inputs for this step are (1) the set of SA2s to form the node set ** P**, and (2) the distance matrix

**generated in Section 3.3 to form the edge set**

*D***. Entries in**

*E***need to be aggregated in order to form the weights associated with each**

*D**e*∈

*E**.*This is because each entry

*Dij*represents the travel time/distance from address

*i*to address

*j*, where

*i*and

*j*∈ S. We need to obtain the weight

*w*(r) between each pair of

*SA2a*and

*SA2b*instead. For example, if

*SA2a*contains 10 addresses and

*SA2b*contains 80 addresses, there will be 80 x 10 = 800 pairs of addresses between

*SA2a*and

*SA2b*. To account for these multiple pairs of addresses, we computed the average new_cost of all of these 800 address pairs using

**. This can be achieved by extracting all entries in**

*D***whose row identifiers correspond to all addresses ∈**

*D**SA2a*and whose columns represent all addresses ∈

*SA2b*, thus producing a sub-matrix

*Dab*. We then calculated the mean value of all entries in

*Dab*to obtain the weight associated with pair

*SA2a*—

*SA2b*.

We repeated the mean value calculation for each pair of SA2s and produced a 62 x 62 small matrix ** D’**, which is essentially the weighted adjacency matrix of the graph

*G*that we wish to partition. Unlike a normal adjacency matrix where each entry represents a binary (namely 1 or 0) adjacency relationship between a pair of graph nodes, each entry of

**represents a continuous adjacency relationship. A larger adjacency value corresponds to a shorter distance between the pair of SA2 nodes. Moreover, we defined a threshold**

*D’**T*such that any adjacency values less than

*T*are re-set to 0, rendering them non-adjacent pairs. This converts

*G*from a fully-connected graph to a partially-connected graph, where only “close-enough” graph nodes have edges running between them. In this step, we set

*T*to be the 99 percentile of all values in

**’ in order to form a sparse graph.**

*D*## 3.5 Graph Partitioning

This step uses the METIS Family of Graph and Hypergraph Partitioning Software to perform graph partitioning on *G*. We applied the multilevel k-way partitioning algorithms (** PartGraphRecursive**) to compute a partitioning solution as per Eq 1. The main objective, as discussed in Section 3.2, is to maximise (or minimise) the travel time (or closeness) between SA2s that are assigned to different polygons. Moreover, we added an additional constraint in order to balance the number of addresses per polygon. This is to ensure the driver workloads across different polygons are similar. This is achieved by (1) setting the node weight as the number of addresses in its corresponding SA2, and (2) passing the node weight list as an extra parameter to the

*PartGraphRecursive*algorithm.

In summary, the input of *PartGraphRecursive* includes (1) the weighted adjacency list of the sparse matrix ** D’**, (2) the node weight for preserving workload balance, and (3) the hyper-parameter — number

*O*of output partitions (i.e. polygons). We initially set

*O*to 27, which was the current number of polygons. The output of

*PartGraphRecursive*includes (1) the total edge cut, namely the total degree of closeness across all polygons, and (2) the partition index for each SA2 ∈ P in

*G.*

Figure 3 (Left) shows the set of 27 polygons produced by *PartGraphRecursive*. Each polygon is given its unique colour, and is composed of one or many SA2s. For polygons with a relatively smaller area, its address density (the darkness) is higher. On the contrary, larger polygons tend to have more spread-out address distributions. This is the consequence of the added constraint that keeps the number of service addresses balanced.

Does our solution work? For driver experience, we want to know if the package demand can be evenly distributed to drivers. For delivery efficiency, we check if the delivery throughput can be further improved under the new polygon solution given that all other input conditions are identical.

## 4.1 Demand balance

This section examines if the total demand is well balanced across multiple polygons in order to maintain a good driver experience.

Figure 4 shows the number of addresses for each one of the 27 polygons. It can be seen that most of the “new” polygons have around 150 addresses. Their distribution appears more uniform compared to that under the current polygon scheme on the right. Concretely, we measure the load balance (distribution uniformity) using the standard deviation (STD) of the number of addresses across all polygons in both schemes. A smaller STD corresponds to a better demand balance. The STD of the current polygon scheme is 69, and the STD of the proposed polygon solution is 37. This suggests that the new polygon scheme achieved an increase of 46% in demand balance compared to the current scheme.

During the graph partitioning step in Section 3.5, we set the node weight as the number of addresses for each SA2 node. The *PartGraphRecursive* algorithm then took as a constraint the node weight list and strived to keep the total weight uniformly distributed across all partitioned polygons.

## 4.2 Travel time

In order to quantitatively evaluate delivery efficiency, we ran a set of TSP route sequencing tasks in order to simulate delivering 4,968 packages under the two different polygon schemes. To obtain robust results, we developed two “simulation” methods to evaluate these two polygon schemes.

In the first simulation, we took as input the demand (i.e. delivery addresses) covered by a given polygon. We used the Google OR-Tools TSP solver to generate a route sequence based on pair-wise travel time estimates stored in the distance matrix *D *(see Section 3.3). The TSP solver uses some (re-configurable) local search strategies to find the global or local optimal solutions from an initial solution. It should be noted that our goal during evaluation is not to find optimal route sequence solutions. Rather, it aims to produce “reasonable” route sequences given demands in the Polygons. As long as the SAME route sequencer is consistently used for simulating routes under both polygon schemes, our evaluation result is useful. The OR-Tools TSP solver returned the time for traversing the output sequence within each polygon, and we repeated this TSP sequence generation for each one of the 20 or 27 polygons.

In the second simulation method, we trained a supervised machine learning model to directly estimate the travel time given coordinates associated with a set of delivery addresses. This is different from the “microscopic simulation” used in the first method. There, we had to literally generated the route sequence to calculate the total travel time by summing over the “actual” travel time from one transit to the next, until reaching the last transit and the depot. In the ML-based second method, we first extracted a list of features from the original set of coordinates. While many of these features are based on a previous study (Liu, He, and Shen 2020) on travel time estimation using ML methods, we proposed four additional features associated with the convex hull, and we found them to be the most important features in correctly estimating the travel time.

We then used the PyCaret AutoML library to train and select the best model ExtraTreeRegressor, which has a *R2* score of 0.8857 for regression. We trained the model using a dataset that contains around 9,000 historical route sequences and used as the ground-truth actual travel time recorded from driver scans. Similar to the TSP solver in the OR-tools, our goal is to compare the total travel time under two polygon schemes rather than predicting absolute travel times with a perfect accuracy.

Table 2 presents cost metrics, namely the total travel time. The total travel time is the travel durations summed over all TSP sequences from 27 polygons. Since the same TSP sequence solver is used on the same set of 4,968 addresses, the only difference that could have explained the slight reduction in travel time lies in how these nearly 5,000 addresses have been allocated into the different sets of polygons. Moreover, when we simulated / estimated the total travel time using two unrelated methods (OR-Tools and ML model), both methods showed that the new polygon scheme is superior to the present one in terms of the delivery efficiency.

## 4.3 Delivery efficiency

Thus far, we intentionally fixed the number (27) of polygon as a constant in order to minimise potential disruptions to Aramex Australia operations and fleet management. On the other hand, our solution can help Aramex Australia operators to experiment with the optimal number of polygons in terms of the tradeoff between delivery efficiency and extra driver workload. Both of them can be calculated from the travel time using Eq. 3 and Eq. 4 respectively.

=delivery efficiency# of packages / (total travel time× # of drivers)

=extra workloadtotal travel time / # of drivers—current travel time / 27

Eq. 2

In Eq. 2, the *number of drivers* is the same as the number of polygons, which is the hyper-parameter — number *O* of output partitions — described in Section 3.5. The* total travel time* is calculated using the simulation method described in Section 4.2 (given the number of polygons 27), and is exemplified in Table 2. We set the number *Oi* of polygons to be 27, 26, 25, 24, 23, 22, 21, and 20 in succession and calculated its corresponding total travel time *tti*. For each pair of *tti *and *Oi*, we used Eq. 2 to calculate. Figure 5 presents all metrics collected from these calculations (also shown in the Appendix). It should be noted that the workload only includes travel time.

As we reduce the number of polygons (drivers), the delivery efficiency — number of packages delivered per hour per driver — also increases, which is a positive outcome. However, we notice that the extra workload per driver (measured in minutes) also increases as the number of driver drops, which is a negative outcome but is expected given the total demand is a constant. Figure 5 provides a quantitative tool for the operators to determine a “sweet spot” where the delivery efficiency is improved at the cost of a reasonable level of additional workload (travel time). For example, if the Aramex team let each driver work for additional 15 minutes, they could potentially increase the delivery efficiency by 29%, increasing from 3.1 to 4.0. Such a “what-if” analysis can be performed by our solution at a large scale for several reasons. First, it has a near-optimal graph partitioning algorithm that takes in as the hyper-parameter the number of polygons. Next, it uses simulations for estimating the travel time. Last, the end to end calculation from graph partitioning to sequence simulation can be completed within an order of seconds.

Our current solution benefits from the fact that Australian government agencies have developed these SAs over the years and made them publicly available. As discussed in Section 3, we were able to exploit the fact that SA2 nodes are relatively similar to the existing Aramex Australia polygons. If this was not the case, it would have been much more challenging to achieve a new polygon scheme that has a minimum disruption for drivers.

To tackle this challenge, we are currently experimenting with an approach based on our reinforcement learning framework as shown in Figure 6. The central idea is that we take the current polygon scheme as the ‘initial solution’, and let an agent explore making small changes to existing polygons, and train it to discern positive changes from negative ones through *rewards* as evaluated via metrics such as travel time, etc. An advantage of this new RL-based approach is that it does not necessarily need good quality BGUs such as SA2s. It can evolve from any previous solutions generated during current (best) solutions.

Aramex Australia used the open source solution to pilot a re-arrangement of their polygon and package allocation scheme and have thus far observed an improvement in delivery efficiency based on two extensive simulations. The solution is available here if you would like to try it yourself and can be run entirely with publicly available data. Moreover, if you have any issues to get the tool running or want to tell us what you think, please contact us via GitHub issues. We are happy to support you.

We would like to thank our Aramex Australia collaborators — *Renee Qian*, *Michael Krason*, *Tom Parkinson*, and *Ruby Wolf *— for providing us with valuable insight in last mile planning. We would like to thank our Amazon colleagues — *Manuel Baeuml* and *Jatin Arora *— for providing useful feedback to improve this blog post.

[ad_2]

Source link