A Simulation-Based ALNS Approach for Last-Mile Delivery with Uncertain Parcel Locker Capacities

2.1 Last-Mile Delivery

Last-mile delivery in city logistics could be viewed as a variant of the general VRP with various additional factors arising from customer requirements, traffic regulations, local policies, etc. The implementation of last-mile delivery mainly relies on home delivery, which requires the direct interaction between customers and deliveryman (Visser et al., 2014). However, the mode of home delivery may encounter delivery failure due to customer absence, and the vehicle utilization in home delivery mode is commonly low in practice (Allen et al., 2018). Hence, various new technologies have been proposed to enhance the last-mile delivery efficiency, such as the application of cargo-bikes, unmanned aerial vehicle (UAV) and robots. Cargo-bikes are used to deliver goods to clients, enabling shorter delivery routes by driving in areas restricted for motorized traffic and being less impacted by traffic conditions (Fikar et al., 2018). Simoni et al. (2020) considered the utilization of robots in last-mile delivery to form a truck-robot delivery system, and showed that the application of robots can improve last-mile delivery performance significantly, especially in some high population density areas. Yan et al. (2023) studied the application of UAVs for last-mile delivery in Taiwan with delivery routing and scheduling design. Cha et al. (2023) further investigated the collaborative delivery using ground vehicles and UAVs for practical logistics delivery purpose.

In addition to the application of new technologies, the design of network is also a nonnegligible issue in last-mile delivery in terms of the deployment of transit stations and front stations especially. For example, Zhou et al. (2018) designed a two-echelon VRP problem for last-mile delivery, in which customers are allowed to take their packages at some transit stations or parcel lockers. Zhao et al. (2021) proposed a location-routing problem for last-mile delivery where simultaneous pickup and delivery are considered and the collected items need to be processed in certain depots. Nieto Isaza et al. (2022) designed a crowdsource last-mile delivery system, which involves a two-stage decision-making process, i.e., the first stage is to make decisions regarding the location of small depots so as to allow flexible crowd-demand matching, and the second stage is to determine the demand allocation of for each deliveryman. Moreover, environmental and social considerations are also inevitable in last-mile city delivery. For example, Trott et al. (2021) proposed a mixed integer linear model for the last-mile delivery, which measures the vehicle emission and its affect to the sustainable development of cities. Caspersen and Navrud (2021) further investigated customer awareness towards the green initiatives and identified that customers prefer environmentally friendly last-mile delivery modes. Du et al. (2022) incorporated the impact of carbon emission into a two-echelon joint delivery problem, aiming to reduce both the operation cost and the adverse environmental impact. Besides, the operations of last-mile delivery could be also influenced by some governmental regulations especially in the context of epidemic (Sandoval et al., 2022).

2.2 Shared Parcel Locker

The application of shared parcel lockers is another promising last-mile delivery mode, which can complement the home delivery mode and improve the delivery efficiency considerably, especially in high population density areas. The work of Punakivi and Tanskanen (2002) was one of the initial studies which introduces the concept of shared reception boxes and estimates that the use of unattended shared reception boxes in last-mile delivery could reduce transportation costs by 55–66% compared to the manned premises with a delivery time window. In addition to the reduction of distribution cost, Edwards et al. (2010) claimed that the application of collection-delivery points was a convenient and environmentally friendly delivery mode, which can also reduce CO2 emissions substantially. Wang et al. (2014) once compared three delivery modes in last-mile delivery, i.e., attended home delivery (AHD), reception box (RB) and collection-and-delivery points (CDPs), and pointed out that the application of CDPs played a significant role in areas with high population density. Grabenschweiger et al. (2021) investigated the last-mile delivery system with locker box stations, in which heterogeneous locker boxes were considered and compensation payment was involved if customer chose to use locker box. Pinchasik et al. (2023) used empirical shipment data and real-world routing optimization procedures to demonstrate that increased parcel locker use can reduce logistics costs, traffic, and emissions from last-mile distribution of a logistics service provider.

The applicability of shared parcel lockers relies on multiple factors. On the one hand, the customer preference and deliveryman perception can affect the usage of shared parcel lockers significantly. Vakulenko et al. (2018) once studied customer perception of shared parcel lockers as a new delivery service and investigated the customer value benefited from the usage of shared parcel lockers. Yu et al. (2021) analyzed the value of shared parcel lockers for customers who choose pickup service and for deliverymen who simplify the delivery process by dropping off their deliveries at shared parcel lockers. Ghaderi et al. (2022) claimed that the application of shared parcel lockers in last-mile delivery can improve delivery efficiency in view of the increase of successful delivery rate. On the other hand, the deployment of shared parcel lockers in last-mile delivery network can influence its effectiveness and efficiency as a new delivery mode. Deutsch and Golany (2018) demonstrated the design of the optimal number, location, and size of parcel lockers can affect the last-mile delivery performance considerably. Jiang et al. (2020) once studied a last-mile vehicle scheduling problem considering both home delivery and locker delivery, in which customers can determine their preferred delivery modes. Lin et al. (2022) studied the location deployment of shared parcel lockers by predicting the likelihood of customers who prefer shared parcel lockers so as to maximize the possible revenues. Rossolov (2023) justified potential Alternative Parcel Location (APL) sites by developing behavioral models for existing choices between home deliveries and post office deliveries, considering time expenditures for the pick-up process and in-store grocery shopping.

2.3 Stochastic Delivery

