Last-mile mail delivery using multimodal green transportation modes

1 Introduction

With the rapid rise of e-commerce and the corresponding surge in industry and academic interests in digital innovation and last-mile parcel delivery, research has predominantly focused on parcel logistics, aiming to reduce last-mile logistics costs and maximize carriers and shippers profits (Sun et al., 2025). However, advances in the effectiveness and convenience of digital communication technologies and platforms have also led to a sharp and sustained decline in traditional mail volumes. The current lack of research focuses on the factors affecting the efficiency of last-mile mail deliveries, presents critical financial challenges for postal carriers and, in many cases, jeopardizes the long-term sustainability of the sector.

Empirical studies indicate that over 90% of people now rely on mobile phones and internet-based platforms for daily communication, significantly reducing the need for personal letters (Shao et al., 2022). Globally, mail-post volumes dropped by approximately 23% between 2015 and 2020 (Universal Postal Union, 2021), with some regions experiencing even steeper declines. In the U.S., the volume of first-class mail fell from 92 billion pieces in financial year 2008 to 46 billion in financial year 2023, a 50% reduction (USPS, 2024). Similarly, the U.K. reported a 9.5% year-over-year decrease in addressed mail deliveries, falling from 8.0 billion in 2021–2022 to 7.3 billion in 2022–2023 (Ofcom, 2023). These figures highlight a systemic transformation in communication habits and underscore the urgency of rethinking the operation of traditional mail services.

Australia Post, the largest carrier in Australia, delivered nearly half a billion parcels in 2022–2023, yet mail volumes have continued to decline, falling by about two-thirds since their 2008 peak (Australian Government, 2023). This sharp volume decrease means that postal carriers are delivering far fewer items per route, even as the delivery network (number of addresses) continues to expand (Crew and Brennan, 2016). Therefore, the distance between delivery points grows, leading to inefficiency and higher cost per mail delivered. Despite this sustained reduction, Australia Post remains committed to maintaining mail delivery and improving social sustainability by reducing incident rates. Replacing motorcycles with electric tricycles for household deliveries to enhance safety and operational efficiency is one of the potential solutions. However, electric tricycles are not viable in all regions due to challenging terrain or high-speed access roads. In these areas, vans are often the only practical option, though they operate with considerably lower productivity, typically 25%–50% less than motorcycles (Australia Post, 2020). Resizing districts, consolidating delivery points, and co-delivering mail and parcels are other suggested solutions (Sandoval et al., 2022). However, the problem remains acute in low-density suburban areas where traditional methods yield low productivity.

Australia Post’s trials indicate that productivity losses are closely linked to local delivery volumes and geographical characteristics (Australia Post, 2020). In low-density areas with minimal parking restrictions, vans can maintain reasonable efficiency by quickly moving between delivery points. Conversely, in dense urban environments with constrained parking, park and loop deliveries, where drivers stop once and deliver to clustered addresses on foot, can outperform van-only approaches. While delivery strategies can be adapted to each of these distinguished situations, this study develops a mail delivery strategy for large geographic regions involving both suburban and rural areas. This heterogeneity necessitates flexible delivery strategies capable of adapting to both high- and low-density operational environments. The current route planning methods employed by Australia Post require improvement and automation, as they are characterized by high complexity, time-consuming processes, heavy reliance on individual planners’ expertise, and a lack of standardization. These limitations hinder the organization’s ability to adapt efficiently and cost-effectively to future network changes. There is therefore a pressing need for an optimizer system to support the last mile network planning for the efficient deployment of vans and integration with other delivery modes.

To address these challenges, the authors, in collaboration with Australia Post, developed a comprehensive Last-Mile Multimodal Delivery Planning and Optimization system. The system is flexible and scalable, integrating vans with green delivery modes, walking and cycling, and optimizing multimodal green vehicle routing problem (MG-VRP) routing and mode allocation based on local conditions, including customer density, parking availability, and mode-specific capacity and speed. The designed multimodal mail delivery system is illustrated in Figure 1, which provides a conceptual view of the different delivery configurations. The workflow spans data ingestion and geospatial preprocessing, followed by optimization powered by adaptive large neighborhood search (ALNS). Within the optimizer system, the ALNS minimizes total costs and delivery time by determining which customers are allocated to direct van service. At the same time, the optimizer identifies customers that can be served by walking or biking. In this case, the van’s optimal parking location is selected, and a closed green-mode delivery loop is constructed beginning and ending at that location. The modeling framework incorporates explicit service times for each customer, representing the on-site time required to complete the delivery task. We consider different service times for van, bike, and walk modes but assume they are uniform across customers served by the same mode. Finally, the optimized route data are transferred to the visualization module of the optimizer system.

The developed optimizer system was implemented in two real-world case studies in Belgrave in the State of Victoria (VIC) and Katoomba in the State of New South Wales (NSW). First, Belgrave was selected as the initial case study region due to its uneven terrain, diverse road network, and more than 6,500 delivery points, providing a unique opportunity for designing a robust solution for the mail delivery network. Katoomba was then selected as the second case study to further validate the optimizer system’s performance and scalability, given its geographical similarity to Belgrave. The Katoomba case study includes more than 20,000 delivery points, illustrating the large-scale nature of the delivery network and demonstrating the practicality of the proposed optimizer system for supporting data-informed decisions in complex last-mile systems. The geographical locations and key statistics of both case studies are presented in Figure 2, while detailed information about the study scope is provided in Section 6.

Despite growing research interest in MG-VRP, a clear gap remains in addressing mail delivery, a segment often overshadowed by parcel focused strategies. Prior studies concentrate on parcel distribution in high density urban settings, often assuming stable or increasing volumes (Kahalimoghadam et al., 2021; Thompson et al., 2024). However, given the steady decline in mail volumes, existing models lack the flexibility and cost efficiency needed to support evolving postal networks. In existing multimodal frameworks, vans typically act as a passive mode, delegating final delivery tasks entirely to walking or bicycle couriers (Martinez-Sykora et al., 2020; Wehbi et al., 2022). This perspective overlooks the operational advantages of using vans not only to deliver mail but also to facilitate walking and cycling delivery. Furthermore, although truck-drone delivery, as a multimodal VRP, has received growing attention in recent literature (Park et al., 2025), it faces substantial barriers to implementation due to limited public acceptance, unresolved legal issues, and persistent safety concerns (Sood et al., 2025). These limitations make drone based last mile solutions unlikely to be feasible in the near future. There remains a notable absence of integrated MG-VRPs in mail delivery. Existing models seldom allocate delivery tasks dynamically across vans, bicycles, and walking modes, nor do they adequately incorporate the implications of declining mail volumes. These limitations form a critical barrier to achieving sustainable and cost-effective mail delivery systems.

To fill these gaps, this study advances the MG-VRP field by developing and implementing an MG-VRP that enhances the economic and environmental sustainability of mail delivery systems under declining mail volumes. Additionally, in the developed model, vans are not only used to facilitate green delivery by walking or biking, but also to deliver letters directly to customers. Finally, the developed letter delivery system does not face public acceptance or legal barriers, as walking and biking delivery do not involve the safety concerns associated with drone delivery. Hence, the main contributions of this study include:

• An MG-VRP is formulated in which mail delivery is performed by vans and two green modes, walking and cycling. The model dynamically determines optimal van parking locations and assigns green delivery modes to clusters of nearby customers within a closed-loop structure.
• An ALNS algorithm is developed to identify optimal van routes, parking locations, and customer clusters best served by green modes, ensuring service time efficiency and CO₂ reduction.
• An optimizer system is developed that unifies data ingestion, route optimization, and visualization within a single platform. The system is scalable and adaptable to evolving demand and operational conditions, enabling delivery planners to design safer, more sustainable, and cost-effective delivery routes.

