The Two-Echelon Unmanned Ground Vehicle Routing Problem: Extreme-Weather Goods Distribution as a Case Study

TL;DR

This paper introduces a new routing problem for unmanned ground vehicles in extreme weather and proposes a superior algorithm to solve it.

Contribution

A novel hybrid algorithm combining Artificial Bee Colony and Wild Horse Optimizer is proposed for the two-echelon UGV routing problem.

Findings

The HABC-WHO algorithm outperformed GA, DWHO, and DABC-FNS on 43 benchmark instances.

The algorithm showed strong solving capability and high precision in extreme-weather material distribution scenarios.

Strategies like large neighborhood search and 2-Opt improved the algorithm's performance.

Abstract

In extreme weather conditions, the use of unmanned ground vehicles (UGVs) for material distribution enhances safety. We introduce a two-echelon unmanned ground vehicle routing problem (2E-UGVRP) and proposes a hybrid Artificial Bee Colony–Wild Horse Optimizer (HABC-WHO) algorithm to solve it. In this approach, the optimal solution obtained from the Artificial Bee Colony algorithm replaces the worst solution of the Wild Horse Optimizer. To further improve the algorithm’s performance, strategies such as large neighborhood search, two-optimization (2-Opt) operation, and satellite subpath crossover are incorporated. The algorithm’s effectiveness is demonstrated through the solution of 43 benchmark instances, with performance comparisons against a Genetic Algorithm (GA), Discrete Wild Horse Optimizer (DWHO), and Discrete Artificial Bee Colony–Fixed Neighborhood Search (DABC-FNS). The results clearly show the significant superiority of the proposed algorithm. Additionally, the algorithm is applied to material distribution by two-echelon UGVs under extreme weather conditions, yielding promising results. Experimental findings indicate that the algorithm exhibits strong solving capability and high precision.

Figures

Full text

1. Introduction

With the continuous growth of e-commerce, urbanization leads to increased transportation requirements in cities, and the last-mile delivery (LMD) is the most expensive part of the supply chain, as high personnel costs incur [1]. Labor cost accounts for the majority of logistics costs. Thus, transportation companies do their best efforts to find a solution for cutting down dispatch costs. Then, the application of unmanned vehicles in logistics can efficiently cut down the dispatch costs and improve the efficiency of transportation.

Recently, several unmanned vehicles have emerged. Starship Technologies, which is an Estonian start-up company, develops autonomous robots for last-mile deliveries. Starship’s robots move along sidewalks and weigh no more than 40 pounds while fully loaded [2]. In 2016, Cainiao launched its first last-mile delivery robot, Xiao G. [3], and Domino’s unveiled its next-generation delivery robot, named DRU (Domino’s Robotic Unit). Dispatch launched its first autonomous delivery vehicle, Carry, which was designed to operate exclusively on bicycle lanes and sidewalks [3]. Piaggio Fast Forward (PFF) positions its products as ‘optimally sized transport solutions’, bridging the gap between drones and automobiles. One of its commercially developed models is Gita [3]. They can be applied to deliver parcels or groceries directly from warehouses or specialized hubs. After robots have finished their dispatch mission, they can autonomously return to their warehouse or hub. For safety, the robots are only permitted to move at pedestrian speed. For reducing the costs and the CO_2_ emissions, improving consumer satisfaction, and building an environmentally friendly economy, an unmanned ground vehicle is an excellent choice for logistics. COVID-19, the global pandemic, has made us realize the importance of contactless delivery services and robotic automation, which will further promote the widespread application of unmanned delivery technology.

Currently, there is extensive literature on parcel delivery using trucks and drones, while studies focusing solely on parcel delivery with unmanned ground vehicles (UGVs) are relatively scarce. The vehicle routing problem (VRP) is an NP-hard optimization problem that aims to determine a set of least-cost delivery routes from a depot to a set of geographically scattered customers, subject to side constraints [4].

In enclosed areas such as residential communities and campuses, UGVs can be used to reduce delivery costs and carbon emissions. To alleviate traffic congestion and enhance delivery efficiency, the concept of two-echelon delivery has been introduced [5]. The 2E-VRP is a variant of the VRP, involving a two-echelon distribution network with a central depot, a set of satellites, and a group of end customers. In the first echelon, the delivery task is typically handled by larger-capacity vehicles, while the second echelon involves smaller-capacity vehicles. Goods are first transported from the depot to the satellites by the first-echelon vehicles, then delivered to the end customers by second-echelon vehicles. Each customer is served by exactly one second-echelon vehicle. The objective is to minimize the total routing cost of both echelons. Figure 1 illustrates an example of the 2E-VRP.