Stochasticity is inevitable in last-mile delivery problems considering the practical delivery scenarios in large cities. Among the various stochastic factors, stochastic travel time and stochastic customer demand are two mostly studied aspects in delivery problems. Laporte et al. (1992) was one of the early studies which views the travel time as a stochastic factor and proposes the concept of stochastic vehicle routing problem, and they designed a chance-constrained programming (CCP) model to solve the proposed problem, which can ensure the “path success” within a specified confidence interval. Miranda and Conceição (2016) studied the stochastic vehicle routing problem with the estimation of the vehicle arrival time and arrival probability within service time windows so as to guarantee the delivery service level. In contrast, the vehicle routing problem with stochastic customer demand refers to the fact that the customer demand is not deterministic in advance (Dror, 1993). The delivery problem with stochastic demand is commonly modeled as a two-stage stochastic programming model, which is comprised of priori optimization strategies and recourse strategies (Goodson et al., 2012). Simulation-based approach is another popular method to tackle the delivery problem with stochastic demand (Liu et al., 2017).

In contrast to delivery problems with stochastic travel time or stochastic customer demand, this study considers a new variant in which uncertainty lies in the available capacities of shared parcel lockers. In practice, residual locker availability fluctuates noticeably, and shortages or under-utilization are frequently observed, which can significantly affect delivery schedules. Hence, the proposed problem not only enriches the theoretical landscape but also yields practically useful managerial insights for locker-aware last-mile planning.

Simheuristics address stochastic routing by combining deterministic search (e.g., ALNS) with Monte Carlo simulation (MCS) to evaluate candidate plans under realizations (e.g., Guimarans et al., 2018; Keenan et al., 2021). Prior applications predominantly model uncertainty at the customer (demand) or arc (travel time) level. Our setting is different as uncertainty acts on a shared facility (parcel-locker residual capacity), which couples multiple routes and customers and can induce execution-phase recourse. Accordingly, our method is problem-tailored as it integrates locker-aware destroy/repair and decoding with overflow penalties (as surrogates for reassignment/rollback), together with a two-phase MCS schedule, i.e., fast in-search screening and deep post-search evaluation of an elite set, and adaptive operator weighting calibrated to reduce overflow-induced add-ons. This integration promotes robustness during the search, not only ex post, and positions the approach within state-of-the-art simheuristic practice for logistics under uncertainty.

Distributional choices in stochastic delivery. In the stochastic-routing literature, bounded continuous distributions are commonly used to model uncertain availability or capacity (e.g., truncated normal, Beta, uniform/triangular variants), while empirical distributions are adopted when field data are available. Our choice of a truncated (bounded) normal reflects the aggregation of many small events and the physical upper bound, and it is interchangeable with other bounded families in our simulation-based evaluation. Hence the method and conclusions are not tied to Gaussianity.

Fuzzy optimization and MCDM (complementary view). Beyond probabilistic models and simheuristics, fuzzy optimization and multicriteria decision-making (MCDM) provide an alternative way to capture epistemic or linguistic uncertainty in logistics planning. Representative studies employ trapezoidal/intuitionistic fuzzy sets with centroid and graded mean ranking operators to aggregate expert assessments and rank policies or configurations; integrated MCDM schemes (e.g., AHP/TOPSIS/VIKOR variants under fuzzy information) are often used for facility siting, carrier/strategy selection, and resource allocation under competing criteria (Celik et al., 2014; Nghiem and Chu, 2024). While our work follows a probabilistic, simulation-based approach to model shared-facility uncertainty (stochastic locker capacity) and optimize routing decisions, fuzzy/MCDM frameworks are complementary as they can encode managerial preferences and imprecise information and are particularly useful when data are scarce or qualitative. In this sense, fuzzy methods and fuzzy-MCDM can inform locker deployment, capacity configuration, or the policy-level tuning that surrounds our operational routing model.

3.1 Problem Description

VRPSTC can be represented as a directed graph G = ( N , A ) , where N = 0 ∪ N P ∪ N C denotes the set of network nodes, in which { 0 } denotes the distribution center, N P and N C represents the set of shared parcel lockers and the set of customers. respectively. Moreover, N C = N C H ∪ N C P , which indicates customers with home delivery service requirements and locker pick-up service, respectively. A = { ( i , j ) | i , j ∈ N , i ≠ j } denotes the set of arcs, and each arc is associated with a delivery cost c i j and travel time t i j .

The distribution center { 0 } possesses a fleet of delivery vehicles K, and the maximum loading capacity of vehicle k is denoted as Q k . For each customer i ∈ N C , a certain delivery demand r i is known in advance. For customer i ∈ N C H , a corresponding time window [ e i , l i ] and service time s i are provided as well. For customer i ∈ N C P , the delivery vehicle needs to deliver the customer order to a certain self-pickup locker. Concerning the settings of the shared parcel locker p ∈ N P , the design capacity Q p and service time s p are provided firstly. However, in practical scenarios, the available capacity of the shared parcel locker is a random variable, as q ~ p , due to unexpected pickup and storage operations from customers. Moreover, the following assumptions are made to facilitate the modelling process of the proposed VRPSTC.

1)

The delivery service mode chosen by each customer is known in advance and the customer order cannot be split.

2)

The specific geographical distribution of shared parcel lockers is known, and lockers operate 24 h a day by default.

3)

The available capacity of shared parcel locker q ~ p is assumed to follow a normal distribution, as q ~ p ∼ N ( μ p , σ p 2 ) , with an additional lower and upper bound settings of [ 0 , Q p ] .