The paper is organized as follows. Section 2 reviews previous studies on multimodal green mail delivery and the application of the ALNS algorithm. Section 3 formulates the van–walk and van–bike delivery models. Section 4 details the developed ALNS as the primary component of the optimizer system, along with its constituent procedures, including local search and removal and insertion operators. Section 5 outlines the optimizer system created for Australia Post. Section 6 reports the results of the MG-VRP, validates the optimizer in a second real-world case study, and extends the system capability through van–bike delivery. Finally, Section 7 concludes the practical value of the optimizer system developed for Australia Post.

2 Literature review

There are limited studies on multimodal green mail delivery. For this reason, the multimodal vehicle routing problem and the green vehicle routing problem are first elaborated separately in a general context, and the focus then narrows to mail delivery within these fields. Finally, we review the related studies that incorporated ALNS for last-mile delivery optimization.

2.1 Multimodal last mile mail delivery

Multimodal last-mile delivery refers to using more than one mode of transport to complete the final segment of delivery. Unlike traditional single-vehicle routes, multimodal systems coordinate these modes in a complementary way to improve efficiency, flexibility, and service performance. The concept originates from multimodal freight transport in long-haul logistics, where goods move across modes such as rail, trucks, and ships. In those systems, the choice of mode is shaped by factors including transport costs, technological advancements, and organizational practices. This same principle is now being adapted to urban and local neighborhood scales, where it supports more efficient and flexible mail distribution in last-mile delivery (Ferrari, 2014; Sanchez Rodrigues et al., 2015) and improves reliability, reduced manual workload, and enabled scalable, low-impact postal operations (Christensen et al., 2021).

Today, multimodal last-mile delivery has evolved into several distinct variants, leveraging combinations such as vans with walking (Thompson and Zhang, 2018; Wehbi et al., 2022), vans with (cargo) bikes (Nguyen et al., 2023), and even trucks deploying drones (Morim et al., 2024) or sidewalk robots (Kahalimoghadam et al., 2026) in experimental programs. Among these, the van-walk delivery model stands out as a strategy in which a van operates as a mobile depot, stopping at strategically selected locations from which deliveries are completed on foot. This system is particularly advantageous in urban environments characterized by pedestrian zones, narrow streets, or delivery time restrictions, where direct vehicle access to each address is inefficient or infeasible. In its traditional form, the delivery driver parks the vehicle and services a cluster of nearby addresses on foot, during which the van is not actively involved in the delivery process (McLeod et al., 2020). In an alternative configuration, dedicated on-foot couriers, or porters, receive parcels from the van and complete the deliveries independently, allowing the vehicle to proceed to the next zone (Wehbi et al., 2022).

The van-bike delivery model, an increasingly adopted multimodal strategy, enhances last-mile efficiency by combining motorized vans with cargo bikes. In most implementations, vans deliver parcels to intermediate hubs or mobile depots, from which bikes complete the final customer deliveries (Michalakopoulou et al., 2024). While this two-echelon structure improves delivery time and reduces emissions through coordinated mode deployment (Hur and Won, 2023), it may require additional infrastructure such as consolidation centers, which entail considerable upfront investment.

While two-mode last-mile delivery has been examined mainly for e-commerce and parcel contexts with rising demand (Janjevic et al., 2021; Janjevic et al., 2019), mail delivery has received far less attention despite distinct operational constraints and declining volumes. Additionally, these modes are often implemented in separate echelons of the delivery network. For example, van delivery is typically assigned to the first echelon, whereas the second echelon is executed on foot (Martinez-Sykora et al., 2020).

However, there are few studies offer notable contributions. Luo et al. (2023), for instance, introduced a collaborative paired-mail-carrier approach to minimize truck usage, though their analysis remains confined to a single mode without addressing declining mail volumes. On the other hand, Laseinde and Mpofu (2017) tackle the critical issues of diminishing mail demand and high costs in rural areas by integrating parcel and mail distribution into a unified last-mile model. Nevertheless, their broader emphasis on postal efficiency overshadows the complexities unique to mail delivery, particularly regarding multimodal integration. More recently, Álvarez-Miranda et al. (2025) redesigned postal districts into a multimodal network to improve workload balance and geographic compactness. However, their approach did not address multimodal transitions, such as van-to-walk or van-to-bike, which are becoming increasingly important for sustainable postal logistics. Blagojević et al. (2023) likewise proposed flexible last-mile models under universal service obligations, but their approach omitted van-to-walk and van-to-bike combinations that are essential for low-impact, adaptive mail delivery.

2.2 Green last mile mail delivery

The green vehicle routing problem (GVRP) seeks routes optimization while satisfying objective function(s), mainly costs and emission externalities. However, in GVRP the delivery modes are sustainable and environmentally friendly to reduce fossil fuel consumption (Pak and Mun, 2024), in alignment with the Paris Agreement and the United Nations Sustainable Development Goals (Kahalimoghadam et al., 2024). Accordingly, CO₂ emissions calculation are incorporated into the GVRP as an objective or constraint and are estimated based on fuel consumption, traffic congestion, or speed (Xue and Zou, 2025). While estimation methods based on the aforementioned factors appear to be more accurate, they require precise data. This requirement reinforces the need for extensive sensors for data collection, which is usually costly (Ryu et al., 2013). As a result, authors often respond by assuming the required data or by developing a separate model to generate it. Both approaches introduce cumulative errors in the final emission estimation and subsequent analyses (Schulte et al., 2024). Hence, in this study, CO₂ emissions are calculated based on vehicle kilometers traveled, which is widely accepted (Lindsay et al., 2011; Mun and Jung, 2025).

Similar to the multimodal vehicle routing problem, in the GVRP the majority of studies have focused on parcel and goods delivery. In a rare application of the GVRP to mail delivery, Choubassi et al. (2016) evaluated the use of cargo cycles for last-mile mail delivery to minimize total routing cost. However, multimodal delivery and walking as an alternative green delivery mode were not considered. In contrast, Azorli and Meira (2022) and Wang et al. (2021) focused exclusively on mail delivery by walking and did not examine cycling mode. Hence, there is a lack in the current literature in evaluating multimodal green mail delivery to optimize service time while reducing CO₂ emissions.

2.3 ALNS

ALNS was first introduced by Ropke and Pisinger (2006) as an heuristic for the pickup and delivery problem with time windows, wherein various destroy and repair operators are adaptively chosen based on their prior effectiveness. ALNS improves robustness by employing segment-based scoring and periodic weight updates, which stabilize operator evaluation while maintaining sensitivity to recent performance, thereby avoiding the dominance of early successes in guiding the search process (Keskin and Çatay, 2016). The success of ALNS is attributed to its exploitation of diverse and extensive neighborhoods, adaptive selection of operators based on historical efficacy, and integration with simulated annealing, which enable thorough solution space exploration and effective avoidance of local optima. Hence, ALNS has found extensive application across numerous VRP variants owing to its flexibility and robust performance (Windras Mara et al., 2022).

In recent studies, ALNS has been effectively applied to multimodal VRPs, addressing complex coordination between delivery modes, typically trucks and drones, through adaptive or hybrid strategies (Peng et al., 2025; Zhang and Zeng, 2025). Its ability to restructure routes and adapt based on operator performance makes it well-suited for dynamic, multi-objective, and constraint-intensive last-mile delivery environments. For example, Mirzapour Al-e-Hashem et al. (2024) presented a hybrid truck–drone delivery model where an electric truck handled parcel distribution and drone coordination, aiming to minimize delivery completion time. A hybrid metaheuristic combining Clarke & Wright, a genetic algorithm, and ALNS was used to generate and refine truck routes and drone assignments. However, the focus remains on parcel delivery, and drone operations differ significantly from the walking and biking modes examined in our study.

