Optimizing routing for autonomous delivery and pickup vehicles in three-dimensional space

1.   Introduction

With the rapid development of the logistics industry, the demand for express delivery is growing significantly creating challenges in the last mile of logistics, such as difficulty in effectively coordinating terminal distribution resources with dispersed, small-batch, and high-frequency consumer demand in urban areas ^[1,2]. Autonomous delivery robots (ADRs) are innovative tools to address this last-mile problem by reducing the freight cost ^[3,4]. ADRs can utilize facilities including staircases and elevator networks to achieve door-to-door delivery across multi-level urban infrastructures ^[5].

The vehicle routing problem (VRP) framework is introduced to solve the last-mile logistics problem using ADRs. By optimizing routes, VRP reduces delivery distances and minimizes customer waiting times ^[6]. Many existing studies have applied VRP into real delivery tasks ^[7,8]. Most focus on trucks deliveries while overlooking the advantages of ADRs, which can offer door-to-door service. A distinguishing feature of ADRs is their handling of goods with 3D characteristics (length, width, and height). Wang et al. considered this feature in their modeling constrains, but this resulted in a complex two-stage model ^[9]. However, for campus deliveries, package sizes are standard. Thus, ADRs containers can be split into smaller cells designed for the standard package sizes. This design can help merchants quickly load goods and allows customers to pick up their items from designated containers ^[10]. In this case, the best way is to transform 3D spatial constraints into one-dimensional matching constraints, thereby reducing model complexity ^[11].

The ability of door-to-door delivery provides a strong foundation for extending the VRP framework into 3D environments in real campus settings ^[12]. The complexity of real-world 3D delivery scenarios, characterized by multi-building distributions and vertical floor transitions, presents significant operational challenges ^[13]. Compared with the previous studies involving two-dimensional environments, usually using Euclidean or Manhattan distances ^[14], we enhance the 3DCityNet architecture through a rigorous 3D spatial parameterization ^[15]. Our algorithm generates a 3D routing matrix that simultaneously accounts for: horizontal street-level networks, vertical circulation systems (elevators/stairwells), and house layout in the same floor.

In this paper, we propose an innovative method to fulfill goods delivery. ADRs are employed in multi-floor buildings to offer door-to-door service in the campus. This paper introduces the Discrete Capacity Vehicle Routing Problem with Simultaneous Delivery-Pickup and Soft Time Windows (DCVRP-SDP-STW) model, which can capture the features of a real-world ADR container. Given the NP-hard nature of the VRP, managing large-scale instances remains a significant challenge ^[16,17]. This paper proposes the Spatiotemporal Multi-Stage Clustering Tabu Ant Colony Search Algorithm with Temporal Pheromones (SMC-TACSATP) to efficiently obtain high-quality solutions. The contributions of this paper are as follows:

● We propose a novel VRP variant that considers both supply and demand features. This model introduces a discrete capacity to account for the physical characteristics of ADRs containers.

● We propose the SMC-TACSATP algorithm for large-scale delivery scenarios. Such a method can provide high-quality solutions within a limited time frame by utilizing a clustering approach to capture spatiotemporal characteristics and proposing an improved ant colony optimization (ACO) algorithm.

● We conduct a case study using map data. An algorithm is developed to calculate 3D distances between points, and the dummy node method is introduced to enable split delivery, which, for the first time, applies VRP into last-mile delivery of goods in urban 3D buildings.

The remainder of this paper is organized as follows: Section 2 reviews related literature. Section 3 presents the optimization model and methods for its simplification, while Section 4 describes the solution algorithm. Section 5 provides a case study, demonstrating the application and advantages of the proposed model and solution algorithm. Finally, Section 6 concludes the paper.

2.   Literature review

One of the main challenges in this paper is how to optimize the VRP under complex constraints ^[18,19]. The first one is modelling the VRP with distinct good sizes in vehicles. A classic variation is the capacitated vehicle routing problem (CVRP), which can be divided into two types. The first type, CVRP, was put forward by Clark and Wright, and assumed identical vehicles with equal capacity ^[20]. The second type heterogeneous fleet VRP, proposed by Quak and de Koster, which considers varying the types of vehicles ^[21]. Lin et al. and Luiz Usberti et al. considered trucks and trailers as part of the fleet ^[22]. In the past decade, VRP involving loading constraints for real containers has gained huge attention regarding the sizes (e.g., length, width, and height) of goods ^[9].