Rationale for the capacity distribution and implementation. Execution-time locker availability is the residual number of free slots after many asynchronous pick-ups and deposits since the last visit. Under mild aggregation effects, this quantity is well approximated by a continuous distribution whose support is naturally bounded by the physical capacity Q p . Accordingly, we model a truncated normal on [ 0 , Q p ] ; in simulation we further round q ~ p to integer slots. This is an engineering approximation that captures the aggregate volatility while respecting capacity bounds. Importantly, our framework is distribution-agnostic: any bounded D p (e.g., Beta scaled to [ 0 , Q p ] , truncated lognormal, or an empirical/bootstrapped D p ) can replace the default without changing the algorithm or evaluation logic. The qualitative insights we report depend mainly on μ p and σ p (Section 5.4), rather than on the specific family.

Execution timing and shared capacity. Routing and locker assignments are planned ex-ante. The available capacity of each parcel locker, q ~ p , is realized at execution time due to asynchronous pick-ups and deposits by third parties. When the realized capacity at a visited locker is insufficient for the parcels assigned to it, pre-specified decoding rules implement reassignment to backup lockers along the route or rollback to home delivery, with the resulting extra mileage/time penalized during simulation-based evaluation. Each locker's available capacity is a shared resource across routes and customers, which creates a facility-level coupling that distinguishes VRPSTC from models where uncertainty is local to customers or arcs.

Figure 2 presents an illustrative example of the collaborative last-mile delivery system. In route 1, customer 1, 2 and 3 choose home delivery mode, and customer 5 and 6 choose the shared parcel locker delivery mode. However, in route 3, due to the uncertainty of the remaining available capacity, the locker 11 is unable to serve both customers 12 and 13 simultaneously. Thus, in this delivery solution, the order of customer 13 is eventually served by the locker 14.

OPEN IN VIEWER

Figure 2.

An illustration of the last-mile delivery with home delivery and parcel locker delivery.

3.2 Mathematical Model

The following symbols and notations are adopted to form the mathematical model. Symbols

Description

Mathematical model

min E q ~ p [ ∑ k ∈ K ∑ j ∈ N ∑ i ∈ N d i j x i j k + ∑ p ∈ N P ∑ i ∈ N C P d i p y i p ]

(1)

Subject to

∑ k ∈ K ∑ j ∈ N x i j k = 1 , ∀ i ∈ N C H

(2)

∑ j ∈ N x i j k = ∑ j ∈ N x j i k , ∀ i ∈ N C H , k ∈ K

(3)

∑ k ∈ K ∑ j ∈ V x i j k = 0 , ∀ i ∈ N C P

(4)

∑ k ∈ K v i k = 1 , ∀ i ∈ N C

(5)

x i j k ≤ v i k , ∀ i ∈ N C , j ∈ V , k ∈ K

(6)

∑ t ∈ N p y i p = 1 , ∀ i ∈ N C P

(7)

∑ i ∈ N C P r i y i p ≤ q ~ p , ∀ p ∈ N P

(8)

∑ k ∈ K ∑ j ∈ N x 0 j k ≤ | K |

(9)

∑ j ∈ N x 0 j k = ∑ i ∈ N x i 0 k , ∀ k ∈ K

(10)

∑ i ∈ N c r i v i k ≤ Q k , ∀ k ∈ K

(11)

∑ k ∈ K v p k ≤ 1 , ∀ p ∈ N P

(12)

∑ k ∈ K ∑ j ∈ N x p j k ≤ 1 , ∀ p ∈ N p

(13)

∑ j ∈ N x p j k = ∑ j ∈ N x j p k , ∀ p ∈ N P

(14)

v p k ≥ 1 − M ( 2 − y i p − v i k ) , ∀ i ∈ N C P , p ∈ N P , k ∈ K

(15)

e i ∑ j ∈ N x i j k ≤ w i k ≤ l i ∑ j ∈ N x i j k , ∀ i ∈ N C H , k ∈ K

(16)

w 0 k = 0 , ∀ k ∈ K

(17)

w j k ≥ w i k + t i j + s i − M ( 1 − x i j k ) , ∀ i , j ∈ N , k ∈ K

(18)

x i j k ∈ { 0 , 1 } , y i p ∈ { 0 , 1 } , z p ∈ { 0 , 1 } , v i k ∈ { 0 , 1 } , w i k ≥ 0

(19)

The objective function (1) is to minimize the sum of total delivery distance for couriers and the self-pickup distance for customers who choose the parcel locker delivery mode. Constraints (2) and (3) indicate that for customers with home delivery mode, the vehicle inflow and outflow are balanced. Constraint (4) means that for customer with locker pickup mode, no vehicles will visit their places. Constraint (5) shows that the demand of each customer is served by only one vehicle. In other words, no order split is allowed for all customers. Constraint (6) represents the relationship between the route variable x i j k and the assignment variable v i k . Constraint (7) enforces that only one parcel locker is assigned for each customer with locker pickup mode. Constraint (8) means that the total allocated number of orders to each parcel locker cannot exceed its remaining available capacity. Constraint (9) and (10) means that the number of vehicles dispatched from the distribution center cannot exceed the fleet size, and all dispatched vehicles must return to the distribution center once the delivery task is completed. Constraint (11) represents the capacity limit of vehicles. Constraint (12) means that maximumly one vehicle is allowed for each parcel locker. And constraint (13) and (14) represents the vehicle inflow and outflow balance for each parcel locker. Equation (15) represents the relationship between v p k , y i p and v i k , which is to maintain the allocation consistency among customers, parcel lockers, and vehicles. Constraint (16) represents the time window settings for each customer with home delivery mode. As for the distribution center, the initial starting time is set as 0 as shown in constraint (17). Constraint (18) indicates the differentiation of service starting time for node i and j provided that they are served by the same vehicle, which can be leveraged to break the possible subtour phenomenon. Constraint (19) are the settings of different decision variables.