ALNS further supports multimodal delivery systems through adaptive removal, insertion, and splitter heuristics, which dynamically explore large solution spaces to identify cost-efficient and high-performing distribution strategies cost (Keskin and Çatay, 2016; Ropke and Pisinger, 2006). This approach significantly improves delivery efficiency and reduces total transportation. For example, Guo et al. (2022) employed an ALNS framework to integrate three transport modes, truck delivery, joint truck–autonomous vehicle delivery, and autonomous self-delivery to balance cost and efficiency. In a multi-echelon multimodal delivery problem, Wang et al. (2024) utilized ALNS to optimize routing, scheduling, and mode selection within space–time networks.

Focusing on bike-based delivery systems, Fontaine (2022) addressed load- and gradient-dependent travel speed variations in cargo bike routing using an ALNS framework. The study demonstrated notable improvements in travel time and route feasibility over traditional methods. Rave and Fontaine (2025) further expanded its application by integrating ALNS into location and routing decisions for micro-hub-based bike delivery systems, proving that accounting for load-dependent travel times substantially enhances both cost efficiency and adherence to time windows. Similarly, Masmoudi et al. (2021) applied an advanced ALNS framework to heterogeneous fleets involving electric and conventional bikes, highlighting its flexibility in mixed-mode and synchronized operations. Collectively, these studies show that ALNS provides a powerful adaptive mechanism for optimizing delivery routes, reducing operational costs, and supporting sustainable urban logistics.

In the study most closely aligned with our research, Le Colleter et al. (2023) tackled the park-and-loop routing problem for mail delivery, wherein vehicles are required to park at designated feasible locations prior to conducting deliveries on foot. Despite utilizing a tailored version of ALNS for problem-solving, the model imposed a simplifying assumption of customers could only be served by walking once the van is parked, thereby excluding the possibility of direct van-based deliveries.

In this study, ALNS was selected to solve the MG-VRP because the proposed framework includes diverse delivery configurations, namely van-only, van–walk, and van–bike operations. ALNS has shown strong flexibility across a wide range of VRP variants and is therefore well suited to this type of problems. Previous studies on VRP have also indicated that exact methods are typically limited to relatively small instances, while larger and more complex problems are commonly solved using heuristic or metaheuristic algorithms because of computational limitations (Fernández Gil et al., 2022; Garside et al., 2024). As the present study considers instances of up to 4,000 nodes, ALNS is considered an appropriate solution approach.

Although the concept of van–walk and van–bike delivery has been examined in earlier studies, no study has formulated and implemented these strategies in a large-scale, real operational network. Some have explored collaboration between vans and green modes, for example an empirical feasibility study of van–walk parcel delivery with wheeled bags in London that did not offer a mathematical formulation or apply optimization methods (Allen et al., 2021). Another study focuses on van–walk operations, estimating dwell-time regressions from observed parcel tours and developing a tour-simulation that varies maximum walking distance to quantify effects on stop consolidation, total tour time, and vehicle-miles traveled, yet it remains empirical and simulation based without a route-optimization algorithm (Dalla Chiara and Goodchild, 2025). Hence, existing studies still lack large-scale, data-driven multimodal delivery systems capable of handling thousands of customers, numerous potential parking locations, and complex terrain across wide geographic areas.

3 Mathematical model of multimodal VRP

Presented herein is a rigorous formulation of a multimodal urban distribution network design problem that explicitly incorporates the operational characteristics of green delivery modes, such as walking and biking. By integrating these features, the model attains enhanced realism and practical applicability, allowing customers to be served either directly by vans or via green modes initiated from strategically selected parking locations determined by the solution algorithm to efficiently serve customer clusters. The network is structured around several fundamental sets: N C, the set of customer nodes, N D, the set of depots, N P, the set of candidate nodes for parking locations, and N, the set of all nodes excluding depots, that is, N = N C ∪ N P. The set W comprises postmen who deliver on foot, while B denotes postmen who deliver by bike. The set V agg comprises all vans included in the modelling. For modeling purposes, vans used solely for direct delivery and not for green-mode operations are represented by the set V, whereas V ′ denotes vans that either deliver directly to customers or enable delivery by walking or cycling. Accordingly, the aggregate van set is V agg = V ∪ V ′.

The transportation modes vary by delivery configuration. For van-walk delivery, Q = V a g g ∪ W; for van-bike delivery, Q = V a g g ∪ B. For each mode q ∈ Q, the parameter t i j q denotes the travel time between nodes i and j, where i , j ∈ N. c q represents the fixed cost associated with mode q, and c q o specifies mode’s operating cost. Additionally, all vehicles and postmen must return to the depots by time T max.

The model employs several decision variables. x i j q is a binary variable indicating whether transportation mode q ∈ Q travels from node i to node j , where i , j ∈ N. y p v is a binary variable signifying whether vehicle v in V ′ parks at spot p ∈ N P. z i j q is a binary variable that equals 1 when the cluster size of customer i ∈ N C serving by a green mode q (in the van–walk model, q ∈ W; in the van–bike model, q ∈ B) does not exceed the capacity of mode; otherwise, z i j q = 0. l j is a variable representing the total travel time of the green-mode loop that starts and ends at node j ∈ N P. s i is service start time at node i ∈ N while a i q is service time at node i ∈ N when served by mode q ∈ Q . We assume that service time for each mode is identical for all customers and identical for all parking nodes, although customer and parking service times differ.

The model incorporates several costs, capacity, and operational parameters that characterize the behavior of each delivery mode. For each mode q ϵ Q, the fixed cost c q represents the cost incurred when a vehicle or courier of that type is deployed, while the operating cost of modes is captured by c q o. The average travel speed v q is used to compute mode-specific travel times across the network. Fleet availability is represented through k d q, which specifies the number of vehicles of type q ϵ Q that can be dispatched from depot d ∈ N D. The load-carrying capacity of each mode is given by η q, denoting the maximum number of mail items that can be transported by a van, bike, or walking courier. To ensure the economic viability of green-mode operations, c q min imposes a minimum cluster size for walk and bike tours. Additionally, the parameter C a p defines the maximum number of mail items a postman can carry while performing a green-mode loop, ensuring that cluster sizes remain physically feasible. A comprehensive description of all sets, parameters, and variables is provided in Supplementary Appendix A.

The mathematical formulation presented below applies to the van–walk delivery configuration. For the van–bike delivery scenario, all sets, parameters, and variables associated with the walking mode should be replaced with their respective counterparts for the bike mode.

Min ∑ i ϵ N ∪ N D ∑ i ϵ N ∪ N D ∑ q ϵ Q t i j q x i j q ⏟ t r a v e l t i m e + ∑ i ϵ N D ∪ V ∑ j ϵ N C ∪ N P ∑ q ϵ Q c q x i j q ⏟ f i x e d c o s t + ∑ i ϵ N ∪ N D ∑ j ϵ N ∪ N D ∑ q ϵ Q c q o d i j q x i j q ⏟ o p e r a t i o n a l c o s t (1)

∑ d , o ∈ N D x d o q = 0 , ∀ q ϵ Q (2)

∑ i ∈ N , i ≠ j , q ∈ Q x i j q = 1 , ∀ j ϵ N C (3)

∑ v ∈ V a g g y p v ≤ 1 , ∀ p ϵ N W (4)

∑ i ∈ N , q = v a g g x i j q ≤ M × y p v , ∀ j ϵ N W ∪ N B , p ϵ N P (5)

∑ j ∈ N c x p j r = y p v , ∀ p ϵ N P , r ϵ N W ∪ N B , v ϵ V ′ (6)

∑ i , j ϵ N x i j q = ∑ j , k ϵ N x j k q , ∀ q ϵ Q \ V ′ (7)