The second challenge is to consider customer-specific time windows in the VRP model. Solomon first considered the VRP with time window (VRPTW) ^[23]. When vehicles arrive before the desired time, they must wait; conversely, arriving late is prohibited. Moreover, Balakrishnan set soft time windows, which allow arriving late at an additional cost, thus adding flexibility to routing decisions ^[24]. Guertin et al. extended this concept by transitioning from static to dynamic scheduling under soft time window constraints ^[25]. Hashimoto et al. contributed a dynamic programming framework for situations involving convex and non-convex penalties ^[26].

Another significant modeling challenge lies in involving both deliveries and pickups. To handle such a problem, Angelelli and Mansini considered the simultaneous pickup and delivery, enabling a vehicle to achieve pickup and delivery tasks within a single trip ^[27]. Goetschalckx and Jacobs-Blecha considered VRP with backhauls, which simplifies the problem by enforcing a strict sequence where all deliveries precede pickups ^[28].

Finding the optimal solution for the VRP is challenging due to its NP-hard nature. One of the foundational methods was introduced by Clark and Wright, which focuses on combining individual routes after ensuring all operational constraints are satisfied ^[20]. Altabeeb et al. also introduced a firefly algorithm to tackle capacity constraints ^[29]. For VRP with time window, Gambardella et al. were the first to successfully apply the ACO algorithm to VRPTW ^[30]. Wang et al. further enhanced this approach by designing a hybrid multi-ant system to solve the VRPTW with service time constraints ^[30]. Zhang et al. integrated time penalties into the updating process of pheromone ^[31]. A typical VRP model, which is similar to our model for ADRs, also addressed this model using ACO ^[32].

The VRP in a 3D environment is a newly-raised problem. Door-to-door delivery belongs to a 3D vehicle routing problem, in which ADRs travel between buildings. Traditional VRP models simply assumed that the nodes are in the same plane and used time or fuel consumption as the weight between nodes ^[33]. To overcome this limitation, Thill et al. put forward 3DCityNet, which integrates an object-oriented GIS model to seamlessly connect indoor and outdoor transportation networks ^[15].

How to solve large-scale optimization is a challenging topic ^[34]. Recent studies show promise in combining machine learning models with the heuristic algorithm to solve the VRP efficiently ^[35]. To solve VRP with time window models, Wang et al. incorporated time information into the distance metric and proposed the 3D k-means approach ^[36]. Erdogan and Miller-Hooks used the DBSACN clustering algorithm to partition nodes and quickly generate an initial solution ^[37].

The following Table 1 lists the literature discussed in this paper, and the methods used.

3.   Problem description

This paper aims to minimize the weighted total distance, number of vehicles, and time window penalties.

The operational area is a campus environment with $|I|$ nodes, comprised by a logistics center and multiple customer nodes. The travel time between node $i \in I$ and node $j \in I$ is denoted by ${g_{ij}}$. Among these nodes, we define $P$ as the subset with pick-up demands and $D$ as those requiring deliveries. There are $|V|$ vehicles in depot to conduct pick-up and delivery tasks. Each customer $i \in I$ sets a time window $[t{e_i}, t{l_i}]$ where $t{e_i}$ is the earliest arriving time permitted and $t{l_i}$ is the latest. If a vehicle's arriving time $t{a_i}$ is not in the time window, a penalty ${w_i}$ will be added to the objective function. It also considers the discrete capacity feature of ADRs. For goods in depot, we use ${\text{|S|}}$ to denote the number of sizes of goods, where the physical size of goods decreases as the index $s$ increases. Considering the containers with different sizes, a vehicle owns ${f_t} \;t$-th-sized containers, where $t \in S$. A $t$-th-sized container can only store 1 $t$-th-sized good. We also allow cross-level storage, which means an $s - 1$-th-sized container is capable of being loaded with $n\;s$-th goods. In order to identify the size of goods, we denote by $cd_i^t$ the goods that need to be delivered in node $i \in I$ with the size $t \in S$, while $cp_i^t$ is for pick-up. In the delivery process, $md_i^t$ is the number of containers used for $t \in S$ size delivery container in node $i \in I$, and $mp_i^t$ is the number of containers used for picking up goods. We optimize the decision variable $X_{ij}^v$, which is 1 if vehicle $v \in V$ passes the link from node $i \in I$ to node $j \in I$, and 0 otherwise. As shown in Figure 1, different vehicle colors represent different tasks, while the colors of the nodes indicate various demand types: red for delivery, blue for pick up, and green for combined delivery and pick up. The numbers denote the sequence of vehicle routing, which is also the output of our model.