Our research focuses on the practical application of two-echelon logistics distribution using UGVs under extreme weather conditions, which is a variant of the 2E-VRP. In this model, the first stage involves large-capacity UGVs delivering goods, while the second stage involves smaller-capacity UGVs for final deliveries. Due to the unsuitability of drones in extreme weather and to prevent risks to human life during the delivery process, UGVs are a more appropriate choice for such conditions. We refer to this model as the two-echelon unmanned ground vehicle routing problem (2E-UGVRP). The contribution and innovation of this paper is the development of a hybrid Artificial Bee Colony–Wild Horse Optimizer (HABC-WHO) to solve the problem. The Wild Horse Optimizer (WHO) [6], as a recently proposed metaheuristic, demonstrates strong local search capabilities when integrated with search strategies, while the Artificial Bee Colony (ABC) algorithm excels in global exploration. By replacing the worst solution in the WHO with the optimal solution from the ABC algorithm, the proposed hybrid algorithm enhances the WHO’s global search ability, mitigates the risk of premature convergence, and improves overall solution quality. Some algorithms are primarily suitable for small-scale problem instances and demonstrate limited effectiveness when tackling large-scale cases, such as linear programming methods. Certain algorithms, like Genetic Algorithms (GAs), are susceptible to premature convergence, while others, such as Particle Swarm Optimization (PSO), often exhibit slow convergence. The algorithm developed in this study effectively mitigates these shortcomings by accelerating convergence while maintaining solution diversity, thereby avoiding premature convergence and enhancing solution accuracy.

The remainder of the paper is organized as follows. Section 2 provides a non-exhaustive literature review. Section 3 presents a detailed problem description. Section 4 provides the mathematical formulation for the 2E-UGVRP. Section 5 explains the solution methodology developed to solve the problem. Section 6 presents extensive computational results. Section 6 and Section 7 present the computational tests and the application to real-world problems, respectively. Finally, Section 8 summarizes the work and highlights future research challenges.

2. Literature Review

The concept of the two-echelon vehicle routing problem (2E-VRP) was first introduced by Jacobsen and Madsen [5] in 1980 through a newspaper distribution system. Crainic, Ricciardi, and Storchi [7] later proposed the concept of a two-echelon delivery network for urban freight distribution, characterized by two levels of platforms. Gonzalez-Feliu et al. [8] were the first to introduce the two-echelon capacitated vehicle routing problem (2E-CVRP). Crainic, Ricciardi, and Storchi [9] refined the definition of the 2E-CVRP, describing it as a two-layer, time-dependent, synchronized, multi-route, multi-warehouse, multi-product, heterogeneous fleet problem with both hard (satellite) and soft (customer) time windows at each echelon. Sluijk et al. [10] explored the two-stage capacitated vehicle routing problem in urban logistics. Perboli, Tadei, and Vigo [11] introduced a mathematical model for the 2E-VRP with a single warehouse, naming it the 2E-VRP. In this problem, vehicles in the first stage transport goods from a warehouse to a satellite located at the echelon boundary. At the satellite, goods are unloaded from the first-stage vehicles and loaded onto smaller second-stage vehicles that meet the delivery requirements. Subsequently, the second-stage vehicles perform deliveries from the satellite to customers. Crainic et al. [12] designed a multi-heuristic algorithm that iteratively solved subpath problems in both stages while adjusting the satellite workload connecting them. In multi-echelon vehicle routing problems, deliveries from the warehouse to the customer are managed by rerouting and consolidating goods through intermediate satellites or distribution centers (cross-docks or distribution centers) [13].

Murray and Chu [14] were the first to study the truck–drone delivery problem, flying sidekick traveling salesman problem (FSTSP), where trucks serve as first-echelon delivery vehicles and drones serve as second-echelon delivery vehicles. Drones are launched from the truck at specified locations to deliver goods to customers. Agatz et al. [15] noted that trucks typically have longer ranges and can carry many parcels, while drones offer faster, cheaper, and more flexible deliveries. However, drone delivery is limited by payload and flight distance, which can be mitigated through the truck–drone delivery problem, improving delivery efficiency. When multiple trucks and drones are involved, the problem is referred to as the vehicle routing problem with drones (VRPD) [14]. Zhou et al. [16] proposed a new variant of the VRPD, where drones performed multiple back-and-forth trips while the paired vehicle stopped at customer nodes. Li et al. [17] introduced a model for homogeneous fleets of trucks and drones used to deliver parcels. In this model, each truck carries the same number of drones, and the characteristics of each truck and drone are identical. The problem is solved using the adaptive large neighborhood search (ALNS) heuristic algorithm.