∑ i , j ϵ N x i j q = ∑ j , k ϵ N x j k q , ∀ q ϵ V ′ (8)

∑ d ϵ N D , j ϵ N x d j q = ∑ j ϵ N , o ϵ N D x j o q , ∀ q ϵ V ′ (9)

∑ j ϵ N x d j q ≤ k d q , ∀ q ϵ Q , d ϵ N D (10)

s i + t i j q + a i q − s j ≤ M 1 − x i j q , ∀ q ϵ Q , i ϵ N , j ϵ N (11)

s i + t i j q + a i q − l j ≤ M 1 − x i j q , ∀ q ϵ V ′ , i ϵ N C , j ϵ N P (12)

s i + t i j q + a i q ≤ M 1 − x i j q + T max , ∀ q ϵ Q \ V ′ , i ϵ N , j ϵ N D (13)

y v p ϵ 0 , 1 , ∀ p ϵ N W , ∀ v ϵ V ′ (14)

x i j q ϵ 0 , 1 , ∀ i , j ϵ N , ∀ q ϵ Q (15)

z i j q ϵ 0 , 1 , ∀ i , j ϵ N R , ∀ q ϵ V ′ (16)

s i ϵ 0 , T max , ∀ i ϵ N (17)

∑ i ϵ N ∑ j ϵ C q x i j q ≤ η q , ∀ q ϵ V ∪ W (18)

∑ j ϵ C q x i j q ≤ z i j q . C a p , ∀ i ϵ N P , q ϵ V ′ (19)

∑ j ∈ C q x i j q ≥ c q min · z i j q , ∀ i ∈ N , q ∈ V ′ (20)

E C O 2 = ∑ q ∈ Q ∑ i ∈ N ∪ N D ∑ j ∈ N ∪ N D e q d i j q x i j q (21)

The objective function in Equation 1 minimizes the total cost of the multimodal delivery system by incorporating travel time, mode-activation costs, and arc-traversal costs. Constraint (Equation 2) prohibits direct travel between depots, ensuring vehicles cannot move from one depot to another without serving customers. Constraint (Equation 3) guarantees each customer is visited exactly once by a single vehicle. Constraint (Equation 4) restricts each parking to a maximum of one aggregate vehicle assignment. Constraint (Equation 5) limits the number of aggregate vehicles assigned to any parking within system-defined bounds. Constraint (Equation 6) activates parking only when it is used for customer service. Constraints (Equations 7, 8) ensure flow conservation at each node for all vehicles and routes. Constraint (Equation 9) requires all vehicles to return to a depot, forming closed routes. Constraint (Equation 10) restricts the number of departing vehicles from each depot to the available fleet size.

Constraint (Equation 11) maintains time consistency by enforcing arrival at each node only after the sum of travel and service times, conditional on route activation. Here, M is a sufficiently large constant used to deactivate a time constraint when an arc is not part of the selected route, ensuring that timing relationships apply only to arcs that are actually used. Constraint (Equation 12) enforces the latest allowable arrival or service time at parking or cluster nodes for specified transportation modes. Constraint (Equation 13) limits the total operational time of each vehicle to the maximum allowed working hours. Equations 14–17 define decision variables, with binary variables representing assignment and routing decisions, and continuous variables indicating service start times. Constraint (Equation 18) ensures that the demand served by each van or walk does not exceed its capacity. Constraint (Equation 19) restricts the number of customers served in a cluster by walking to the capacity limit when the stop point is a parking. Constraint (Equation 20) enforces a minimum cluster size by requiring that, if a green mode is assigned to a cluster, the number of customers served must be at least ​c q min, ensuring operational efficiency for walking or biking deliveries. In the van–bike delivery model, Constraints (Equations 18–20) are adjusted to account for the unique operational characteristics and capacity limitations of the bike mode. Finally, emissions based on vehicle kilometer traveled and transportation mode are calculated in Equation 21. In this equation, e q denotes the tailpipe emission factor (kg CO₂/km), which is zero for green modes, while it is equal to 0.19 kg CO₂/km for vans.

4 Solution methodology

The solution methodology begins with data preparation, then constructs an initial set of van-only routes and finally refines these routes and extends them to cooperative delivery. Section 4.1 introduces the data preparation framework that underpins the optimizer, while Section 5 describes how this step is implemented in the prototype system. Section 4.2 builds an initial van-only solution using a parallel Clarke–Wright savings (CWs) heuristic, followed by 2-opt (Helsgaun, 2000) and 3-opt (Li et al., 2016) local search to shorten tours and remove crossings. Section 4.3 develops an ALNS tailored to large delivery instances, seeded with the Clarke–Wright solution. The ALNS iterates destroy–repair cycles that combine multiple removal and insertion operators with post-repair local search based on 2-opt and Lin–Kernighan moves and adapts operator weights online using a simulated-annealing acceptance scheme to balance diversification and intensification. Sections 4.3.1–4.3.6 introduce cooperative delivery in three steps. First, density-based spatial clustering of applications with noise (DBSCAN) (Li et al., 2016) is employed to identify walkable customer clusters. Next, the algorithm constructs and refines closed walking loops using Clarke–Wright and 2-opt. Finally, a cost-based rule determines whether each customer is served by a van or a green mode, considering travel speeds, service times, distances, and the time required for mode switching. The same framework is subsequently applied to van–bike delivery. An overview of the proposed methodology is presented in Figure 3.

4.1 Data preparation

Prior to the optimization process, several preparatory steps are required to compute input parameters and define feasibility conditions for both arc-level routing and cluster-level assignments. These steps ensure the model’s computational tractability and adherence to operational constraints. For each pair of nodes i , j ∈ N and each transportation mode q ∈ Q, the travel distance d i j q ​is calculated using the Haversine formula, which accounts for the Earth’s curvature by estimating the great-circle distance between two geographic coordinates. The formula relies on two intermediate values: the angular distance a, Equation 22, which incorporates the latitude and longitude differences in radians, and the central angle c, Equation 23, which is used to compute the final distance, d i j q, in Equation 24 based on Earth’s radius. where Δ φ = φ j − φ i​ and Δ λ = λ j − λ i represent the differences in latitude and longitude (in radians), φ and λ denote the latitude and longitude of nodes i and j, respectively, and R is the Earth’s radius (typically 6,371 km). Once the distance is computed, the travel time for each mode q is derived as d i j q = t i j q * v q, where v q is the average speed for mode q (e.g., 30 km/h for vans, 5 km/h for walking).

a = s i n 2 Δ λ 2 + cos φ i · ⁡ cos φ j · s i n 2 Δ λ 2 (22)

c = 2 · arctan ⁡ 2 a , 1 − a (23)

d i j q = R · c , ∀ i , j ∈ N ∪ N D , q ∈ Q (24)

To ensure computational efficiency and practical applicability, two feasibility matrices are introduced to restrict the optimization model to valid routing and clustering configurations. These matrices incorporate spatial and modal constraints, thereby reducing the solution space and eliminating infeasible decisions. The arc feasibility matrix defines whether a given arc i , j can be traversed by mode q ∈ Q. The binary parameter δ i j q​ equals 1 if arc i , j is accessible by mode q, based on physical, regulatory, or infrastructure conditions (e.g., the absence of walkways between certain nodes), and 0 otherwise, as defined in Equation 25. To ensure vehicles and postmen only traverse permitted arcs, the routing variable x i j q, is constrained by δ i j q in Equation 26.

δ i j q = 1 , if a r c i , j is usable by mode q 0 , otherwise , ∀ i , j ∈ N ∪ N D , q ∈ Q (25)

x i j q ≤ δ i j q , ∀ i , j ∈ N , q ∈ Q (26)