3.1.   Notation

Table 2 presents the notations in this model.

3.2.   Mathematical model

● P1: Delivery-Only Model

In this scenario, all goods are stored at the depot and must be delivered to various customer locations. This model is applicable in situations where a logistics center serves multiple customers with no return goods.

● P2: Pick-Up-Only Model

All goods are collected from customer locations and transported back to the logistics center. This model is applicable in situations where customer locations are primarily suppliers or the focus is on returning goods to a central location.

● P3: Combined Delivery and Pick-Up Model

This model combines the scenarios of P1 and P2, where goods are both delivered from the depot to customers and collected from customers back to the logistics center. It is a more complex situation where vehicles handle both deliveries and pick-ups simultaneously.

3.2.1.   Delivery model

The objective function (1) aims to minimize the weighted total distance, number of vehicles, and time window penalties.

Subject to:

Constraints (2) and (3) define the decision variables. Constraint (4) ensures that all vehicles leave from a fixed depot. Constraint (5) requires that all vehicles return to the depot. Constraint (6) indicates that the vehicles cannot travel to the same place. Constraint (7) calculates the number of activated vehicles. $\sum\limits_{v = 1}^{|V|} {X_{00}^v}$ represents the vehicles that were not activated, as they returned directly to the depot after leaving. Constraints (8) and (9) require that each node can only receive one service from one vehicle per application. It is important to note here that Constraints (8) and (9) are not incompatible with split delivery of discrete containers, as the dummy node technique will be introduced below, and this step will deal with discrete goods as multiple virtual nodes so that they can be delivered in a manner that satisfies Constraints (8) and (9) to satisfy a customer's split delivery. Constraint (10) ensures that vehicles must leave the same node after visiting it. Constraint (11) is used to restrict the time to reach node $j$. Constraint (12) defines the penalties incurred when a vehicle fails to meet the customer's requested time window. Constraint (13) ensures that the number of goods of different sizes loaded on the robot while starting from depot cannot exceed the capacity of the different size containers of the robot. ${n^{s - t}}$ represents how many $s - th$ goods the $t - th$ containers can accommodate. Thus, the left side of the inequality (13) calculates how many equivalent $s - th$ goods have been loaded onto vehicle $v$. The right side of the inequality (13) calculates the maximum amount of equivalent $s - th$ goods that the vehicle can load.

3.2.2.   Pick up model

The model in which all goods are requested to be delivered to the depot is formulated as follows:

Constraints (2)–(12), and

Constraint (15) ensures that the number of goods of different sizes loaded on the robot while coming back to depot cannot exceed the capacity of containers. The principle is consistent with Constraint (13).

3.2.3.   Delivery and pick up model

The model in which all goods are sent from or to the depot, usually expressed as One-to-Many-to-One Single Vehicle Pickup and Delivery Problem, is formulated as follows:

Constraints (2)–(12), and

If the delivery and pick up tasks are mixed, the complexity of the problem will increase dramatically. It is necessary to check whether all nodes meet the capacity constraint. Constraints (17) and (18) calculate the number of goods to delivery and the number of goods to be picked up before this node. Figure 2 describes the concrete explanation of Constraint (19), which determines whether the capacity constraints for different sizes are satisfied for each node, by calculating the quantity of goods that still need to be delivered and the goods received earlier, summing them up to obtain the current load.

3.3.   Revised model

3.3.1.   Model linearization

The non-linear Constraint (12) caused by calculating the time penalty function incurs a great computation challenge; as such, series of linearization methods are proposed as follows. Since this is a general linearization technique, we omit ($\omega$) in all decision variables for simplicity of notation. Constraint (12) can be divided into two parts. The question is transformed into two subproblems.

Constraint (21) is replaced by Constraints (23)–(28) to achieve model linearization.