The literature [18] also explores variants of the 2E-VRP, including the two-echelon stochastic vehicle routing problem with loading bay reservation strategy, which addresses parking and loading inefficiencies and is solved using a hybrid genetic algorithm. Zhong et al. [19] proposed an Artificial Bee Colony Genetic Algorithm (ABCGA) to solve standard 2E-VRPs with up to 51 customers. References [19,20,21,22] investigate the two-echelon vehicle routing problem for electric vehicles (2E-EVRP). Additionally, the two-echelon electric vehicle routing problem (2E-EVRP) has been studied, where charging or battery-swapping stations expand the transportation network. Crainic [23] and Marques [24] addressed the use of satellites as automatic charging stations, with secondary deliveries considered as multiple trips in the 2E-VRP. Vincent et al. [25] proposed a new 2EVRP variant incorporating time windows, intermediate facilities, and occasional drivers, termed the two-echelon vehicle routing problem with time windows, intermediate facilities, and occasional drivers (2E-VRPTW-IF-OD). A hybrid adaptive large neighborhood search (HALNS) algorithm was used for solving this variant. Lehmann and Winkenbach [26] developed a hybrid approach integrating exact methods with adaptive large neighborhood search for the time-constrained two-echelon routing problem. Dellaert et al. [27] introduced a variant of the 2E-VRP with customer hard time windows, resulting in synchronization issues between first- and second-echelon vehicles. Zamal et al. [28] developed an adaptive large neighborhood search algorithm to solve the two-echelon vehicle routing problem with pickups and deliveries. Lehmann and Winkenbach [26] designed a solution for the simultaneous pickup and delivery problem in a multi-trip 2E-VRP, while Mhamedi et al. [29] applied a branch-and-bound algorithm to solve a two-stage vehicle routing problem with time windows. They introduced effective deep dual-optimal inequalities and applied known valid inequalities. Dellaert N et al. [30] proposed a two-stage vehicle routing problem with time windows and multiple products. Zhou et al. [31] used a variable neighborhood tabu search algorithm for solving the 2E-VRP with time windows and simultaneous pickup and delivery. Kim et al. [32] explored the two-echelon hetero-collaborative routing problem (2E-HCRP), which facilitated parcel transfers between couriers and customer nodes, optimizing last-mile deliveries. Al Theeb et al. [33] proposed a heuristic solution based on greedy randomized adaptive search procedure to address the multi-objective two-echelon vaccine distribution problem. Grangier et al. [34] studied the two-echelon multiple-trip vehicle routing problem with satellite synchronization (2E-MTVRP-SS), solving it with an adaptive large neighborhood search method. Marques et al. [35] proposed an improved branch-price-and-cut algorithm to solve the two-stage vehicle routing problem with capacity constraints.

Anderluh et al. [36] considered a distribution system involving cargo bikes in the second echelon, with a real-world instance based on Vienna, Austria. Kergosien et al. [37] examined logistics problems at a hospital complex in Tours, France, where the first echelon involved vehicle transport between hospital units, and the second echelon focused on constructing delivery routes within individual buildings.

Currently, two main approaches are used to solve the 2E-VRP: exact algorithms such as Mixed-Integer Linear Programming (MILP) [2,14] and dynamic programming [38], and metaheuristic algorithms, including the Horsefly Algorithm [39], adaptive large neighborhood search (ALNS) [17,32,40], variable neighborhood search [31], tabu search [31,41], Genetic Algorithms [18,19,37], and Artificial Bee Colony (ABC) [19,42]. Currently, no algorithm outperforms all others across all problem instances. Table 1 provides an overview of the main solution methods for 2E-VRP.

The previous literature on two-echelon delivery problems has primarily considered truck–drone or truck–unmanned vehicle combinations. Our research, for the first time, introduces the two-echelon unmanned ground vehicle routing problem (2E-UGVRP), where unmanned ground vehicles (UGVs) are utilized in both delivery echelons. A novel solution approach, the hybrid Artificial Bee Colony–Wild Horse Optimizer (HABC-WHO), is proposed for this problem. This paper addresses a research gap and provides new insights into unmanned logistics optimization.

3. Problem Description

The two-echelon unmanned ground distribution network involves two fleets of vehicles that handle transfers at different stages [45]. In the first echelon, unmanned vehicles (first echelon) collect goods from the central depot and distribute them to several intermediate facilities known as satellites, where the goods are unloaded and transferred to another type of unmanned vehicle for further transportation to customers (second echelon). The central depot has a limited number of unmanned vehicles, each capable of servicing multiple satellite facilities. In the first echelon, if the delivery volume to a satellite exceeds the capacity of a single unmanned vehicle, multiple first-echelon vehicles may service that satellite. In the second echelon, since each customer’s demand is less than the capacity of the second-echelon vehicles, each customer is served by exactly one unmanned vehicle.