In addition, a cluster-feasibility matrix is constructed during data preparation to identify which customers can be served from each parking location under each green mode. The binary parameter θ p j q equals 1 if customer j ∈ N C lies within the maximum allowable walking or biking time from parking location p ∈ N P (i.e., t p j q ≤ T max) and if the resulting cluster size satisfies the mode-specific minimum and maximum limits; otherwise, θ p j q = 0. This precomputed parameter is used exclusively in the clustering and initial-solution construction steps (Sections 4.1–4.3) to filter infeasible parking–customer–mode combinations before optimization, ensuring that green-mode segments are formed only in spatially and operationally viable contexts.

4.2 Initial solution construction

To generate the initial solution, the algorithm employs the parallel CWs heuristic. This classical approach begins by assuming that each customer j ∈ N C is served by an individual direct route from the depot and back, resulting in an initial total cost of ∑ j ∈ N C c 0 j + c j 0, where c 0 j and c j 0 denote the travel costs from the depot to customer j and the return trip. The parameter S i j denotes the savings obtained by connecting customers i and j directly rather than serving them through separate depot trips, where i , j ∈ N C. The heuristic then computes the savings for merging customer pairs i , j ∈ N C into a single route, as defined in Equation 27.where S i j​ denotes the savings obtained by connecting customers i and j directly rather than visiting them via separate depot trips, and c i j​ is the travel cost between them. The savings are ranked in descending order, and feasible customer pairs are iteratively merged based on the highest available savings, subject to route feasibility and capacity constraints. This process continues until no further beneficial merges are possible.

S i j = c 0 i + c j 0 − c i j (27)

The initial solution generated by the CWs heuristic provides a rapid means to obtain feasible delivery routes and enables the identification of customer clusters suitable for assignment to alternative transport modes, such as walking and biking. However, due to its greedy and sequential nature, the CWs approach often produces routes that are locally inefficient, with superfluous detours, route crossings, or suboptimal visit sequences. To overcome these local inefficiencies, local search algorithms, specifically 2-opt and 3-opt, are employed to systematically improve each route. The purpose of incorporating these methods is to refine the initial solution by removing redundant connections and optimizing the order of customer visits, which is particularly crucial in large-scale, multimodal urban distribution settings where local errors can compound and adversely affect system-wide performance.

The mathematical formulation of the 2-opt algorithm is as follows. Given a route r = v 0 , v 1 , … , v n , v 0, where v 0​ denotes the depot and v 1 , … , v n represent the sequence of customer nodes, the 2-opt procedure examines all possible pairs of non-adjacent edges within the route. For each valid pair i , j, with 0 < i < j < n, let the current edge be v i − 1 , v i and v j , v j + 1, where the indices are selected such that the two edges do not overlap. Subsequently, the algorithm constructs a new route by reversing the sequence of nodes between v i​ and v j​, yielding the following modified route as v 0 , … , v i − 1 , v j , v j − 1 , … , v i , v j + 1 , … , v n , v 0. This reordering effectively eliminates any potential crossing between the selected edges, potentially reducing the total route length. The incremental change in route cost associated with this operation is computed as Equation 28.where c a , b denotes the travel cost or distance between nodes a and b. A negative value of Δ indicates an improvement in the objective function. In such cases, the algorithm replaces the original sequence with the improved one. This process is repeated exhaustively for all eligible pairs i , j until no further improvement can be obtained, at which point the route is considered 2-optimal.

Δ = c v i − 1 , v j + c v i , v j + 1 − c v i − 1 , v i + c v j , v j + 1 (28)

The 3-opt algorithm generalizes this approach by simultaneously considering three distinct edges within the route. For a given route r = v 0 , v 1 , … , v n , v 0, the algorithm identifies three breakpoints, corresponding to indices 0 ≤ a < b < c < n, and removes the edges v a , v a + 1, v b , v b + 1 and v c , v c + 1. The removal of these edges partitions the route into four segments, presented in Equation 29.

A : v 0 , … , v a B : v a + 1 , … , v b C : v b + 1 , … , v c D : v c + 1 , … , v n , v 0 (29)

The 3-opt procedure systematically explores all non-trivial reconnections of these segments, considering both their original and reversed orientations. Each reconnection generates a new candidate route, and the cost associated with each configuration is given by the sum of the new connecting edges. For a candidate reconnection, the change in route cost is evaluated by Equation 30.where c n e w k denotes the cost of the k − t h newly created edge in the candidate reconnection. If any reconnection yields Δ < 0, the route is updated accordingly. This process is applied iteratively across all valid edge triplets and reconnection patterns until no further cost reductions are attainable, at which point the route is deemed 3-optimal.

Δ = c n e w 1 + c n e w 2 + c n e w 3 − c v a , v a + 1 + c v b , v b + 1 + c v c , v c + 1 , (30)

4.3 Adaptive large neighborhood search (ALNS)

Following the local improvement phase, the solution is further refined using the ALNS metaheuristic. ALNS is well suited to complex, highly constrained VRP variants such as multimodal mail delivery, where routes combine van-walk and van-bike segments. The method explores the solution space through iterative destroy and repair steps that balance search intensification and diversification. In the destroy phase, elements of the current solution are selectively removed, including stops, sub routes, or transfer chains involving specific modes. The repair phase reinserts the removed elements with heuristics that account for modal constraints, including transfer locations, access restrictions, and vehicle compatibility. A central feature is adaptive operator selection, which dynamically prioritizes removal and insertion strategies based on observed performance.

4.3.1 Removal operators

To promote diversification and avoid local optima, the ALNS framework integrates four removal heuristics: Random Removal (RR), Cluster Removal (CR), Closest Nodes Removal (CNR), and Peripheral Nodes Removal (PNR). Let R be the set of routes and N r​ the customer nodes in route r. Let N C = ⋃ r ∈ R N r be the set of all served customers. Given a destruction rate d ∈ ( 0 , 1 ] , define κ = ⌊ d · ∣ N C ∣ ⌋ as the number of nodes to remove in a destroy step.

Random Removal (RR): RR selects a set of customers for deletion based on a specified degree of destruction (Fontaine, 2022). By removing randomly chosen nodes from the current routes, it creates space for repair operators and injects stochasticity that helps diversify the search. In this research, the RR heuristic selects a subset of κ customers uniformly at random from the complete set of served customer nodes N C​. Let R be the current set of routes and N r​ the set of nodes in route r ∈ R. The set of nodes selected for removal is denoted as N R R and is defined in Equation 31.

N R R ⊆ ⋃ r ∈ R N r w h e r e ∣ N R R ∣ = κ (31)

This method introduces stochastic perturbation to escape local minimum and promote exploration of new regions in the solution space.

Cluster Removal (CR): The CR heuristic targets geographically concentrated customer nodes identified using the DBSCAN clustering algorithm. Is a widely used unsupervised machine learning algorithm designed to identify clusters of arbitrary shape in spatial datasets. It operates by grouping together points that are closely packed; while marking as outliers those points that lie in low-density regions ccc. Let ϵ denote the maximum distance threshold and let m i n P t s represent the minimum number of nodes required to form a cluster. The set of nodes selected for removal, denoted N c l ​, consists of at least κ customer nodes identified within one or more clusters, as defined in Equation 32. This approach facilitates the removal of densely located customers, thereby enabling re-clustering and the construction of more efficient sub-tours.

N C R = D B S C A N N C ; ϵ , m i n P t s w h e r e ∣ N C R ∣ ≥ κ (32)

This method allows the removal of densely located customers, thus opening space for re-clustering and more efficient sub-tour formation.