Model positioning and distinctive features (interpretive note). The formulation captures a mechanism that is underrepresented in stochastic VRPs: uncertainty acts at a shared facility, i.e., the residual capacity of parcel lockers, rather than at the customer or arc level. Because the available capacity of a locker is a common pool jointly affected by multiple routes and customers, the shared-capacity constraint couples their feasibility conditions. When shortfalls occur upon execution, this coupling naturally induces recourse, namely reassignment to alternative lockers or rollback to home delivery. The model therefore jointly decides vehicle routing and locker assignment while coordinating the two delivery modes under a shared facility constraint. In terms of timing, plans are generated before capacity is realized; the realization of locker availability occurs at execution time, and decoding rules with overflow penalties act as proxies for operational recourse. This simulation-based evaluation (Monte Carlo simulation with overflow penalties) provides a risk-aware service-level surrogate without committing to a specific probability threshold. The above mechanism distinguishes the formulation from classic stochastic-demand or stochastic-travel-time VRPs, where uncertainty is local to customers or arcs and does not create a facility-level pooling constraint.

In this research, ALNS is adapted to solve the deterministic problem, which is one of the most well-acknowledged algorithms to solve various vehicle routing problems (Dayarian et al., 2016; Keskin et al., 2021; Windras Mara et al., 2022). Figure 4 illustrates the flowchart of the ALNS algorithm. Initially, a simple and efficient heuristic method is used to construct an initial solution, thereby starting the search process of the ALNS algorithm. In the following searching phase, the initial solution is first destroyed and then repaired using selected operators, leading to an updated solution. Then, the updated solution is evaluated for updating the global optimal solution and the current solution. Meanwhile, relevant operators are evaluated and scored, and the weights of these operators are updated after a certain number of iterations. The above process is repeated until the termination condition is met.

OPEN IN VIEWER

Figure 4.

The flowchart of the ALNS algorithm.

4.1.1 Initial Solution Generation

The ALNS starts from a randomized initial plan S 0 . The node set contains home-delivery customers N C H , parcel-locker customers N C P , and lockers N P . We generate S 0 in two steps:

Step 1 – Route construction (randomized). Create a list by shuffling N C H ∪ N C P ∪ N P . Split it into routes by scanning the list and opening a new route when the vehicle capacity would be violated. Within each route, keep the order induced by the scan (this yields a simple but diverse starting solution).

Step 2 – Decoding under deterministic locker capacity. Given route k with sequence R k , let C k P L = { i ∈ R k ∩ N C P } denote PL customers on the route in their visiting order, and P k = { p ∈ R k ∩ N P } the lockers appearing on the route in order. For each locker p ∈ P k set available capacity to its deterministic value Q p . Traverse C k P L in order; for each PL customer i, assign i to the next locker ahead on the route (in P k ) that still has remaining capacity. If none has enough capacity, mark i as unserved and apply an overflow penalty (or, equivalently, a rollback to home delivery). After decoding, remove unused lockers and unattached PL customers from the route representation. The same decoding procedure is used in simulation by replacing Q p with the realized capacity q ^ p .

Illustrative cases (Figure 5). Let Route 1 contain PL customers { 5 , 6 , 7 , 8 } and lockers { 9 , 10 } .

(a)

If Q 9 ≥ r 5 + r 6 + r 7 + r 8 , only locker 9 is used; locker 10 and PL nodes are removed after decoding.

(b)

If Q 9 < r 5 + r 6 + r 7 + r 8 but Q 9 + Q 10 ≥ r 5 + r 6 + r 7 + r 8 , then customers are split across lockers 9 and 10 according to the route order.

(c)

If Q 9 + Q 10 < r 5 + r 6 + r 7 + r 8 , leftover demand is unserved and penalized (e.g., customer 8).

(d)

If a route contains only HD customers and lockers (no PL customers), lockers are removed since they provide no service on that route.

OPEN IN VIEWER

Figure 5.

Initial solution and decoding under locker capacity. (Decoding assigns PL customers to the next locker with remaining capacity along the route; unused lockers/PL nodes are removed. The same rules apply in MCS with realized q ^ p ).

We use the same decoding logic within Monte Carlo simulation by sampling q ^ p and repeating the assignment.

4.1.2 Destroy Operators

Let current solution S contain route set K . Let | S | be the number of customer/locker nodes in S. Let the removal quota be q = ρ | S | with ρ ∈ ( 0 , 1 ) (Table 1 uses ρ = 0.1 ). Removed nodes are stored in list R .

(a)

Random removal. Uniformly pick q nodes from S and push to R .

(b)

Random-route removal. Uniformly pick a non-empty route k ∈ K ; remove all its nodes.

(c)

Shaw removal. Similarity between nodes i , j is ϕ ( i , j ) = α d ^ i j + β r ^ i j + γ e ^ i j , where d ^ i j = d i j max d , r ^ i j = | r i − r j | max r , e ^ i j = | e i − e j | max e are normalized distance/demand/earliest-time differences (Table 1 uses α = 6 , β = 4 , γ = 5 ). Seed with a random node i; iteratively remove j = arg min j ∈ R ϕ ( i , j ) until | R | = q .