In the two-echelon unmanned ground vehicle routing problem, we define as arc (i,j) the direct route connecting node i to node j. If both nodes are satellites or if one is the depot and the other is a satellite, we define the arc as belonging to the first-echelon routing, while if both nodes are customers or if one is a satellite and the other is a customer, the arc belongs to the second-echelon routing. The first echelon connects the depots to the satellites; the second echelon is where the cargo is delivered from the satellites to the customers.

We denote the central depot as $[eqn]$ , the satellite facilities as $[eqn]$ , and the customers as C. The first-echelon fleet departs from the central depot, services the satellite facilities, and returns to the depot. The second-echelon fleet departs from the satellite facilities, serves the end customers, and returns to the satellite facilities. The capacity of each vehicle is fixed, with the first-echelon vehicles having a larger capacity than the second-echelon vehicles. The objective is to minimize the total transportation cost while meeting the vehicle capacity constraints and ensuring service to all customers. All customer demands are known, and since each customer’s demand is smaller than the capacity of the second-echelon vehicles, they are served by a single second-echelon vehicle.

In our model, we temporarily do not consider the fixed cost of vehicle transportation, assuming it to be constant. We focus on minimizing the total length of the transportation routes. Therefore, the total transportation cost consists of two components: the cost from the central depot to the satellites and the cost from the satellites to the end customers.

4. Mathematical Formulation

Let an undirected graph $[eqn]$ be given, where the node set N is partitioned as $[eqn]$ . $[eqn]$ denotes the central depot, $[eqn]$ represents the set of $[eqn]$ satellites, and $[eqn]$ denotes the set of $[eqn]$ customers. The edge set E is defined as $[eqn]$ .

Let $[eqn]$ and $[eqn]$ denote the number of vehicles available in the first and second echelons, respectively. The parameter $[eqn]$ represents the maximum number of second-echelon delivery routes that can be initiated from satellite k. Vehicle capacities in the first and second echelons are represented by $[eqn]$ and $[eqn]$ , respectively. The demand of each customer $[eqn]$ is denoted by $[eqn]$ , and $[eqn]$ is the cost associated with traversing arc $[eqn]$ .

Each satellite $[eqn]$ is associated with a per-unit handling cost $[eqn]$ , which includes loading and unloading operations. The binary variable $[eqn]$ equals one if satellite k is activated in the distribution network, and zero otherwise. The variables $[eqn]$ and $[eqn]$ represent the volume of goods transported through arc $[eqn]$ in the first and second echelons, respectively. Let $[eqn]$ denote the number of unmanned ground vehicles traversing arc $[eqn]$ in the first-echelon network. The binary variable $[eqn]$ equals one if arc $[eqn]$ is included in a second-echelon routing path originating from satellite $[eqn]$ , and zero otherwise. The binary variable $[eqn]$ equals one if customer $[eqn]$ is assigned to satellite k, and zero otherwise.

The model to minimize the total cost of the system may be formulated as follows:

$$[eqn]$$

Let $[eqn]$ represent the fixed cost for each satellite’s participation in the operation. To simplify the model, we omit the fixed cost, and thus Function (1) can be simplified to (2):

$$[eqn]$$

subject to

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Constraint (3) limits the number of vehicles used by the first-echelon fleet to the total number of vehicles available at the central depot. Constraint (4) specifies that if k $[eqn]$ , each first-echelon route must begin and end at the depot, meaning the unmanned vehicles in the first-echelon fleet start their route at the central depot and return there upon completing deliveries. Constraint (5) restricts the number of second-echelon unmanned vehicles to the total available in the second-echelon fleet. Constraint (6) imposes capacity restrictions on the satellite depots. Constraint (7) ensures vehicle flow balance at the satellite depots, maintaining consistency in the number of second-echelon vehicles entering and leaving each satellite. Constraint (8) requires that the total volume dispatched from the central depot equals the total demand of all customers. Constraint (9) enforces volume balance at the satellite nodes, ensuring that the inbound volume matches the outbound volume, and prevents the formation of sub-tours. Constraint (10) ensures that the total delivery volume on the first-echelon routes does not exceed the cumulative capacity of all first-echelon unmanned vehicles. Constraint (11) similarly ensures that the total delivery volume on the second-echelon routes does not surpass the total capacity of the second-echelon unmanned vehicles. Constraints (12) and (13) stipulate that routes must be completed before the vehicles return to their respective depots. Constraints (14) and (15) ensure that customer j is assigned to satellite k only if satellite k is responsible for delivering goods to customer j, thereby allocating the customer to the service area of the corresponding satellite. Constraints (16) and (17) ensure that each customer is serviced by exactly one second-echelon unmanned vehicle. Constraint (18) mandates that each customer is assigned to only one satellite. Constraint (19) establishes that a second-echelon unmanned vehicle can only depart from a satellite after the first-echelon vehicle has delivered goods to that satellite. Constraints (20) through (24) allow multiple first-echelon unmanned vehicles to service a single satellite depot if necessary.