Closest Nodes Removal (CNR): CNR identifies pairs of customers that are nearest to one another across different routes, then removes them up to the prescribed destruction level (Mirzapour Al-e-Hashem et al., 2024). By disrupting tightly coupled cross-route pairs, this operator can unlock more efficient inter-route exchanges. Here, CNR identifies κ nodes that are closest to others on different routes. Let d i j q​ be the travel distance between nodes i and j under mode q. The algorithm selects node pairs i , j, where i ∈ N r 1, j ∈ N r 2 and r 1 ≠ r 2​, such that the cumulative pairwise distance is minimized. The resulting set of nodes selected for removal is denoted by N C N R ​, as defined in Equation 33. This approach aims to disrupt strongly interconnected inter-route structures, potentially enabling more effective route reconfiguration during the reconstruction phase.

N C N R ⊆ ⋃ i ∈ N r 1 , j ∈ N r 2 , r 1 ≠ r 2 i , j where ∣ N C N R ∣ = κ (33)

Peripheral Nodes Removal (PNR): The PNR heuristic aims to eliminate customer nodes located at the peripheries of existing routes, specifically those positioned at the beginning or end of each route (Zhang and Zeng, 2025). For each route r ∈ R , let i r s​ and ​i r e represent the first and last customer nodes visited, excluding the depot. The algorithm constructs the removal set N P N P​ by iteratively selecting nodes from these endpoints with the shortest distances to the depot, until the desired destruction level κ is achieved. Equation 34 represents this process.

N P R P ⊆ ⋃ r ∈ R i r s , i r e where ∣ N P R P ∣ = κ (34)

4.3.2 Insertion operators

The insertion phase reconstructs complete routes by reinserting removed customers while respecting feasibility limits such as maximum travel time.

Greedy insertion: Greedy insertion evaluates all admissible positions for a node and selects the position that yields the smallest increase in route cost, thereby rebuilding solutions efficiently (Ropke and Pisinger, 2006). The incremental cost of placing customer i between consecutive nodes a and b in a route, using delivery mode q, is calculated by Equation 35

δ i = c a , i , q + c i , b , q − c a , b , q (35)

Regret-based insertion: Regret-based insertion prioritizes customers whose best positions are fragile. For each unrouted customer i, compute the two lowest feasible greedy insertion costs over all admissible positions and modes (Masmoudi et al., 2021). The regret value is defined in Equation 36. At each iteration, rank customers by regret in descending order, insert the top customer at its current best position, and then recompute costs and regrets to reflect the updated routes, enforcing capacity and time constraints. Tie-breaks can use the smaller best-cost or the shorter detour.

ρ i = δ i , 2 − δ i , 1 (36)

K-regret insertion: The k-regret insertion heuristic extends the regret-based strategy by evaluating the quality of multiple alternative insertion positions, rather than only comparing the best and second-best options (Windras Mara et al., 2022). For each customer i, the algorithm computes the k lowest feasible greedy insertion costs, denoted δ i , 1 , δ i , 2 , … , δ i , k, where δ i , 1 is the minimum among them. The regret value ρ i is then calculated as Equation 37.

ρ i = δ i , 2 + δ i , 3 + … + δ i , k − δ i , 1 (37)

This metric quantifies the expected loss incurred if the current best insertion position becomes unavailable in future iterations. By prioritizing customers with higher regret values, the heuristic aims to reduce the risk of unfavorable late insertions, which can significantly degrade overall solution quality. The k-regret method thus provides a more informed and globally aware insertion strategy that balances immediate cost minimization with long-term flexibility in route construction.

4.3.3 Adaptive selection of removals and insertions

In each iteration of the ALNS framework, one removal and one insertion operator are selected using an adaptive probability mechanism. The selection probability of each operator is proportional to a weight that reflects its historical performance during recent iterations. These weights are periodically updated based on the quality of solutions generated by the operators. Specifically, scores are awarded when an operator contributes to a new global best solution, improves the current solution, or results in an accepted but non-improving solution. The accumulated score is normalized by the number of times the operator was applied, and the weights are adjusted using a reaction factor to balance past and recent performance.

This adaptive selection process offers several advantages. First, it enables the search to dynamically favor operators that are well-suited to the current search landscape, without requiring manual tuning of operator preferences. Second, by maintaining a diverse pool of operators and probabilistically selecting among them, the algorithm preserves a balance between intensification and diversification. This is especially beneficial in multimodal delivery problems, where different removal and insertion strategies may be more effective in different phases of the search or under varying operational constraints (e.g., van-to-walk vs. bike-to-cluster transitions). Importantly, ALNS adapts operator usage over time: heuristic weights are refreshed every 100 iterations over a total iteration horizon of 3,000 iterations, enabling the framework to emphasize high-performing operators while still allowing exploration of less frequently used ones. This dynamic weighting mechanism ensures that the search maintains both exploitation of effective heuristics and diversity in the solution space, thereby enhancing robustness and convergence quality.

To formalize this process, let t denote the current iteration. The removal operator R k ∈ R is selected with probability p R k t, computed as shown in Equation 38. This probability is proportional to the operator’s current weight w R k t, normalized over all n r removal operators:

p R k t = w R k t ∑ j = 1 n r w R j t . (38)

Once selected, the removal operator R k is applied to the current solution S, removing a subset N rem ⊆ N C of size κ, resulting in the partial solution S ′, as given in Equation 39:

S ′ = R k S , κ (39)

Similarly, the insertion operator I l ∈ I is chosen with probability p I l t, defined analogously in Equation 40:

p I l t = w I l t ∑ j = 1 n i w I j t . (40)

The selected operator I l then reconstructs a complete solution by reinserting the removed nodes from N rem into S ′, producing a new candidate solution S new, as shown in Equation 41:

S n e w = I l S ′ (41)

This mechanism ensures that operator selection is guided by recent empirical performance while maintaining flexibility and diversity in the search.

4.3.4 Local search algorithms

To further refine solutions, local search is selectively applied after certain destroy–repair cycles, particularly when a new best solution is found. Given the potential computational cost of local search, it is not applied in every iteration but rather used strategically to intensify promising solutions. The first local improvement strategy involves applying the 2-opt heuristic (Croes, 1958), which removes intersecting edges and reconnects segments to reduce total travel cost. Due to its simplicity and speed, 2-opt offers effective gains with minimal overhead.

For deeper local optimization, a variant of the Lin–Kernighan heuristic is employed. This method performs a sequence of edge exchanges within the route, enabling exploration of significantly larger neighborhoods compared to 2-opt. We adopt an efficient implementation similar to the one proposed by (Helsgaun, 2000), which incorporates adaptive move selection and sensitivity analysis to guide the search. Together, these local search components provide strong intensification without overwhelming the overall computational time and are especially beneficial in complex multimodal instances where route geometry significantly impacts feasibility and efficiency.

4.3.5 Cooperative delivery

To incorporate walking as a viable alternative to vehicle service, we first extract stabilized vehicle routes and apply the DBSCAN clustering algorithm to identify spatially dense groups of customer nodes. These clusters represent candidate areas for foot delivery, where proximity supports efficient walking tours. Each cluster is then used to generate a closed walking loop, initially constructed with a CWs heuristic and subsequently refined using 2-opt local search. A mode selection rule is applied to choose between vehicle or walking service for each cluster, based on total cost incorporating travel speed, service time, distance, and transfer penalties. This approach facilitates cooperative delivery while maintaining operational consistency and optimizing resource allocation.

4.3.6 Acceptance criterion

After each destruction and repair step, the newly generated solution S n e w​ is compared with the current solution S. If the objective function value of S n e w​ is lower than that of S, the new solution is always accepted. Otherwise, S n e w​ is accepted with a probability determined by the simulated annealing criterion. This can be mathematically represented by Equation 42.

S = S n e w if f S n e w < f S o r U < ⁡ exp − f S n e w − f S T S otherwise (42)