Constraint (20) is first divided into two subproblems. When $\mu _i^e = 1$ and $\nu _i^e = 0$, Constraint (25) reduces to $w_{ci}^e \leqslant 0$ and Constraint (26) is redundant. Together with Constraint (23), we have $w_{ci}^e = 0$ and Constraint (24) reduces to ${t_{ei}} \leqslant {t_{ai}}$. This indicates that vehicle arrives at $i$ node later than the earliest time. The penalty for being early is equal to 0. When $\mu _i^e = 0$ and $\nu _i^e = 1$, Constraint (26) reduces to $w_{ci}^e \leqslant \alpha ({t_{ei}} - {t_{ai}})$ and Constraint (25) is redundant. Together with Constraints (23) and (24), we have $w_{ci}^e = \alpha ({t_{ei}} - {t_{ai}})$ and ${t_{ei}} \geqslant {t_{ai}}$. This represents that vehicle arrives at $i$ node earlier than the earliest time. The penalty for being early is equal to $\alpha ({t_{ei}} - {t_{ai}})$. When $\mu _i^e = 0$ and $\nu _i^e = 0$, we have $w_{ci}^e = \alpha ({t_{ei}} - {t_{ai}})$ and $w_{ci}^e = 0$. This represents the arrival of the vehicle just at the earliest allowed arrival time.

For Constraint (22), via the same linearization method, Constraint (22) can be substituted by Constraints (29)–(34).

3.3.2.   Method for calculating the minimum number of activated vehicles

The number of activated vehicles plays a vital role in objective function. During the solving process, solutions with a higher number of activated vehicles will be discarded due to the dominant cost of vehicle activation. To reduce the search space and simplify the objective function, this paper transforms the traditional optimization problem into the following bi-level optimization problem. Constraint (35) gives the upper model objective function.

The result $vn$ of optimizing the upper equation will be input into the lower model as the number of vehicles to optimize the lower equation.

The bi-level programming model is complicated to solve. However, in this paper, the number of startup vehicles can be uniquely determined in the upper model. This process is well-suited for P1 and P2. However, for P3, where certain cells can be used iteratively during the delivery process, we can only obtain a relatively approximate minimum number of vehicles. The method to calculate the minimum number of activated vehicles $vn$ is shown in Algorithm 1.

4.   Solution algorithm

The DCVRP-SDP-STW is a typical NP-hard problem. To address the large-scale delivery VRP challenge, we propose an improved ACO algorithm as depicted in Figure 3. Dummy nodes are introduced to achieve split deliveries, and a complete distance graph is generated to consider the 3D feature of buildings. As a key component of the improved ACO algorithm, the SMC algorithm clusters the nodes based on both time window distribution and spatial location, which reduces computing time. To improve the solution diversity of ACO, we put forward three mechanisms.

4.1.   Spatiotemporal multi-stage clustering

We propose a SMC method, integrating spatiotemporal dimensions for more effective classification of customer points. As illustrated in Figure 4, the SMC approach divides the problem into two stages. In the first stage, we focus on temporal information, specifically the time window of each customer. By employing Euclidean distance, initial clustering is performed to divide the delivery process into multiple periods. This process creates a distance matrix, as described in Eq (41), to cluster the time feature. The second stage involves clustering based on spatial features, such as $e$(enclosure), $f$(floor), and $d$(door) information, which describe the customer's spatial characteristics. The distance function between two nodes is computed according to Eq (42), ensuring that both spatial proximity and structural features are considered. The proposed spatiotemporal SMC method can reduce computing time by effectively decomposing large-scale routing problems into smaller, manageable sub-problems.

4.2.   TACSATP algorithm

This section introduces the details of the TACSATP algorithm. Algorithm 2 gives the pseudo-code. The TACSATP algorithm is an ACO-based metaheuristic algorithm. It employs multiple ants to iteratively construct solutions by dynamically managing capacity, cluster constraints, and node selections via a tabu table, while guiding search efficiency through spatiotemporal pheromone updates. The three subsections introduce three mechanisms used in this algorithm. To address the bi-level optimization problem, we propose a hierarchical solution update strategy that sequentially evaluates: 1) the number of activated vehicles as the primary criterion, and 2) the objective function value defined in Eq (36) as the secondary criterion.