5. Solution Method

We first assign customers to satellites and propose a hybrid Artificial Bee Colony–Wild Horse Optimizer (HABC-WHO) using K-means clustering to solve the 2E-UGVRP.

5.1. Construction of Initial Solutions

The initial solutions are constructed using a clustering approach. Specifically, all customers are grouped using the K-means clustering algorithm, and each customer is assigned to the nearest satellite. The primary paths are then generated through a random permutation to form the initial solution.

5.2. The Principle of Wild Horse Optimizer (WHO)

The Wild Horse Optimizer (WHO) [6] is a novel metaheuristic algorithm proposed by Iraj Naruei and Farshid Keynia in 2021, designed to solve optimization problems in continuous systems. The algorithm is inspired by the foraging, mating, and group leadership behaviors observed in wild horses. A wild horse herd consists of stable family groups, each led by a stallion (the leader), along with one or more mares and their offspring. The stallion leads the group members in search of suitable habitats. The basic principles of the Wild Horse Optimizer include the following five components:

(a) Creation of initial population, grouping of herds, and selection of leaders: An initial population is randomly generated and grouped. Each group is led by a stallion, with the remaining members consisting of mares and offspring [6].

(b) Foraging and mating behaviors of wild horses: Foals spend most of their time foraging within the group. In the foraging process, the stallion is considered the center of the foraging area, and group members move and explore around the leader within varying radii. A foal may leave group i and join a temporary group, while another foal may leave group j and join the same temporary group. If these two foals are of opposite sex and not closely related, they can mate. The offspring produced from this mating must leave the temporary group and join another group (e.g., group k), forming the next generation. This process of leaving, mating, and breeding is repeated by all foals within the herd [6].

(c) Leadership of the group by the stallion: The stallion leads its group members to more favorable habitats. If the current group occupies a dominant position in a given area, they continue to stay there. However, if another group holds a more dominant position in the same area, the current group must relocate. The habitat represents the current optimal solution [6].

(d) Communication and selection of leaders: Leaders are selected based on fitness. If a group member has a better fitness value than the current leader, the leader’s position is exchanged with the member. The leader’s position represents the best solution for the group [6].

(e) Saving the optimal solution: The fitness values of the leaders of all groups are compared, and the group leader with the best fitness value represents the global optimal solution. This solution is then stored [6].

5.3. Solution Encoding and Decoding

The paths for the first and second stages are encoded separately. In the first stage, the number of satellites is limited, so we use discrete encoding; in the second stage, we use continuous encoding. Suppose the first stage consists of 3 satellite stations, with the encoding format shown in Figure 2. The second-stage path represents the corresponding customers, which are decoded into integer paths using the Min-position Matched Value (MPMV) method [46]. The first- and second-stage paths are separated by ‘0’, where the numbers before and after ‘0’ represent different vehicles delivering goods. The path in Figure 2a, $[eqn]$ , indicates that two primary unmanned vehicles are involved in the delivery. The first vehicle delivers to satellite stations 2 and 1, while the second vehicle delivers to satellite stations 1 and 3. Figure 2b shows the customer codes corresponding to the satellite deliveries, where the customers belonging to S1 are served by two vehicles.

The WHO is an optimization algorithm designed for continuous problems, and the solution it represents is continuous. Therefore, to solve the 2E-UGVRP, the solution represented by WHO needs to be discretized and decoded to represent the customer delivery path. The decoding method used here is the Min-position Matched Value (MPMV) [46] method. After decoding the scheme in Figure 2 using MPMV, the corresponding delivery scheme is shown in Figure 3.

5.4. Wild Horse State Updates

The positions of stallions and foals are defined as $[eqn]$ , $[eqn]$ , where $[eqn]$ represents the ith individual, and $[eqn]$ represents the population size; n is the number of customers.

The foal grazing behavior is carried out according to (25):

$$[eqn]$$

The $[eqn]$ is the leader position, the optimal individual in this group. $[eqn]$ is the current position of the group member. R is a random number within [−2, 2], which mainly controls the angle between individuals and the leader. The adaptive mechanism Z is calculated according to (26)–(28). The symbol ‘×’ denotes multiplication.

$$[eqn]$$ $$[eqn]$$ $$[eqn]$$ $$[eqn]$$