(d)

Worst removal. Node-removal gain Δ f ( i ; S ) = f ( S ) − f ( S ∖ { i } ) , remove i * = arg max i Δ f ( i ; S ) and repeat until | R | = q

(e)

Worst-distance removal. For node i with predecessor p and successor n, Δ d ( i ) = d p , i + d i , n − d p , n , remove iteratively the node with largest Δ d ( i ) until q nodes are removed.

OPEN IN VIEWER

Table 1.

Parameter Settings.

Parameter description	Value
Maximum number of iterations ( N )	5000
Number of iterations for weight update ( m )	100
Adaptive weighting response factor ( θ )	0.1
Cooling rate ( μ )	0.998
New global optimal solution score ( σ 1 )	10
Better solution score ( σ 2 )	4
Worse solution but accepted score ( σ 3 )	2.5
Removal ratio ( α )	0.1
First weight for similarity calculation ( λ 1 )	6
Second weight for similarity calculation ( λ 2 )	4
Third weight for similarity calculation ( λ 3 )	5
Smoothing parameter ( λ )	0.1
Number of simulation times for fast evaluation ( N s 1 )	500
Number of simulation times for extensive evaluation ( N s 2 )	3000
Size of solution list ( M )	10

4.1.3 Repair Operators

Candidate positions that violate vehicle capacity or HD time windows are discarded. Time is checked by forward-propagating service start times; lockers have no time window by default.

Let S ⊕ ( i , p o s ) denote inserting node i at position pos on some route; insertion cost is the objective increment Δ f ( i , p o s ; S ) = f ( S ⊕ ( i , p o s ) ) − f ( S ) .

(a)

Greedy insertion. For each i ∈ R , compute min p o s Δ f ( i , p o s ; S ) ; pick i * = arg min i min p o s Δ f ( i , p o s ; S ) , insert at its best position; repeat.

(b)

Regret-2 insertion. For each i, let Δ 1 ( i ) ≤ Δ 2 ( i ) be the best and second-best feasible increments. Regret value Γ 2 ( i ) = Δ 2 ( i ) − Δ 1 ( i ) . Pick i * = arg max i Γ 2 ( i ) and insert at its best position. (Regret- k is analogous.)

(c)

Distance-greedy insertion. For each i, for each arc ( p , n ) , define Δ d ( i ; ( p , n ) ) = d p , i + d i , n − d p , n ; select the pair i , p o s with minimum Δ d .

(d)

Greedy-with-noise. Perturb greedy costs: Δ f ~ ( i , p o s ) = ( 1 + θ ε ) Δ f ( i , p o s ) with ε ∼ U [ − 1 , 1 ] and θ = 0.1 , then apply rule (a).

(e )

Distance-greedy-with-noise. Apply the same perturbation to Δ d .

After reinsertion, decoding assigns PL customers to lockers along the route with remaining capacity; overflow triggers the penalty term.

4.1.4 Adaptive Mechanism

Let D and R be destroy/repair operator sets; each operator o has weight w o and cumulative score π o initialized as w o = 1 , π o = 0 .

Route wheel section. At iteration t , Pr ( o chosen ) = w o ∑ o ′ w o ′ .

Scoring. After applying a pair ( d , r ) and obtaining S ′ from current S, assign an event score σ to both d and r based on the outcome:

σ = { σ best = 10 f ( S ′ ) < f ( S * ) ( new global best ) ; σ imp = 4 f ( S ′ ) < f ( S ) ( improvement ) ; σ acc = 2.5 f ( S ′ ) ≥ f ( S ) but accepted ( SA ) ; 0 otherwise ,

and update π d ← π d + σ , π r ← π r + σ . Count usages u o within the current segment of ψ iterations (Table 1: ψ = 100 ).

Weight update at the end of segment. With smoothing factor λ ∈ ( 0 , 1 ) (Table 1: λ = 0.1 ),

w o ← ( 1 − λ ) w o + λ π o max { 1 , u o } , ∀ o ∈ D ∪ R .

Then reset π o ← 0 , u o ← 0 to encourage diversification.

The above setting is the canonical ALNS adaptation scheme; our scores ( 10 , 4 , 2.5 ) , ψ = 100 and λ = 0.1 fall within commonly used ranges and match Table 1.

4.1.5 Acceptance and Stopping Criteria

SA acceptance. With Δ = f ( S ′ ) − f ( S ) and temperature T t

p acc = { 1 Δ < 0 exp ( − Δ / T t ) Δ ≥ 0

T t = T 0 ⋅ α t with T 0 the initial temperature and α = 0.998 (Table 1).

Stopping. Stop at t = 5000 (Table 1); a time limit can be used equivalently.

4.2 Monte Carlo Simulation

ALNS searches for high-quality deterministic plans, while MCS evaluates their performance under stochastic locker-capacity realizations. In this work, every newly accepted solution during the ALNS search is evaluated by MCS rather than only ex-post selection.

Two-phase MCS schedule. Each newly accepted solution is first evaluated with a small number of Monte Carlo replications to obtain a fast estimate; promising ones are inserted into an elite set. After the ALNS loop terminates, solutions in the elite set undergo a larger number of replications for deep evaluation. The final plan is selected by expected cost E ^ [ cost ( S ) ] , with tiebreakers preferring fewer/shorter overflow-related penalties. This budgeted schedule balances runtime and stability and steers the search toward locker-aware robust plans.

Capacity sampling and decoding (plug-in). In each replication m, for each locker p, we draw one realization q ^ p ( m ) ∼ D p , shared across all routes within the same replication. By default, D p is a truncated (bounded) normal on [ 0 , Q p ] (i.e., clip a Gaussian to the physical capacity) and the realized capacity is rounded to integer slots; however, the framework is distribution-agnostic and any bounded alternative (e.g., Beta scaled to [ 0 , Q p ] , truncated lognormal, or an empirical/bootstrapped D p ) can be used without changing the algorithm. Decoding then assigns parcel-locker customers to the next locker with remaining capacity along the route; any shortfall triggers overflow penalties.

Realization-based objective and interpretation. For each replication m, the execution cost of plan S is computed as cos t ( m ) ( S ) = T exec ( m ) ( S ) + Π m ( S ) , Π m ( S ) = ∑ p ∈ N P ∑ k π Δ p , k ( m ) ( S ) , where T exec ( m ) ( S ) is the travel time under realized capacities, Π m ( S ) is the overflow-induced recourse penalty, Δ p , k ( m ) ( S ) is the capacity shortfall at locker p on route k under replication m, and π monetizes reassignment/rollback as a proxy for the delay or unmet lock allocation. The Monte Carlo estimator E ^ [ cost ( S ) ] is therefore our primary realization-based robustness indicator. Accordingly, STC1 evaluates the deterministic plan S det under realization (i.e., S det with MCS), while STC2 is the best uncertainty-aware plan selected from the elite set by E ^ [ cost ( ⋅ ) ] .

The details of MCS implementation are presented as follows.

 MCS implementation

  Notations

   N s Number of simulations, indexed using m

   S The current solution

   K Number of routes in S, indexed using k

   N k Number of nodes in route k, indexed using i

   n i k The i th node of route k

   T C Total cost

 1) Initialize m = 0 , k = 0 , i = 0 , T C = 0

 2) While m ≤ N s

   // sample once per locker (shared across routes)

   For each p ∈ N P , do q ^ p m ∼ D p //e.g., truncated N ( μ p , σ p 2 ) on [ 0 , Q p ]

   While k ≤ K

    While i ≤ N k

     If n i k is a locker node p, then use q ^ p //no re-sampling inside the replication

      i = i + 1

   Decode the route k //apply reassignment/rollback if overflow occurs; add penalties)

   Calculate the cost of route k

   T C = T C + c o s t S

   m = m + 1

 Output the stochastic performance of solution S as E ( c o s t S ) = T C N s

5.1 Instance Generation and Parameter Settings

Test instances are derived from the classical Solomon benchmark dataset as (Solomon, 1987) with additional modifications to suit the features of this research problem. In general, given N customer nodes, the number of customers with home delivery and self-pickup requirements are set as N C P = 0.4 * N and N C H = 0.4 * N , and the left are treated as the shared parcel lockers as N P = 0.2 * N . Moreover, the design capacity of one shared parcel locker is set as 100 by default. It is noted that the settings of ALNS parameters may affect the algorithmic performance considerably. Hence, some preliminary experiments are conducted in order to find the appropriate parameter settings. Table 1 summarizes the relevant parameter settings.

Anchoring of mean/variance and bounds. The parameter grids for locker availability used in Sections 5.3–5.4 were selected to span low/medium/high regimes of mean availability together with representative levels of dispersion, so as to reflect typical utilization patterns and peak-time volatility. The ranges are anchored by an independent field data collection on parcel-locker usage and by standard operational constraints (e.g., physical upper bounds Q p and integer slots). Our goal is to cover realistic operating regimes rather than to fit a single instance-specific distribution.

5.2 Evaluating the Performance of ALNS