4.2.1.   Ant colony optimization with tabu table

In this section, we discuss the combination of ACO and TS. To alleviate the issue of ACO easily falling into local optima, we introduce a tabu table during the search process to maintain solution diversity by storing previously visited solutions. Unlike strict prohibition mechanisms, the tabu table in TACSATP implements a soft penalty mechanism, allowing some flexibility in revisiting solutions under certain conditions. Algorithm 3 illustrates the process of searching methods with tabu tables.

In the ACO process, the next node is selected using a roulette wheel method to introduce randomness, ensuring diverse path exploration. After selecting the next node, the algorithm checks the tabu table: if the node is new, its corresponding edge is initialized with an iteration countdown before release; if it already exists, the countdown is decremented by one. If the node remains restricted, the selection repeats until a valid node is found. This approach balances solution diversity and feasibility, prevents premature convergence by avoiding repetitive path revisits, and enhances overall search efficiency.

4.2.2.   Improved search mechanism

In the TACSATP algorithm, two mechanisms are designed to generate routes.

Mechanism 1: It adapts the search process to the clustering method by changing the range of the candidate set. During route generation, a node is randomly selected from multiple clusters. Nodes within the same cluster are then selected until the cluster is exhausted, a process called intra-cluster node selection. After that, the process returns to global node selection to determine the next node from a different cluster. This approach allows the ACO to focus on node traversal within clusters, significantly reducing computational complexity.

Mechanism 2: To address the issue of potential capacity waste with discrete capacity containers, this mechanism introduces a lookahead feature. When a delivery robot's large capacity container cannot be fully utilized, the robot would return to the depot, leaving the small capacity container underutilized. The lookahead mechanism allows the search to consider multiple steps forward if capacity constraints are not met, thus improving the overall utilization rate of the robot's containers.

In Figure 5, the relationship between two mechanisms is reflected in how they interact during the search process: Mechanism 1 affects the route construction by structuring the node selection order, whereas Mechanism 2 affects decision-making by improving capacity utilization.

4.2.3.   Spatiotemporal pheromones based on prior knowledge

In the process of constructing the route selection, we propose the search method of long- and short-time pheromone. The specific search method is shown in Eq (43), where the time pheromone is calculated as Eq (44). In this way, the time and distance penalty relationship can be jointly considered during the search process, allowing for the capture of spatiotemporal relationships.

To consider the local information and a priori information in the search process, our approach is divided into two parts. The first part captures points with the current time window penalty, referred to as short-time pheromones, shown in Eq (45). $z_{ij}^l$ represents the time from i node to j node for the l ant. The second part introduces a priori knowledge from previous iterations to capture long-term time penalty relationships named global time pheromone. The global time pheromone computation process is illustrated in Figure 6. In order to give ambiguity to vehicle arrivals to improve the performance of the solution, this paper constructs the possibility of vehicle arrivals in Eq (46).

After deriving the arrival density function, integrating over the range of $t{e_j}$ and $t{l_j}$, as in Eq (48), the probability of the vehicle arriving on time can be computed, denoted as $\Delta {\omega _j}$. Inspired by the pheromone mechanism in the ACO algorithm, an iterative process also exists for the time pheromone:

In this way, a memorizing and forgetting process for global time information is created.

5.   Experiment study

5.1.   Experiment setting

To demonstrate the proposed models and algorithms, the express order information and map of Southeast University, Nanjing, China is used. The selected building information includes four buildings, A, B, C, and D, each with six floors and ten rooms per floor, with a total of 240 rooms. The specific map information of the park can be seen in Figure 7.

The generated order information and time distribution are shown in Table 3 and Figure 8.

5.2.   Data preprocessing

5.2.1.   Creating distance graph

As shown as Figure 9, the Euclidean distance and Manhattan distance are not applicable in a 3D environment, because ADRs can move between floors including climbing stairs or using an elevator. In this paper, the following algorithm is proposed, which assumes the complete graph of distances is $G = (V, E)$, where $V = \{ {v_1}, {v_2}, ..., {v_n}\}$ corresponds to customer nodes $I = \{ {i_1}, {i_2}, ..., {i_n}\}$, with $n$ being the number of customers. Each node records its corresponding 3D spatial information $({e_i}, {f_i}, {d_i})$, where ${e_i}$ is the coordinate, ${f_i}$ indicates the floor number, and ${d_i}$ represents the corresponding door number. The purpose of the algorithm proposed in this paper is to calculate the weight matrix $W$.