Let P be a binary vector composed of 0 and 1, both $[eqn]$ and $[eqn]$ are random vectors uniformly distributed within the range [0, 1]; $[eqn]$ is a scalar random number in the range [0, 1]. The variable $[eqn]$ denotes the set of indices in the random vector $[eqn]$ that satisfies $[eqn]$ . The symbol $[eqn]$ represents the dot product operation. $[eqn]$ is an adaptive mechanism that defines a coefficient linearly decreasing from 1 to 0, with $[eqn]$ denoting the current iteration number.

A random number Rand is generated. If Rand exceeds the crossover probability $[eqn]$ , the position of the foal is updated according to Equation (26); otherwise, the update is performed based on Equation (30).

$$[eqn]$$

The $[eqn]$ represents the position of individual q in group i after it returns to group i after outliers, and similarly, the $[eqn]$ represents the position of individual z in group j after it returns to group j. $[eqn]$ also represents the offspring individual produced by the mating of individual q in group i and individual z in group j. The positions of the two parents are represented inside the brackets in (30). The end denotes the last foal in the group. Crossover means that the qth foal from group i mates with the zth foal from group j, and the resulting offspring becomes the pth horse in group k. The value of the pth foal is obtained by calculating the average of the qth and zth foal.

The leader position is updated as described in Equation (31).

$$[eqn]$$

where $[eqn]$ denotes the next position of the leader in group i, while $[eqn]$ represents the current global best solution. $[eqn]$ is the current position of the group i leader, and Z is an adaptive parameter calculated by (32). R is a random number within [−2, 2], and $[eqn]$ = 3.14.

The fitness function is represented by (32), where fitness is the individual’s fitness, and $[eqn]$ represents the path length:

$$[eqn]$$

The pseudo-code of WHO is shown in Algorithm 1. Algorithm 1 Pseudo-code of WHORequire: Objective function $[eqn]$ Ensure: The best solution $[eqn]$ 1:Initialize the wild horse population( $[eqn]$ ), $[eqn]$ //Initialize the population.2: $[eqn]$ =13:Create groups and select leaders 4:Find the best horse as the solution // Find the global optimal value 5:while  $[eqn]$  do6:       Calculate $[eqn]$ by Equation (32) // Calculate the population fitness 7:       Calculate $[eqn]$ by Equation (29) 8:       for  $[eqn]$  do 9:          Update the position of the stallion by Equation (31) 10:        for  $[eqn]$  do 11:              Update the position of the foal by Equation (25) 12:        end for 13:        Select leaders // Find the optimal values for each group 14:        Find the best leader // Find the global optimal value 15:       end for 16:end while17:Return global best solution $[eqn]$

5.5. Hybrid Artificial Bee Colony–Wild Horse Optimizer (HABC-WHO)

To improve the accuracy of the solution, we propose a hybrid optimization algorithm named the Hybrid Artificial Bee Colony–Wild Horse Optimizer (HABC-WHO). This hybrid approach is designed to leverage the complementary strengths of both algorithms. The core idea of HABC-WHO is to enhance the global search capability of the Wild Horse Optimizer (WHO) by integrating the exploitation ability of the Artificial Bee Colony (ABC) algorithm [47]. Specifically, during the iteration process, the worst individual in the WHO population is replaced with the best solution obtained from the ABC component. This mechanism allows the algorithm to avoid premature convergence and maintain population diversity, thereby improving its convergence speed and solution accuracy.

The Artificial Bee Colony (ABC) algorithm, originally proposed by Karaboga from Erciyes University in Turkey [47], is a swarm-based optimization technique inspired by the intelligent foraging behavior of honeybee colonies. By embedding the ABC’s intelligent local search into the WHO framework, the proposed HABC-WHO aims to strike a better balance between exploration and exploitation, making it more effective in solving complex optimization problems.

The HABC-WHO algorithm procedure is shown in Algorithm 2.