Where f . denotes the objective function, T is the temperature parameter, and U ∼ U 0 , 1 is a random variable sampled from the uniform distribution on the interval [0, 1]. If the newly accepted solution improves upon the incumbent best solution S best, the best solution is updated accordingly, as shown in Equation 43:

S best = S , if f S < f S best S best , otherwise (43)

5 Last-mile multimodal delivery planning and optimization system

This section provides a step-by-step guide from data and software preparation to run the developed Last-Mile Multimodal Delivery Planning and Optimization system. Firstly, the data preparation process is explained. This is followed by a detailed explanation of the software settings. Finally, the designed user interface and how it integrates all components of the optimizer system will be explained. These steps, illustrated in Figure 4, are described in detail in the remainder of this section.

The primary inputs for the model are the demand points, provided by Australia Post, are in shapefile format, which is compatible with ArcGIS Pro software. The data typically contains approximately 50 attributes, offering a detailed overview of the case study. For instance, they indicate how many delivery tours (or rounds) are required to service all customers in the case study before optimization. Longitude and latitude of the demand points are also essential attributes used in the optimization process.

The origin–destination matrix is generated from the demand-point coordinates using OpenStreetMap (OSM). However, several preprocessing steps are required before the matrix can be built. To facilitate this, a Python model, Origin–Destination Matrix Builder, was developed to call routing services based on the OSM API and compute travel distances. Because Python interacts more effectively with Excel files than with shapefiles in this workflow, the selected demand data must first be converted into Excel format. ArcGIS Pro is used for this conversion and for defining the case study area.

The main optimizer, an ALNS algorithm, was implemented in C++20 and compiled using Microsoft Visual Studio 2022. However, the results, which are the multimodal routes visualizing van segments alongside walk or bike segments, are displayed in ArcGIS Pro. Hence, a Python interface is required to execute the optimizer and pass the results to ArcGIS Pro for visualization. Therefore, Microsoft Visual Studio 2022 is used to support both C++ development and the creation of Python bindings, acting as the bridge between the C++ optimizer and ArcGIS Pro.

Although ArcGIS Pro is compatible with Python, it is recommended to use the Python version embedded in ArcGIS Pro (arcgispro-py3) when developing integrations. Therefore, the ArcGIS Pro Python environment is added to Visual Studio to ensure compatibility and avoid version conflicts. To execute the C++ optimizer model from Python, Pybind11 is used. Since ArcGIS Pro installs its own Pybind11 package inside the arcgispro-py3 environment, it is critical to compile the C++ extension using the exact Pybind11 version and headers provided by ArcGIS Pro to ensure compatibility. These components, C++ optimization engine, Python interface, and ArcGIS Pro visualization, are integrated into a single standalone package, forming a complete decision-support system that has been delivered to Australia Post.

6 Results and applications of multimodal mail delivery

Computational experiments were conducted to assess the performance of the proposed algorithm. The model was executed on a desktop computer equipped with an Intel Core i7 processor operating at 2.8 GHz and 16.0 GB of RAM.

6.1 Van-walk delivery

The Belgrave distribution area in VIC is considered as the first case study in this research, while in the case study there are approximately 15,700 delivery addresses, we focus the analysis on the Belgrave service area covering postcodes 3,160 and 3,158, where 6,524 green delivery points were identified and modeled as network nodes. To assess scalability and robustness, a spectrum of instance sizes was constructed from the Belgrave network, ranging from small subsets to large, near-full cases. The cooperative van-plus-walk methodology was solved on each instance, enabling comparison of routing structure, travel effort, and computational behavior across scales within the same geography.

Table 1 reports the results of applying the developed ALNS to the Belgrave case study for different instance sizes. As can be seen, when the instance size increases, the optimization system performs better in terms of CO₂ reduction and delivery time savings. This improvement occurs because, as the number of demand points increases, the potential clusters that can be served by walking also increase. Since walking, as a green delivery mode, is not only emission free but also faster, the delivery process is accelerated as the total number of stops, the need to find parking, and the time required to get out of and back into the vehicle are reduced.

TABLE 1. Performance of van-walk multimodal delivery across scales.

Figures 5, 6 illustrate the delivery routes for the 200- and 300-customer instances, where bold lines indicate walking segments and thin lines represent van segments. The 200-customer instance exhibits fewer and shorter walking segments compared with the 300-customer instance. In addition, walking segments reduce local zigzagging movements by vans, which explains the delivery time savings reported in Table 1 when the MG-VRP is employed, even when the total distance is similar or slightly higher.

Figure 7 summarizes the weight evolution of destroy operators for the 3,000-customer instance. The destroy operator weights are updated after each cycle based on the improvement achieved. At the beginning, weight fluctuations are more pronounced, while as the iterations progress, the weights move toward a steady-state phase. This behavior indicates that, in the early stages, the solution methodology explores the search space more broadly to identify high-quality solutions. At the end of the iteration limit, the closest node destroy operator attains the highest weight, indicating that it outperforms the other destroy operators. Similarly, Figures 8, 9 present the corresponding trajectories for the repair and local-search phases, where the regret repair and 2-opt operators outperform the other operators, respectively.

6.2 Reflecting customers density

One of the significant challenges faced by Australia Post is the declining demand for mail delivery, primarily driven by digitalization. Since the ratio of green customers to total nodes in the network plays a critical role in last-mile delivery (LMD) planning and route optimization, a fine-tuned function has been integrated into the main model to account for customer density. To ensure this density adjustment is accurately reflected in the input data, extensive effort was dedicated to designing and developing this feature in a way that minimizes any impact on other components of the model. Another essential feature incorporated into the model is the random selection of nodes as customer locations. This addition reflects the real-world variability of mail and small parcel demand. The function has been designed with high flexibility, enabling easy modification of customer density to accommodate diverse scenarios.

The developed model demonstrates substantial flexibility, allowing instant adjustment of customer density parameters, yet this parameter must be calibrated to reflect contemporary mail delivery practices observed by Australia Post. To achieve this alignment, an extensive sensitivity analysis was undertaken to identify the density rate that most accurately represents real operational conditions. Four instance sizes were examined, small, medium, large, and extra-large, corresponding to 300, 600, 1,200, and 2,400 nodes, respectively. The density parameter was initially set at 100% and then progressively reduced to 50, 35, and 25%. However, due to the limited number and low spatial concentration of customers in the 300 and 600-node instances under the 25% scenario, these cases were excluded from further analysis. In total, the study encompassed 14 experimental settings, enabling a comprehensive assessment of the model’s responsiveness and robustness across a wide range of density conditions.

These scenarios were compared by measuring key metrics, including the number of routes, the number of vans and walking segments, the total service time, and the computational time. The results of the analysis are summarized in Table 2. The results confirm that higher density and larger instances correspond to more customers being served by walking and, consequently, greater CO₂ emission reductions. Moreover, Figures 10–12 report the number of walk segments for a 1,200-node instance at 35%, 50%, and 100% density, respectively.

TABLE 2. Impact of varying customer density and instance size on routing design and optimizer performance.

Two field validations were conducted in Belgrave to assess the van-walk model against observed operations. In the first, a pre-defined route of 127 customers drawn from 1,200 demand points yielded a predicted total time of 3 h, while the observed duration was 3 h and 20 min; the gap was attributed largely to limited local network familiarity, indicating that the model’s predictions closely reflect real conditions. In the second, routes covering 270 customers out of 1,080 demand points, with no walk segments, produced a 6-h prediction versus 6 h and 30 min observed, comprising roughly 2 h driving and 4 h and 30 min service time over 80 km at an average speed of 40 km/h. Adjusting the assumed van speed from 45 km/h to 40 km/h increases the predicted time by 14 min, bringing the estimate into closer alignment with observations. The 1-min per-mail service assumption implies about 4.5 h of service, consistent with measured times. Relevant route and walk-segment visuals are provided in Supplementary Appendix B.