It can be observed that in distance calculations in a 3D environment, three primary distinct cases exist: ① From depot to customer, ② From customer to customer between different buildings, and ③ From customer to customer within the same building. Figure 9a) illustrates the movement of the delivery vehicle between buildings for these three cases. Figure 9b) provides a more intuitive representation by placing this scenario within a 3D environment. By applying Algorithm 4, the corresponding complete distance graph can be obtained in Figure 9c).

5.2.2.   Adding dummy node

Due to the discrete capacity window feature, dummy nodes are introduced. This paper generates a distance complete graph in 3D environment $G = (V, E, W)$, where vertex set $V$ is a combination of the customer set ${\text{V }} = {\text{ }}\left\{ {{v_1}, {\text{ }}{v_2}, ...., {v_n}} \right\}$ and depot ${v_0}$. For each customer node ${v_i}$, if ${v_i}$ owns ${m_i}$ demand orders, the customer node ${v_i}$ is then augmented to $\{ v_i^1, v_i^2, ...., v_i^{{m_i}}\}$. $\{ v_i^1, v_i^2, ...., v_i^{{m_i}}\}$ corresponding to a 3D location information remaining consistent with the 3D spatial information of the node ${v_i}$ in $G$. After that, the corresponding generated graph $G' = (V', E')$ is input into the weight calculation algorithm, resulting in a complete graph $G' = (V', E', W')$, where the vertex set ${\text{V' }} = {\text{ }}\left\{ {v_1^1, {\text{ }}v_1^2, ...., v_n^{{m_n}}} \right\}$. The number of vertexes $V'$ is $\sum\nolimits_{i = 1}^n {{m_i}}$. Thus, after generating the augmented complete graph of distances, the distance between any two nodes can be directly indexed during the algorithm implementation.

5.3.   Optimization results

In this paper, we set $|S| = 2$, as the discrete containers are divided into two types: large and small, corresponding to large goods and small goods; the large container can accommodate $n = 4$ small containers. We generate customer demands for different models with different customer numbers under ratios of large to small goods (Ratio = sum of large goods number/ sum of small goods number). For example, if the total number of customer demands is set to 100 with a large-to-small ratio of 1:1, we randomly assign 50 large goods and 50 small goods across the customers. The corresponding augmented complete distance graph's cumulative order numbers range from 200 to 1500. We also give the baseline. The characteristic of baseline route planning is that vehicles always first serve customers near the depot such as in the order of T5A101, T5A102, T5A103, etc., resulting in a lower distance cost. In conducting the experiments, this paper uses the following parameter settings as Table 4.

The complete 45 optimization results can be viewed in Table 6 in the appendix. Results show that the algorithm achieves large reductions in total costs across all scenarios compared to the baseline, with improvement ratios ranging from 0.23 to 0.94. The algorithm adeptly balances distance and time costs to optimize the delivery process. And the number of activated vehicles can often reach the minimum. Regarding P3, it adjusts constraints to integrate delivery and pick-up operations, resulting in a significant reduction in activated vehicles. It indicates the effectiveness of the search mechanism.

To highlight the heterogeneity among models, Figure 10 compares their performance under a scenario with 200 customers and a goods size ratio of 1:2. P1 and P2 show similar results across all metrics, while P3 stands out with significantly higher total cost and computation time, but fewer required vehicles—reflecting its use of a lookahead search mechanism. Cost breakdown shows that distance costs are relatively consistent across models, whereas time costs dominate, especially in large-scale instances. Notably, P3 achieves the lowest distance cost due to more efficient container utilization—each container handles both delivery and pickup, reducing routing demands. The higher computation time for P3 further confirms its greater algorithmic complexity.

5.4.   Ablation experiment

We conduct the comprehensive analysis of several solution algorithms for the problem P3. The baseline is an unoptimized algorithm. The algorithms analyzed include standard ACO and several optimized variants, which are parts of the SMC-TACSATP ablation studies. All algorithms, except the baseline, use the SMC method to enhance computational speed, hence they are not specifically labeled.