The initial solutions are constructed using a clustering approach. All customers are clustered using k-means clustering, and assigned to the nearest satellite. The primary path is initialized by randomly arranging the customers to form the initial route. Algorithm 2 Pseudo-code of HABC-WHORequire: Objective function $[eqn]$ Ensure: The best solution $[eqn]$ 1:Initialize the wild horse population( $[eqn]$ ), $[eqn]$ //Initialize the population 2:Iter=1 3:Create groups and select leaders 4:Find the best horse as the solution // Find the global optimal value 5:while  $[eqn]$  do6:       Calculate $[eqn]$ by Equation (32) // Calculate the population fitness 7:       Calculate TDR by Equation (29) 8:       for  $[eqn]$  do 9:          Update the position of stallions by Equation (31) 10:          Large neighborhood operator //Apply large neighborhood operations for stallions 11:        for  $[eqn]$  do 12:              Update the position of foals by Equation (25) 13:              Large neighborhood operator //Apply large neighborhood operations for stallions 14:        end for 15:        Find the $[eqn]$ of each group 16:        if $[eqn]$  then 17:              Exchange $[eqn]$ , $[eqn]$ //The best foal in the group is exchanged with the stallion 18:        end if 19:        Find the best leader // Find the global optimal value 20:         $[eqn]$ decoding 21:       end for 22:       Assign the wild horse population to the ABC population 23:       Perform large neighborhood operation during the employed bee phase. 24:       Update the bee colony information. 25:       The best individual of ABC replaces the worst individual of WHO 26:       2-Opt operator 27:       Satellite subpath crossover operator 28:       if $[eqn]$ < $[eqn]$  then 29:        Exchange $[eqn]$ , $[eqn]$ 30:       end if 31:end while32: Return global best solution $[eqn]$

5.6. Search Strategy

To improve the convergence speed of the algorithm, we introduce different search strategies, including the large neighborhood search strategy, 2-opt strategy [48], and satellite subpath crossover strategy.

5.6.1. Large Neighborhood Search

The basic Wild Horse Optimizer has a relatively slow convergence rate. To accelerate convergence, we introduce a large neighborhood search strategy [36], which includes destroy and repair operators. The destroy operator randomly removes several customers from a route and places them into a customer pool. The repair operator then reinserts the removed customers back into the satellite routes using a roulette wheel selection mechanism. This large neighborhood search strategy is applied separately to the routes in both the first and second stages.

5.6.2. Two-Optimization (2-Opt) Operation

To enhance local search capabilities, we incorporate the two-optimization (2-Opt) [48] search strategy. This strategy is applied to each satellite route by disconnecting two edges within a route and reconnecting them in a different configuration, as illustrated in Figure 4. If the new configuration yields a lower cost than the original, the new connection is adopted.

5.6.3. Satellite Subpath Crossover Strategy

The subpaths between two satellites are crossed. A point is selected on each satellite’s subpath, dividing each subpath into two segments. Then, the two segments of each satellite’s subpath are exchanged. The satellite path crossover operation is shown in Figure 5. Figure 5a shows the paths of two satellites before crossover, and Figure 5b shows the paths of the two satellites after crossover.

6. Computational Tests

6.1. Parameter Settings for HABC-WHO

When the population size was set to $[eqn]$ = 50, the stallion proportion ( $[eqn]$ ) could be selected as either 0.2 or 0.5 to ensure an equal number of colts in each group. When $[eqn]$ , the number of foals in each group exceeded that of stallions, thereby promoting diversity among the colts within the group. Therefore, $[eqn]$ was set to 0.2 in these experiments. Extensive experimental evaluations revealed that the algorithm achieved optimal performance when the mating probability (PC) was set to around 0.13. Thus, $[eqn]$ was chosen as 0.13 in the experiments.

The following were the relevant parameter settings for HABC-DWHO:

• Population size: popN = 50.
• Stallion ratio: PS = 0.2.
• Mating probability: PC = 0.13.
• Maximum number of iterations: MI = 300.
• Number of stallions: NS = popN × PS = 10.
• Number of foals: Nf = popN × (1 − PS) = 40.

The computational configuration was as follows:

• OS: Windows 10 (×64).
• CPU: Intel Core i5-11400 (2.60 GHz).
• RAM: 16 GB.
• Language: Matlab 2016B.

The parameter settings for the three comparison algorithms were as follows: For the GA, the population size was set to 50, with a crossover probability of 0.8 and a mutation probability of 0.2. The maximum number of iterations was set to 300. For the DWHO, the population size was also set to 50, with a stallion ratio of 0.2 and a mating probability of 0.13, while the maximum number of iterations was kept at 300. For the DABC-FNS, the population size was set to 50, the limit parameter was set to 50, and the maximum number of iterations was set to 300.

6.2. Experimental Verification

Here, we used standard instances to verify the results of our algorithm. The results were compared with those obtained by a Genetic Algorithm (GA), a Discrete Wild Horse Optimizer (DWHO) [49], and a Discrete Artificial Bee Colony–Fixed Neighborhood Search (DABC-FNS) [42], as shown in Table 2, Table 3 and Table 4.