6.3 Expanding the multimodal delivery model applications

This section develops a van–bike delivery model that addresses limitations of the van–walk methodology by leveraging higher travel speeds to improve practicality in lower-density areas. Methodological refinements increase flexibility and enable seamless transfer to additional case studies with geographic characteristics comparable to Belgrave, replicating the workflow from data preparation to advanced analysis. These enhancements optimize current performance and demonstrate external validity, strengthening robustness and applicability across diverse real-world contexts.

6.3.1 Van-bike implementation

Key parameters distinguishing walking delivery include service time, van–walk transitions, minimum and maximum customers per segment, and walking speed. These were adapted for bike delivery through literature review, model testing, and consultation with Australia Post. Segment size was derived from carrying capacity and delivery mix. With a 39 kg basket limit, mail weights of 250 g and 500 g, and an 80–20 split of small to large mails alongside parcels, the maximum customers per bike segment was set to 40, and a minimum of 4 was imposed to maintain operational viability.

Bike speed was selected via sensitivity analysis. A 25% customer ratio sampled 600 of 1,200 demand points, the subset was fixed, and a new origin–destination matrix was computed to remove location effects. To address heuristic variability, three runs were executed for each speed, and key indicators, including bike segment count, customers served by bikes, service time, and time saved versus van-only delivery, were summarized in Table 3 with averages. The results support an optimal bike speed of 11 km/h. Table 4 compares the final parameter sets for van–walk and van–bike configurations.

TABLE 3. Evaluation of van-bike delivery model performance against various bike speeds.

TABLE 4. Identified parameters for van-bike delivery and comparison with van-walk parameters.

6.3.2 Validating the optimizer system

Katoomba, NSW, illustrated in Figure 13, was selected as the comparative case study due to its geographic and demographic similarity to Belgrave, enabling validation of the developed optimizer system. Australia Post supplied a shapefile with more than 20,000 demand points, the full area was too large for a single analysis, so demand points were stratified by postcodes and delivery routes. Following consultation with an Australia Post expert, a representative subarea was defined that specifies the number of demand points, delivery rounds, and relevant postcodes.

A preliminary analysis was performed to assess the feasibility of van-walk and van-bike delivery for the full study area, using a 25% customer rate. The case study contains 920 customers, and the results are summarised in Table 5. As indicated, van-bike delivery generates more green-mode segments and yields greater absolute and percentage reductions in service time and CO₂ emissions compared with van-walk delivery.

TABLE 5. Van-walk and Van-bike delivery results in Katoomba.

Similar to the previous section, to eliminate the impact of random customer selection and heuristic solutions, 25% of the nodes were randomly selected, and an origin-destination matrix was developed for them. The model was then run three times for both van-walk and van-bike delivery methods. The results are presented in Table 6. In this experiment also van-bike delivery supper passed van-walk delivery in both saving time and reducing CO₂. The locations of green delivery modes, including walking and bike segments, are depicted in Figures 14, 15.

TABLE 6. Performance of the multimodal delivery is Katoomba case study.

7 Conclusion

An adaptive optimizer system was developed and delivered to Australia Post to help them in addressing the inefficiency challenge in last-mile mail delivery due to declining demand. This decision-support system designs the mail delivery procedures, including optimum routes and park locations, which is formulated as an MG-VRP. Mail delivery starts from a depot by vans which serve some customers, while the optimizer system intelligently obtains the optimum clusters to be served by green modes. This capability was supported density-based analysis, combined with randomized customer selection, adaptive operator weighting within the ALNS. These features, alongside transferring optimized solutions to ArcGIS Pro for visualization and operational interpretation, transform the optimization system into a practical expert system, enabling decision makers to design mail delivery procedures. The structure of the designed optimizer system is discussed in Section 5.

Evidence from the Belgrave case study confirms that the MG-VRP yields measurable efficiency gain as instance size increases. In the Belgrave case study, the maximum reduction in delivery service time was 21.67% for the 4,000-customer instance, and the maximum reduction in CO₂ emissions was 11.36% for the 2,000-customer instance, as reported in Table 1. Sensitivity analyses show that performance depends heavily on customer density, which was explicitly modeled and evaluated across 12 scenarios ranging from 300 to 2,400 demand points as presented in Table 2. At lower densities, walking segments declined or disappeared, and multimodal benefits diminished. In contrast, at higher densities, more walking segments were activated, leading to greater service time and emission savings, as illustrated in Table 2 and Figures 10–12.

External validity was explored through the generalization of van–bike operations and by testing the optimization system in Katoomba, in the state of NSW. In Katoomba, van–bike consistently enabled more green delivery segments and consequently lower service time and CO₂ emissions than van–walk. In Katoomba, the van–bike delivery mode enabled more green delivery segments and consequently lower service time and CO₂ emissions than the van–walk. As reported in Table 5, the van–bike configuration outperformed the van–walk configuration, achieving a greater reduction in service time, 8.5% compared with 5.1%, and a greater reduction in CO₂ emissions, 14.46% compared with 7.86%. This finding supports the use of bicycles in areas less suited to walking due to lower density or challenging geography. From a managerial perspective, van–walk is effective in dense neighborhoods, as it reduces service time. In more spread-out areas, van–bike is preferable, mainly due to the higher speed of bicycles compared to walking. The integration of the density function, mode-capacity rules, and operator-weight adaptation within the ALNS provides planners with a practical toolkit to select delivery modes, define segment sizes, and choose parking or transfer points within operational limits.

Beyond academic contribution, this research culminated in an expert decision support system that has been delivered to Australia Post for operational use. The system is now integrated into their daily decision-making framework, enabling planners to evaluate routing strategies under various demand densities and modal configurations. The optimizer’s capability to generate GIS-based route visualizations, coupled with its flexibility to incorporate scenario-specific parameters, makes it a practical asset for improving last-mile efficiency. The adoption of this optimization system by Australia Post marks a successful translation of intelligent research into real-world expert system deployment.

Despite promising results, the study has some limitations. Parameters for speeds, service times, and transitions were fixed, and congestion or curbside interference was not modeled. Average speeds were assumed, although both case studies involve heterogeneous road types with different speed limits. Parking feasibility was pre-screened rather than dynamically optimized. Labor rules, safety, and multi-depot coordination were not explicitly incorporated. In addition, emissions were estimated based on total distance traveled, which may be refined using microscopic models that incorporate vehicle and route characteristics. The origin–destination matrix was derived from OSM data; although freely available, higher-resolution datasets, such as address-level demand, census-based population data, and detailed road networks, could improve accuracy.

Accordingly, future work could extend the framework to include time-dependent travel, real-time congestion, stochastic demand, richer labor and curb constraints, multi-echelon parcel–mail co-delivery, and improved emissions or cost modeling to enhance its strategic applicability. Finally, a digital twin could be developed by connecting the optimizer to real-time data (Thompson et al., 2026), enabling online deployment and broader decision-support capabilities.

Generative AI statement

The author(s) declared that generative AI was not used in the creation of this manuscript.

Any alternative text (alt text) provided alongside figures in this article has been generated by Frontiers with the support of artificial intelligence and reasonable efforts have been made to ensure accuracy, including review by the authors wherever possible. If you identify any issues, please contact us.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/ffutr.2026.1802453/full#supplementary-material

Supplementary material Supplementary material at the source

Footnotes

Edited by: Ching-Chiao Yang, National Kaohsiung University of Science and Technology, Taiwan

Reviewed by: Leila Nayibe Ramírez Castañeda, Universidad Libre, Colombia

Roberto Ulises Estrada Meza, University of Guadalajara, Mexico