● TACSATP: Fully optimized variant that uses all methods proposed in this paper

● ACSATP: A variant that only removes the tabu table

● ACSA: A variant that removes the temporal pheromones and the tabu list

● TACOTP: A variant that removes the searching mechanism proposed in this paper

● ACOTP: A variant that further removes both the searching mechanism and the tabu list

● ACO: Standard implementation of the ACO algorithm

● Baseline: Unoptimized routing generated based on the proximity principle

Table 5 shows the advantages of SMC-TACSAP by varying the size-to-goods ratio for a fixed number of 300 customers. The metrics compared include the objective function value (OV), the number of used vehicles (VN), the ratio of large to small goods (Ratio)—which represents the proportion of large goods relative to small goods—and the overall capacity utilization (AC). Its integration of time pheromones and a custom search mechanism significantly improves performance, outperforming all other variants in 11 out of 16 cases. Compared to the baseline and standard ACO, TACSATP reduces costs by over 80% in some settings, proving its strong efficiency and optimization capability in complex delivery environments.

In Figure 11, we fix the large-to-small goods ratio at 1:2. P3 experiments range from 50 to 500 customers, showcasing SMC-TACSATP's advantages. For comparative experiments, Baseline, ACSA, TACOTP, and TACSATP are selected. It can be observed that the indicators increase with the number of customers. To eliminate this increasing trend, Figure 12 calculates the mean values for each metric, including the average cost and the average number of activated vehicles per customer. By analyzing these two figures, it can be observed that TACSATP consistently achieves a lower value.

The proposed model and algorithms demonstrate substantial advantages in complex 3D delivery environments by accurately modeling multi-floor, multi-order scenarios through the use of virtual nodes and a customized distance calculation method. Specifically, our model effectively addresses the simultaneous pickup and delivery problem within a 3D environment, reducing overall operational costs while reliably meeting customer demands.

Key mechanisms such as the SMC and spatiotemporal pheromones informed by prior knowledge contribute directly to performance improvements. By effectively decomposing complex routing problems into manageable temporal and spatial sub-problems, the SMC approach optimizes computational efficiency and route effectiveness. Additionally, the proposed TACSATP algorithm leverages ACO integrated with a tabu table and improved search mechanisms, significantly enhancing the exploration of feasible solutions and preventing premature convergence.

Extensive experimentation demonstrates notable cost reductions across various customer scales and ratios of goods, achieving savings ranging from 23% to 94%. Particularly, the integrated P3 model utilizes a forward-looking search mechanism that effectively coordinates delivery and pickup tasks, thereby optimizing container utilization, minimizing vehicle deployment, and reducing travel distances.

6.   Conclusions

This paper presents a novel approach to solve the last-mile delivery problem in 3D environment by addressing the DCVRP-SDP-STW for electric ADRs. The improved ant colony search algorithm significantly enhanced solution quality. The case study indicates that the algorithm achieved optimization improvements ranging from 23% to 94% compared to the pre-optimization state. These results highlight the effectiveness of integrating advanced algorithms with 3D modeling in complex delivery systems.

This paper explores the large-scale route optimization within 3D environment. Future research might focus on refining operational models, particularly by integrating dynamic consumer behavior data and real-time adjustments in delivery strategies. Additionally, investigating the interaction between smart delivery robots and consumer preferences, as well as the potential for differentiated pricing mechanisms, could lead to more efficient and cost-effective logistics operations. By incorporating machine learning and advanced heuristic methods^[38,39], future studies can further enhance the quality of solutions in complex delivery systems, contributing to the development of more intelligent and autonomous logistics networks.

Use of AI tools declaration

The authors declare they have not used Artificial Intelligence (AI) tools in the creation of this paper.

Acknowledgments

This paper is supported by the Youth Program of National Natural Science Foundation of China (No. 72301066), Fundamental Research Funds for the Central Universities, China (No. 2242024RCB0048), Natural Science Foundation of Jiangsu Province (No. SBK2024021394), Wuxi Science and Technology Development Fund Project (No. K20231014), National Undergraduate Training Programs for Innovation and Entrepreneurship (NO. 202410286135Z).

Conflict of interest

The authors declare there is no conflict of interest.

Appendix