The Genetic Algorithm (GA) follows an evolutionary framework that balances exploration and exploitation strategies. Exploration is facilitated through crossover operators such as order crossover (OX) and partially mapped crossover (PMX), while exploitation is enhanced via mutation and local search techniques, including the 2-Opt heuristic and large neighborhood search (LNS). The DWHO is derived from the Wild Horse Optimizer (WHO) and incorporates three local search strategies—swap, reverse, and insertion operations—along with the largest-order-value (LOV) decoding technique [49]. The DABC-FNS framework is based on the discrete Artificial Bee Colony (ABC) algorithm with fixed neighborhood search strategies, including random single-point insertion, random sequence insertion, and random insertion. Moreover, a 2-opt strategy is applied within a fixed neighborhood to further refine local search performance [42].

The proposed algorithm was compared with the GA, DWHO, and DABC-FNS in solving the standard benchmark problems from sets 2, 3, and 4. As shown in Table 2, for the set 2 instances, the proposed algorithm significantly outperformed the GA and DWHO. In comparison with the DABC-FNS, the proposed method was superior in 10 instances, equivalent in 5 instances, and slightly worse in only 2 instances. In Table 3, the proposed algorithm again demonstrated superior performance compared to the GA and DWHO, with better results in 11 instances, equivalent performance in 4 instances, and slightly inferior results in 3 instances when compared to DABC-FNS. Furthermore, as indicated in Table 4, the proposed algorithm outperformed the other three algorithms in solving the set 4 instances, showcasing its superior solution capability. In Table 2 and Table 3, the negative sign indicates that the obtained solution was better than the best known solution.

To further verify the effectiveness of the algorithm proposed in this paper, we used the Wilcoxon’s rank sum test, a nonparametric statistic test, and the significance level $[eqn]$ was set to 0.05. The symbols ‘+’, ‘=’, and ‘−’ denote that HABC-WHO was significantly better than, similar to, and significantly worse than its competitor on an instance, respectively, according to the p value. The comparison results are shown in Table 5, Table 6 and Table 7. It can be seen from Table 5, Table 6 and Table 7 that the solving ability of the algorithm proposed in this paper was significantly superior to that of the three algorithms being compared. From the convergence curves in Figure 6, it can be similarly observed that the proposed algorithm demonstrates superior solving capability compared to the three benchmark algorithms.

6.3. Time Complexity Analysis

The maximum number of iterations for the proposed algorithm is $[eqn]$ , and the population size of the HABC-WHO algorithm is m. The number of customers is n. The time complexity of population initialization is $[eqn]$ , while both the herding and mating operations have a time complexity of $[eqn]$ . The fitness evaluation requires $[eqn]$ , the large neighborhood search operation has a complexity of $[eqn]$ , the subpath crossover operation is O(n/2), and the computation of the global best solution is O(n). Additionally, the Artificial Bee Colony (ABC) algorithm has a time complexity of $[eqn]$ . Therefore, the overall time complexity of the proposed algorithm is $[eqn]$ , which simplifies to $[eqn]$ .

7. Application to Real-World Problem

In a university campus scenario under extreme weather conditions, it is essential to ensure continuous supply of daily necessities to students. Due to the extreme weather, traditional delivery methods and drone delivery are not feasible, so unmanned ground vehicles (UGVs) are employed. The case involved one distribution center, two satellites, and fifty terminal receiving users. The geographic distribution of the case is shown in Figure 7. For computational convenience, the Euclidean distance was used to measure distances between points. The HABC-WHO achieves an optimal solution of 51.85, demonstrating significant improvement over competing methods. Table 8 summarizes the comparative performance of the four algorithms on this instance. The results demonstrate that the proposed algorithm consistently outperforms the three competing methods across all comparison metrics.

8. Conclusions

Previous studies on two-echelon delivery problems primarily considered truck–drone or truck–unmanned vehicle combinations. This study, for the first time, introduced the two-echelon unmanned ground vehicle routing problem (2E-UGVRP), where unmanned ground vehicles (UGVs) were utilized in both delivery echelons. A corresponding mathematical model was formulated, and a novel solution approach, the hybrid Artificial Bee Colony–Wild Horse Optimizer (HABC-WHO), was proposed. This work addressed a research gap and provided new insights into unmanned logistics optimization. The HABC-WHO proposed in this paper demonstrated enhanced capability in solving the 2E-UGVRP problem. By incorporating strategies like large neighborhood search, 2-Opt local search, and satellite subpath crossover, the algorithm was shown to outperform others in solving 43 benchmark instances. Additionally, the application to a case study under extreme weather conditions revealed that the proposed algorithm’s global search strategy was significantly better than that of the compared algorithms, highlighting the effectiveness and superiority of the approach. Future research may focus on multi-objective extensions, such as incorporating time window constraints and customer satisfaction levels. In future work, we plan to explore the potential of this algorithm for solving other discrete combinatorial optimization problems, further expanding its applicability.