The rapid development of e-commerce has increasingly exposed the limitations of traditional manual delivery in solving the “last-mile” problem. In urban environments, Unmanned Aerial Vehicles (UAV drones) offer a promising alternative, operating free from ground-level congestion and traffic restrictions, thereby alleviating the pressure on couriers. A prerequisite for efficient delivery is a rational, optimized task allocation scheme that maximizes operational efficiency before the drones are dispatched.
The problem of multi-UAV task allocation for urban logistics is a combinatorial optimization challenge subject to multiple objectives and constraints. Constructing a suitable mathematical model is the foundational step in addressing this problem. This paper proposes a balanced multi-UAV task allocation model optimized by a Multi-strategy Improved Hippopotamus Optimization (MIHO) algorithm. The model defines economic cost, delay penalty, and safety risk incorporating load balance as its optimization objectives. To solve this model effectively, we enhance the original Hippopotamus Optimization (HO) algorithm by employing Latinized Stratified Sampling for population initialization, integrating a Sine-Cosine strategy to bolster global search capabilities, and fusing a Simulated Annealing mechanism to refine solution precision and avoid premature convergence. Simulation results validate that our approach effectively meets urban logistics demands, with the MIHO algorithm demonstrating superior performance in convergence accuracy and stability compared to several state-of-the-art algorithms.
1. Problem Description and Model Formulation
Consider an urban logistics operation where a single distribution center (depot) serves a set of customer locations. A fleet of heterogeneous, rechargeable multi-rotor UAV drones is available at the depot. Each customer point has a known location and delivery demand. Every UAV drone takes off from the depot, sequentially serves its assigned customer points, and must return to the same depot upon completion. The core problem is to assign customers (tasks) to UAV drones and determine their service sequence to optimize the defined objectives while adhering to operational constraints.
1.1 Decision Variables
Let $N$ be the number of available UAV drones, denoted as $U = \{U_1, U_2, …, U_N\}$. Let $M$ be the number of delivery tasks (customer points), denoted as $T = \{T_1, T_2, …, T_M\}$. We define a virtual return task $T_0$ (or $T_{M+1}$) representing the return to the depot. The binary decision variable $x_{ij}$ is defined as:
$$
x_{ij} =
\begin{cases}
1, & \text{if task } T_j \text{ is assigned to UAV } U_i \\
0, & \text{otherwise}
\end{cases}
\qquad (1)
$$
1.2 Objective Functions
The total objective function $F$ is a weighted sum of three normalized sub-objectives: Economic Cost ($C$), Delay Penalty ($S$), and Safety Risk ($K$).
$$
\min F = \alpha_1 C + \alpha_2 S + \alpha_3 K \qquad (2)
$$
where $\alpha_1, \alpha_2, \alpha_3$ are weight coefficients satisfying $\alpha_1 + \alpha_2 + \alpha_3 = 1$. Each sub-objective is detailed below.
1.2.1 Economic Cost ($C$)
The economic cost comprises delivery (energy) cost $C_1$, management cost $C_2$, and fixed operational cost $C_3$.
Delivery Cost: Accounts for energy consumption during loaded and unloaded (return) flight.
$$
C_1 = \sum_{i=1}^{N} \sum_{j=1}^{M+1} (\lambda_{ij} D_{ij} + \mu_{ij} D_{i0}) \cdot x_{ij} \qquad (3)
$$
where $\lambda_{ij}$ and $\mu_{ij}$ are the energy cost per unit distance for loaded and unloaded flight for $U_i$, respectively. $D_{ij}$ is the Euclidean distance for $U_i$ to fly to task $T_j$, and $D_{i0}$ is the distance from its last task location back to the depot.
Management Cost: Cost related to UAV drone oversight and policy compliance.
$$
C_2 = \sum_{i=1}^{N} \sum_{j=1}^{M+1} \eta_i \cdot \frac{L_{ij}}{v_i} \cdot x_{ij} \qquad (4)
$$
where $\eta_i$ is the management cost coefficient for $U_i$, $L_{ij} = D_{ij} + D_{i0}$ is the total distance for a trip involving task $T_j$, and $v_i$ is the speed of $U_i$.
Fixed Operational Cost: Cost of depot infrastructure and basic operations.
$$
C_3 = \sum_{i=1}^{N} \varphi_i \cdot x_{ij} \qquad (5)
$$
where $\varphi_i$ is the fixed cost for deploying $U_i$.
The total economic cost is: $C = C_1 + C_2 + C_3$.
1.2.2 Delay Penalty ($S$)
A soft time window $[0, \kappa_j]$ is defined for each task $T_j$, where $\kappa_j$ is the latest acceptable delivery time. A penalty is incurred if service is completed after $\kappa_j$.
$$
S = \sum_{i=1}^{N} \sum_{j=1}^{M} S^i_{\tau_j} \qquad (6)
$$
where $S^i_{\tau_j}$ is the delay penalty for $U_i$ executing task $T_j$:
$$
S^i_{\tau_j} =
\begin{cases}
0 \cdot x_{ij}, & 0 < \tau_{ij} \le \kappa_j \\
\zeta \left[ \tau_{ij} – \kappa_j \right] \cdot x_{ij}, & \tau_{ij} > \kappa_j
\end{cases}
\qquad (7)
$$
Here, $\tau_{ij} = L_{ij}/v_i + \tau^p_j$ is the total time for $U_i$ to complete task $T_j$ and return, $\tau^p_j$ is the unloading time at $T_j$. $\zeta$ is a dynamic penalty coefficient that increases linearly with iteration $t$: $\zeta = \zeta_0 \cdot t / t_{max}$, encouraging timely solutions as the optimization progresses.
1.2.3 Safety Risk ($K$)
Safety risk is modeled considering UAV drone reliability and load balance. Higher load, longer distance, and imbalanced task allocation increase risk.
$$
K = (1 – K_r) \cdot (1 + \phi \cdot \varpi) \qquad (8)
$$
Reliability ($K_r$): Defined as the ratio of the sum of individual UAV drone reliabilities to the sum of their maximum ranges. The reliability of $U_i$ completing its assigned tour is modeled as a function of remaining range.
$$
K_r = \frac{ \sum_{i=1}^{N} K^i_r }{ \sum_{i=1}^{N} L_{i(max)} } = \frac{ \sum_{i=1}^{N} \prod_{j=1}^{M+1} \left(1 – \frac{L_{ij}}{L^j_{i(remain)}} \right) \cdot x_{ij} }{ \sum_{i=1}^{N} L_{i(max)} } \qquad (9)
$$
where $L_{i(max)}$ is the maximum range of $U_i$, and $L^j_{i(remain)} = L_{i(max)} – \sum_{m=1}^{j} L_{im} \cdot x_{im}$ is its remaining range after completing $j$ tasks.
Load Imbalance ($\varpi$): Measured by the coefficient of variation of the loads carried by all utilized UAV drones. Balanced loads reduce systemic risk.
$$
\varpi = \frac{\sigma_{load}}{\mu_{load}} = \frac{ \sqrt{ \frac{1}{N-1} \sum_{i=1}^{N} \left( q_i – \frac{1}{N}\sum_{i=1}^{N} q_i \right)^2 } }{ \frac{1}{N}\sum_{i=1}^{N} q_i } \qquad (10)
$$
where $q_i = \sum_{j=1}^{M} x_{ij} \cdot q_j$ is the total load of $U_i$, and $q_j$ is the demand weight of task $T_j$.
$\phi$ is a weighting coefficient for the load imbalance factor. Finally, each sub-objective is normalized to [0,1] using min-max normalization before being combined in Eq.(2).
1.3 Model Constraints
The model is subject to the following constraints:
Task Completion Constraint: All delivery tasks must be performed exactly once.
$$
\sum_{i=1}^{N} \sum_{j=1}^{M} x_{ij} = M \qquad (11)
$$
$$
\sum_{i=1}^{N} x_{ij} = 1, \quad \forall j = 1,2,…,M \qquad (12)
$$
Maximum Range Constraint: A UAV drone’s total flight distance cannot exceed its maximum range.
$$
\sum_{j=1}^{M+1} x_{ij} \cdot L_{ij} \le L_{i(max)}, \quad \forall i \qquad (13)
$$
Payload Capacity Constraint: A UAV drone’s total carried weight cannot exceed its maximum payload.
$$
q_i = \sum_{j=1}^{M} x_{ij} \cdot q_j \le q_{i(max)}, \quad \forall i \qquad (14)
$$
2. Multi-Strategy Improved Hippopotamus Optimization (MIHO) Algorithm
The Hippopotamus Optimization (HO) algorithm is a novel metaheuristic inspired by the foraging, defensive, and evasive behaviors of hippopotamuses. It demonstrates good local search capability and solution precision. However, to better solve the complex urban logistics task allocation problem, we propose three strategic improvements, resulting in the MIHO algorithm.
2.1 Latinized Stratified Sampling (LSS) for Population Initialization
Instead of random initialization or standard Latin Hypercube Sampling (LHS), we employ Latinized Stratified Sampling (LSS) to generate the initial population of hippopotamus individuals (solutions). LSS first divides each dimension’s normalized range [0,1] into $N$ equal strata. Rather than sampling randomly within each stratum, LSS deterministically selects the midpoint of each stratum. These midpoints are then randomly permuted across individuals for each dimension. This method guarantees perfect stratification and maximum minimum distance among samples, significantly improving population diversity and spatial coverage, providing a better starting point for the optimization of UAV drones task allocation.
2.2 Sine-Cosine Strategy for Position Update
To enhance the global exploration ability and help the algorithm escape local optima, we integrate strategies from the Sine Cosine Algorithm (SCA) into the position update formulas of HO.
In the early exploration phase (males competing for dominance), the original update $x^{M}_{ij} = x_{ij} + y_1 (D_{hippo} – I_1 x_{ij})$ is modified by incorporating sine and cosine functions:
$$
x^{M}_{ij} = x_{ij} + \frac{r_{s1}[\sin(r_{s2}) + \cos(r_{s2})]}{2} (D_{hippo} – I_1 x_{ij}) \qquad (15)
$$
In the later development phase (escaping predators), the update $x^{HE}_{ij} = x_{ij} + r_1 [ b^l_{j,local} + s_1 (b^u_{j,local} – b^l_{j,local}) ]$ is modified using a sine function for stable local exploitation:
$$
x^{HE}_{ij} = x_{ij} + r^*_{s1} \sin(r_{s2}) [ b^l_{j,local} + s_1 (b^u_{j,local} – b^l_{j,local}) ] \qquad (16)
$$
The parameters $r_{s1}$ and $r^*_{s1}$ are adaptive step-size coefficients, and $r_{s2}$ is a random angle, defined as:
$$
r_{s1} = \omega_{max} – (\omega_{max} – \omega_{min})\frac{t}{t_{max}}, \quad r_{s2} = 2\pi \cdot \text{rand}(), \quad r^*_{s1} = \omega_{max} \times \frac{t}{t_{max}}
$$
These sinusoidal perturbations introduce periodic oscillations, enabling a more thorough search of the solution space for optimal allocation of UAV drones.
2.3 Simulated Annealing (SA) Inspired Acceptance Criterion
To refine solution accuracy in the later stages and prevent premature convergence, we incorporate the Metropolis acceptance criterion from Simulated Annealing. After the standard HO position updates, a neighborhood perturbation is applied to the current solution $x_{curr}$ to generate a new candidate $x_{new}$. The new solution is accepted if it is better, or with a probability $p_{SA}$ if it is worse:
$$
p_{SA} = \begin{cases}
1, & F(x_{new}) \le F(x_{curr}) \\
\exp\left( -\frac{F(x_{new}) – F(x_{curr})}{T_{SA}} \right), & \text{otherwise}
\end{cases} \qquad (17)
$$
The temperature $T_{SA}$ cools down exponentially: $T_{SA} = T_0 \cdot \rho^{\,t-1}$, where $T_0$ is the initial temperature and $\rho$ is the cooling rate. This mechanism allows the algorithm to occasionally accept inferior solutions, facilitating escape from local optima and leading to more precise final solutions for the UAV drones scheduling problem.
3. Performance Validation of the MIHO Algorithm
3.1 Benchmark Function Tests
To verify the enhanced search capabilities of MIHO, we compare it against the original HO, Whale Optimization Algorithm (WOA), an improved Dung Beetle Optimizer (MIDBO), and two other HO variants (IHO1, IHO2) on a set of 9 benchmark functions, including unimodal, multimodal, hybrid, and composite types. Population size is 30, and maximum iterations are 500. The results, averaged over 30 independent runs, are summarized below.
| Function | MIHO | HO | WOA | MIDBO | IHO1 | IHO2 |
|---|---|---|---|---|---|---|
| F1 (Unimodal) | 1.47e+09 ± 3.66e+08 | 7.35e+10 ± 8.52e+09 | 1.08e+11 ± 1.10e+11 | 1.35e+11 ± 5.03e+10 | 1.52e+11 ± 1.13e+10 | 7.02e+10 ± 7.41e+09 |
| F2 (Multimodal) | 1.27e+03 ± 8.30e+01 | 1.15e+04 ± 1.68e+03 | 2.16e+04 ± 3.23e+03 | 2.55e+04 ± 1.36e+04 | 3.04e+04 ± 4.69e+03 | 1.13e+04 ± 2.16e+03 |
| F3 (Multimodal) | 6.01e+02 ± 4.77e-01 | 6.54e+02 ± 1.07e+01 | 6.72e+02 ± 1.34e+01 | 6.42e+02 ± 1.38e+01 | 6.59e+02 ± 1.37e+01 | 6.52e+02 ± 9.16e+00 |
| F4 (Hybrid) | 4.29e+08 ± 1.14e+08 | 1.65e+10 ± 3.67e+09 | 2.93e+10 ± 8.26e+09 | 2.45e+10 ± 7.70e+09 | 6.20e+10 ± 1.07e+10 | 1.45e+10 ± 4.50e+09 |
| F5 (Hybrid) | 2.04e+03 ± 1.10e+01 | 2.14e+03 ± 3.30e+01 | 2.24e+03 ± 6.42e+01 | 2.16e+03 ± 5.17e+01 | 2.17e+03 ± 3.02e+01 | 2.15e+03 ± 3.00e+01 |
| F6 (Hybrid) | 1.75e+06 ± 6.50e+05 | 4.45e+08 ± 3.15e+08 | 3.02e+09 ± 1.09e+09 | 2.01e+09 ± 1.54e+09 | 9.23e+09 ± 2.73e+09 | 4.28e+08 ± 3.10e+08 |
| F7 (Composite) | 3.39e+03 ± 7.90e+01 | 5.04e+03 ± 2.69e+02 | 5.31e+03 ± 2.73e+02 | 4.72e+03 ± 2.14e+02 | 4.81e+03 ± 1.59e+02 | 5.22e+03 ± 3.10e+02 |
| F8 (Composite) | 2.90e+03 ± 1.28e-04 | 3.10e+03 ± 8.87e+01 | 3.09e+03 ± 9.04e+01 | 3.01e+03 ± 4.78e+01 | 3.07e+03 ± 9.19e+01 | 3.11e+03 ± 9.83e+01 |
| F9 (Composite) | 5.53e+06 ± 6.35e+06 | 1.36e+09 ± 5.11e+08 | 2.75e+09 ± 1.21e+09 | 1.85e+09 ± 1.41e+09 | 7.51e+09 ± 3.27e+09 | 1.04e+09 ± 5.38e+08 |
The results clearly show that MIHO achieves the best (lowest) average fitness and standard deviation on all nine test functions, demonstrating its superior convergence accuracy, robustness, and ability to escape local optima compared to the other algorithms. The convergence curves further confirm that MIHO converges faster and to better values.
4. Simulation and Analysis for Multi-UAV Drones Task Allocation
4.1 Simulation Setup and Parameter Tuning
We consider a scenario with a depot at coordinates [125, 100, 25], $N=8$ available heterogeneous UAV drones, and $M=20$ delivery tasks. The parameters for the UAV drones and tasks are listed below. The algorithm parameters are: population size $n=30$, max iterations $t_{max}=200$, objective weights $\alpha_1=0.3, \alpha_2=0.3, \alpha_3=0.4$ (emphasizing safety). Other parameters like initial penalty $\zeta_0$, load balance weight $\phi$, and SA parameters $T_0, \rho$ were tuned via grid search.
| UAV $U_i$ | $\eta_i$ | $\lambda_{ij}$ | $\mu_{ij}$ | $\varphi_i$ | $v_i$ (m/s) | $q_{i(max)}$ (kg) | $L_{i(max)}$ (km) |
|---|---|---|---|---|---|---|---|
| U1 | 0.3 | 0.25 | 0.10 | 60.5 | 10 | 4.5 | 6.0 |
| U2 | 0.3 | 0.30 | 0.15 | 70.0 | 12 | 5.0 | 4.5 |
| U3 | 0.3 | 0.20 | 0.05 | 46.7 | 12 | 3.5 | 4.0 |
| U4 | 0.3 | 0.22 | 0.07 | 76.6 | 15 | 5.0 | 4.2 |
| U5 | 0.3 | 0.26 | 0.11 | 56.9 | 14 | 4.0 | 4.5 |
| U6 | 0.3 | 0.24 | 0.09 | 81.7 | 13 | 6.0 | 4.7 |
| U7 | 0.3 | 0.20 | 0.05 | 64.2 | 11 | 4.5 | 5.0 |
| U8 | 0.3 | 0.27 | 0.12 | 78.1 | 15 | 5.5 | 4.3 |
| Task $T_j$ | Coordinates | Time Window $[0, \kappa_j]$ (s) | Unload Time $\tau^p_j$ (s) | Demand $q_j$ (kg) |
|---|---|---|---|---|
| T1 | [75, 450, 50] | [0, 100] | 30 | 0.5 |
| T2 | [225, 400, 50] | [0, 180] | 60 | 1.0 |
| T3 | [1000, 750, 60] | [0, 200] | 60 | 1.0 |
| … | … | … | … | … |
| T20 | [1075, 425, 30] | [0, 570] | 66 | 1.1 |
4.2 Task Allocation Results and Analysis
The MIHO algorithm successfully solves the model, producing the following optimal task allocation scheme. UAV U8 is not utilized, and the tasks are balanced among the other seven UAV drones.
| UAV Drone | Task Execution Sequence | Flight Time (s) | Total Distance (m) | Total Load (kg) |
|---|---|---|---|---|
| U1 | T2 → T8 → T16 | 517.36 | 1746.03 | 4.5 |
| U2 | T9 → T10 → T4 | 502.62 | 2852.01 | 3.7 |
| U3 | T6 → T18 → T20 | 460.06 | 2068.59 | 3.4 |
| U4 | T11 → T1 → T17 | 351.51 | 1870.88 | 3.2 |
| U5 | T5 → T19 → T12 | 466.04 | 2409.90 | 3.8 |
| U6 | T13 → T15 | 266.16 | 802.52 | 2.9 |
| U7 | T3 → T14 → T7 | 549.56 | 2817.90 | 3.4 |
| U8 | — | — | — | — |
This allocation respects all constraints: each task is done once, and for each UAV drone, the total distance is within $L_{i(max)}$ and the total load is within $q_{i(max)}$. The load imbalance coefficient $\varpi$ for this plan is 0.432, indicating a reasonably balanced workload across the fleet of UAV drones, which contributes to lowering the overall safety risk.
4.3 Comparative Performance Evaluation
To further validate the effectiveness of MIHO for the specific problem of multi-UAV drones task allocation, we compare it against the same five algorithms (HO, WOA, MIDBO, IHO1, IHO2). Each algorithm was run 30 times independently on the urban logistics model. The performance metrics are the best, average, and standard deviation of the normalized total objective function $F^*$.
| Algorithm | Best $F^*$ | Average $F^*$ | Std. Dev. of $F^*$ |
|---|---|---|---|
| MIHO | 0.0380 | 0.0689 | 0.0159 |
| HO | 0.0870 | 0.1309 | 0.0217 |
| WOA | 0.1080 | 0.1787 | 0.0267 |
| MIDBO | 0.0823 | 0.1622 | 0.0299 |
| IHO1 | 0.1150 | 0.1644 | 0.0217 |
| IHO2 | 0.0859 | 0.1325 | 0.0210 |
MIHO outperforms all competitors across all three metrics. Compared to the original HO algorithm, MIHO reduces the best fitness by 56.32%, the average fitness by 47.36%, and the standard deviation by 26.73%. This demonstrates that the proposed improvements—LSS initialization, Sine-Cosine strategy, and SA acceptance—collectively endow the algorithm with a stronger ability to find higher-quality, more stable, and more reliable task allocation plans for a fleet of UAV drones operating in a complex urban logistics environment.
5. Conclusion
This paper addresses the multi-UAV task allocation problem in urban logistics delivery. We formulated a balanced task allocation model that simultaneously minimizes economic cost, delay penalty, and a safety risk measure incorporating load balance among UAV drones. To solve this complex model effectively, we proposed a Multi-strategy Improved Hippopotamus Optimization (MIHO) algorithm. Key improvements include Latinized Stratified Sampling for diverse population initialization, a Sine-Cosine strategy to enhance global exploration and escape local optima, and a Simulated Annealing acceptance criterion to refine solution accuracy. Extensive tests on benchmark functions confirm the superior search capability of MIHO. Simulation results for the urban logistics scenario demonstrate that the MIHO algorithm successfully generates feasible, efficient, and balanced task allocation plans for the UAV drones, significantly outperforming the original HO and other metaheuristic algorithms in terms of solution quality and robustness. The proposed model and algorithm provide an effective decision-support tool for optimizing UAV-based last-mile logistics operations.