The performance of the ALNS is firstly examined and validated in comparison with the commercial optimization solver Gurobi 9.1.2 over a number of test instances. In order to facilitate the solving process, the time limit of Gurobi is set as 1800 s. Moreover, for each instance, the proposed ALNS is applied 10 times, and the average objective value and CPU times are collected. Table 2 summarizes the comparison result of ALNS and Gurobi over instances with 10, 20, and 30 nodes. It is noted that for instance with 10 nodes, both ALNS and Gurobi can find the optimal solutions within a short time. When the scale of instances reaches 20 nodes, Gurobi can find feasible solutions for most cases with the time limit setting of 1800 s. In contrast, the computational time of ALNS for solving instances with 20 nodes is around only 20 s, which distinctly dominates that of Gurobi. For instances with 30 nodes, although Gurobi can still find feasible solutions for some instances, ALNS can find better solutions over Gurobi using only 50 s averagely. Gurobi solves the deterministic surrogate (locker capacity fixed at E [ q ~ p ] and serves as a small-scale reference for solution quality and runtime; for larger instances we focus on ALNS because the deterministic MIP escalates quickly, and deterministic solvers do not optimize the stochastic objective used in deployment.

Instances with larger scales are also tested to further assess the performance of the proposed ALNS. As shown in Table 3, the devised ALNS is capable of solving instance with 40, 50, and 60 nodes using 113, 235, and 402 s averagely. By contrast, Gurobi cannot find any feasible solutions within 1800 s. Moreover, it is noticed that the ALNS possesses a relatively stable performance in terms of a small gap between the best objective value and the average objective value. Thus, it is concluded that the proposed ALNS performs well in terms of both computational efficiency and solution quality, especially when the scale of instances becomes large.

OPEN IN VIEWER

Table 3.

The Performance of ALNS Over Large-Scale Instances.

Scales	40			50			60
Instances	Best	Average	CPU (s)	Best	Average	CPU (s)	Best	Average	CPU (s)
C102	589.48	598.36	83.52	689.43	719.96	168.47	811.66	842.07	278.83
C103	581.49	599.59	80.03	696.89	716.44	163.82	811.26	851.80	270.29
C105	579.60	601.59	77.77	701.75	722.51	166.82	813.41	859.35	281.36
C106	586.43	598.54	86.70	703.26	725.51	166.40	822.61	873.88	280.29
C107	580.14	595.56	84.78	694.94	720.99	165.19	822.04	865.11	275.94
R102	624.69	646.40	94.37	883.75	943.10	191.52	882.92	929.74	321.30
R105	614.22	654.88	97.62	831.05	949.62	199.00	854.76	891.66	294.00
R106	619.25	646.15	92.21	874.94	935.49	194.82	826.39	896.55	286.42
RC103	1025.91	1055.79	86.46	865.19	955.69	152.90	1138.64	1352.88	263.39
RC104	1005.93	1044.71	84.02	939.83	982.85	150.70	1171.12	1297.66	258.76
RC105	997.03	1037.75	85.99	892.19	986.81	158.51	1260.05	1354.99	261.17
RC107	987.92	1042.57	84.56	927.92	972.77	152.57	1152.08	1325.45	264.75
RC108	1000.88	1063.58	84.03	834.59	952.27	149.97	1221.53	1349.80	264.11
Average	753.31	783.50	86.31	810.44	868.00	167.75	968.34	1053.15	276.97

In this sub-section, the available capacity of each locker q ~ p is randomly generate following a normal distribution as q ~ p ∼ N ( 45 , 100 ) , which is consistent with the initial parameter settings. It is noted that due to the uncertain available capacity of the shared parcel locker, the reassignment between customers and the shared parcel locker may occur in case of insufficient available capacities. In this case, the application of MCS is implemented to evaluate the stochastic performance of a given solution. As described previously, when the ALNS reaches its termination criterion, ten promising solutions are stored as potentially optimal solutions. Then, a comprehensive and extensive MCS procedure is applied to further assess these solutions accurately, and eventually the solution with best performance is determined as the best stochastic solution. Figure 6 illustrates a detailed comparison between best deterministic solution and the best stochastic solution. Noted that the distribution routes between deterministic and stochastic scenarios vary significantly, which intrinsically attributes to the uncertain available capacity of the shared parcel lockers. In addition, the matching between customers and the shared parcel locker changes dramatically as well.

OPEN IN VIEWER

Figure 6.

The comparison of optimal solutions in deterministic and stochastic scenarios.

Table 4 summarizes a performance comparison between deterministic and stochastic scenarios over some test instances, in which DTC1 represents the objective value of the optimal deterministic solution, STC1 indicates the expected objective value of the deterministic solution with MCS evaluation, and STC2 represents the expected objective value of the optimal stochastic solution. Note that the solution of DTC1 is acquired using the expected available capacity as a deterministic input, and this corresponding solution is further assessed through MCS so as to obtain the expected objective value in STC1. Simulation results show that the value of STC1 exceeds that of DTC1 to the extent of 73.51% on average, which indicates that the ignorance of uncertainty regarding the available capacity of the shared parcel locker could lead to a huge operational cost deviation. The difference between STC1 and STC2 suggests that the optimal stochastic solution and the optimal deterministic solution might be also slightly different, which is practically consist with the principle of sim-heuristic framework as the optimal deterministic solutions are likely to the promising solutions in stochastic scenarios. Indeed, the optimal stochastic solution and the optimal deterministic solution in RC101 is the same, and the slight difference between STC1 and STC2 regarding RC101 is due to the randomized MCS procedure. Figure 7 visualizes the difference between STC1 and STC2 over these instances. Across instances, executing the deterministic plan (DTC1 → STC1) yields a substantially higher expected cost due to locker shortfalls, while the simheuristic solution (STC2) consistently narrows this gap by reducing reassignment/rollback incidence.

OPEN IN VIEWER

Figure 7.

The comparison between deterministic and stochastic solutions.

OPEN IN VIEWER

Table 4.

Performance Comparison Between Deterministic and Stochastic Solutions.

# of nodes	Instances	DTC1	STC1	STC2	Gap (%)
40	C101	770.83	1015.33	985.83	2.91%
	C102	689.73	1071.91	877.73	18.12%
	C106	679.08	1291.97	873.19	32.41%
	C107	709.97	1389.92	934.13	32.79%
	R101	767.25	1118.52	1022.43	8.59%
	R102	845.75	1807.12	1111.05	38.52%
	R105	828.13	1521.55	1095.54	28.00%
	R106	853.46	1393.26	1043.47	25.11%
	RC101	1291.49	2017.72	2002.16	0.77%
	RC102	1405.50	2404.14	2153.22	10.44%
	RC103	1429.87	2608.86	2278.05	12.68%
	RC105	1373.32	2441.88	2100.54	13.98%
	RC106	1351.60	2466.98	2132.72	13.55%
Average		999.69	1734.55	1431.54	18.3%

Beyond the deterministic objective (DTC1), the gap to STC1 quantifies the expected execution penalty caused by locker shortfalls when the deterministic plan runs under capacity realizations (Section 4.2). The reduction from STC1 to STC2 captures the robustness gain obtained by screening plans across realizations within the simheuristic and selecting by E ^ [ cost ( ⋅ ) ] . Hence, Table 4 already provides realization-based evidence: by curbing capacity shortfalls (monetized via overflow penalties), the uncertainty-aware plan narrows the plan–execution gap and lowers expected recourse cost, consistent with the reassignments visible in Figure 6.

In this subsection, the impact of the available capacity of the shared parcel locker is further investigated in terms of different parameter settings of the mean and variance. As described previously, the design capacity of the shared parcel locker is set as 100 by default, and normal distribution is assumed to simulate its remaining available capacity, as q ~ ∼ N ( μ , σ 2 ) . Table 5 presents the detailed outcome regarding different settings of μ and σ . Firstly, it is noticed that when the mean available capacity is small, e.g., μ = 40 and μ = 50 , the total cost incurred is significantly over that in case of μ = 60 , 70 , and\; 80 , which is mainly due to the either the irregular matching between customers and lockers or the insufficient capacity violation. Secondly, it is observed that in case of small available capacity, as μ = 40 and μ = 50 , the total cost increases slightly along with the enlargement of its corresponding variance from σ = 5 to σ = 15 . However, such an observation does not repeat when the expected available capacity becomes large. In other words, the total cost fluctuates without any specific regular patterns in case of μ = 60 , 70 and 80 . The above phenomenon may attribute to the penalty settings in case of the available capacity violation as the violation probability with small expected available capacity is considerably higher than that with large expected available capacity. Figure 8 further presents a straightforward illustration regarding the total cost changes with different mean and variance settings.

OPEN IN VIEWER

Figure 8.

Total cost changes with different fluctuations of available capacities.

OPEN IN VIEWER

Table 5.

Sensitivity Analysis Regarding the Variance of the Available Capacity of Shared Parcel Lockers.

	μ = 40			μ = 50			μ = 60			μ = 70			μ = 80
Instance	σ = 5	σ = 15	σ = 25	σ = 5	σ = 15	σ = 25	σ = 5	σ = 15	σ = 25	σ = 5	σ = 15	σ = 25	σ = 5	σ = 15	σ = 25
C101	1266.9	1669.4	1910.2	798.73	926.25	1538.4	777.33	756.62	915.55	746.42	781.32	774.79	711.51	724.82	741.16
C102	1139.0	1574.0	1756.5	683.42	846.02	1070.8	676.21	693.12	820.61	651.21	700.44	715.49	620.20	649.94	688.13
C106	1149.3	1593.8	1747.9	705.77	859.69	1054.9	694.70	692.58	791.03	688.32	666.81	684.70	646.97	683.88	667.82
C107	1172.2	1616.6	1786.7	666.62	846.20	1099.4	664.08	727.90	753.39	680.72	678.43	701.41	635.85	661.48	665.13
R101	1052.91	1635.43	1746.09	814.75	1078.95	1245.98	828.44	919.49	1073.91	807.47	813.71	913.13	757.32	816.49	789.78
R102	1065.77	1642.39	1872.47	889.21	1015.59	1432.34	778.46	820.07	995.30	770.48	852.93	923.33	702.21	833.60	892.72
R105	1082.02	1660.21	2077.47	887.01	1082.10	1281.00	828.36	817.51	1007.83	782.47	760.95	950.91	745.79	793.36	794.21
R106	1088.61	1648.27	2016.57	875.12	1053.38	1313.42	818.40	930.72	943.41	727.57	836.12	886.82	788.93	806.61	817.74
RC101	2470.6	2876.1	2890.8	1402.0	1648.3	2111.31	1257.14	1286.64	1789.74	1187.28	1263.31	1355.11	1182.52	1213.19	1228.63
RC102	2444.6	2843.0	2971.9	1432.7	1925.1	2125.9	1224.50	1427.58	1637.45	1172.14	1331.57	1406.83	1257.40	1259.10	1242.69
RC103	2502.7	2999.7	3281.6	1468.2	2077.1	1979.5	1359.94	1402.20	1873.29	1301.71	1216.11	1470.06	1310.16	1226.31	1262.48
RC105	2842.0	3123.2	3208.0	1435.4	1803.8	2314.4	1345.07	1316.65	1561.50	1236.94	1285.67	1442.27	1215.17	1186.02	1306.08
RC106	2691.0	2879.6	3009.0	1557.9	1858.2	2138.4	1397.89	1362.67	1576.26	1313.41	1239.62	1349.32	1230.57	1199.37	1248.48
Average	1964.2	2352.8	2506.9	1127.8	1421.1	1714.7	1044.0	1073.9	1302.0	997.5	1018.1	1099.9	978.9	978.2	1005.6

Across instances, mean availability dominates the cost response: when the expected locker capacity is tight, total cost increases markedly relative to cases with ample mean availability (Table 5; Figure 8). In contrast, the impact of variance is conditional on the mean. Under tight means, higher variance systematically raises cost because shortfalls (and thus reassignment/rollback) become more frequent. When the mean is ample, variance has a much weaker and less systematic effect. Operationally, this pattern suggests that raising baseline availability (e.g., adding cells or redistributing loads) is the first-order lever, while variance mitigation (e.g., visit timing or micro-rebalancing) becomes relevant only when means cannot be lifted. These findings complement the plan–execution evidence in Section 5.3: reducing the frequency of capacity shortfalls directly limits recourse actions and narrows the plan–execution gap between deterministic plans and stochastic realizations. These mean-dominant patterns persist across bounded distributional families with matched ( μ , σ ) (e.g., truncated normal vs. Beta), since the shortfall probability, and thus recourse frequency, is driven primarily by the mean.