Research on Vehicle-UAV Integrated Routing Optimization Problem to Deliver Medical Supplies

1. Introduction

The integration of technology within the healthcare sector has experienced pronounced advancement, mainly through the adoption of Unmanned Aerial Vehicles (UAVs) for the distribution of medical supplies to geographically isolated or inaccessible regions. Within the domain of emergency system planning, the imperative for a multifaceted and dependable transportation infrastructure capable of disseminating emergency supplies is accentuated. This infrastructure endeavors to harness the distinct benefits offered by various modes of transport in terms of experience and extensive coverage [1]. To ameliorate the proficiency of emergency materials transport, the employment of intelligent, high-technology distribution apparatuses, including UAVs, emerges as indispensable. The evolution of 5G, big data analytics, and the Internet of Things has facilitated the transition of UAVs from exclusively military utilization to civilian applications, endowing them with unique distributional benefits such as no-contact delivery, extended range capabilities, and utility in disaster-stricken contexts [2].

Empirical research has yielded optimized distribution schemes for UAV-based logistics, underscoring the pivotal role these systems play within emergency logistics by enabling rapid and precise supply distribution during crises. The endeavor to integrate vehicle-UAV routing systems principally targets reducing delivery times from healthcare facilities to strategically positioned vehicles within road networks, thereby ensuring efficacious coordination with terrestrial vehicle convoys. This methodology seeks to mitigate a substantial fraction of healthcare facilitie’s operational expenditures, with the objective of reallocating saved resources towards the amplification of patient care services [3,4,5].

The practical deployment of UAVs within the logistics sector has witnessed considerable advancement. In 2013, Amazon’s proclamation of its UAV delivery service marked a pivotal moment, with a successful geographical positioning system (GPS)-guided supply deliveries by 2016 [6]. In 2015, France Post achieve a fully autonomous package delivery over a 9-mile distance, thereby showcasing the transformative potential of UAV technology within logistics [7]. Subsequent to receiving certification for the operation of UAVs beyond visual line of sight (BVLOS), UPS unveiled new prospects for autonomous deliveries across extensive areas [8,9]. Alphabet’s Wing has been at the forefront of UAV delivery system innovation, aiming at augmenting efficiency, expedition, and environmental sustainability, with operations spanning the United States, Australia, and Finland. Wing’s operational model caters to a wide array of consumer requirements, encompassing food and pharmacy products [10]. Zipline has specialized in the healthcare segment, leveraging UAVs to dispatch essential medical supplies to remote locales. Through its operations in Rwanda and Ghana, Zipline has revolutionized medical supply logistics, facilitating timely deliveries of blood, vaccines, and other critical supplies [11]. The Parcel-copter project by DHL further epitomizes the convergence of UAV technology with traditional logistics frameworks aimed at surmounting geographical and infrastructural barriers in delivery services [12].

These developments underscore the transformative potential of UAVs in refining healthcare logistics and establishing a resilient and efficient transport system for both emergency and routine medical supplies. This article delves into the intricacies of integrated vehicle-UAV delivery systems, with the objective of enhancing healthcare access through the timely and efficient distribution of crucial medical resources. The emphasis on developing a routing system that simultaneously minimizes delivery times and costs while elevating the efficiency and resilience of healthcare logistics contributes significantly to the academic discourse by advancing the understanding of spatiotemporal routing within UAV-to-ground vehicle networks.

2. Literature Review

2.1. Current Research Status of Emergency Material Distribution

After an emergency occurs, emergency rescue supplies must be safely and quickly delivered. Various provinces and cities have issued multiple guidance opinions on emergency material support, and domestic and foreign scholars have also been paying close attention to research in the field of emergency materials. In 2022, Jung and Kim [14] designed two different types of decision variables, operation sequence, and current battery level, to transport emergency packages to remote islands using UAVs combined with wind direction and speed conditions. Hu Zhongjun [15] studied the problem of emergency material dispatch and transportation allocation after urban flood disasters in 2018. Deng Xiuqin [16] studied the establishment of an emergency logistics distribution and material support system under COVID-19 in 2020. Julie [17] studied the dynamic allocation of emergency supplies in 2020 and proposed a plan for post-disaster emergency supplies allocation and vehicle distribution paths. Qin Jin [18] studied the emergency support of people’s livelihood materials in significant emergencies in 2023. Yu et al. [19] proposed an emergency material distribution route optimization problem with the goal of minimizing the total time, which cannot be predicted before the delivery of rescue demand sequences. They provided three online strategies and their competitive ratios. Wu et al. [20] proposed an online strategy and proved the lower bound of competition ratio for the optimization problem of emergency material distribution that incurs time costs if immediate service is not available. Su et al. [21] established an online selection model for emergency material distribution paths in response to the situation where each demand point on the network sends out emergency material distribution service requests with delivery time requirements. However, the distribution center cannot know the time and location of the request in advance.

2.2. Current Research Status of Vehicle-UAV Integrated Delivery Path Optimization

In 2014, AMP Company disclosed its research on the joint distribution system of "electronic delivery vehicles+UAVs.". Wohlsen [22] elaborated on AMP's vision of using trucks for delivery while UAVs can deliver independently, and trucks can support UAV charging and safe landing. Yurek and Ozmutlu [23] investigated several logistics companies deploying UAVs with the aim of improving efficiency and shortening delivery time. She studied the integrated service of customers by installing a UAV on the roof of a vehicle, which is considered a supplement to the vehicle and is defined as the Traveling Salesman Problem UAV problem for modeling. Freitas et al. [24] extended the TSP problem by considering constraints such as the flight time limit of UAVs and the package not exceeding the payload of UAVs. Sacramento [25] defines a Vehicle Routing Problem (VRP) as a situation where two trucks each carry a UAV, and the UAV serves as direct access to a single customer. As a variant of TSP-D, a model is constructed to reduce fuel consumption costs by reducing truck operating time. Cavani [26] mainly considers the synchronization between trucks and multiple UAVs in the problem of finding the shortest duration path for trucks and UAVs to serve all customers. He constructed a compact mixed integer linear programming model. Murray's [27] research deploys multiple UAVs from trucks to serve long-distance customers, assuming that trucks can launch UAVs from different locations but can only launch UAVs carrying a single package at a time and recycle them at different service points such as warehouses, in order to shorten the total delivery time between UAVs and trucks under dynamic collaboration as the goal and establish a mathematical model. Under the "vehicle+multi UAV" joint distribution model, Liu Wusheng [28] uses UAVs as the decision-making leader to allocate routes in three stages. The UAV can be extended to deliver multiple demand points in a single service, and a model with the goal of minimizing delivery distance is established. Leibo et al. [29] adopted the new mode of "truck UAV" combined distribution under the background of COVID-19 to implement material distribution in the epidemic area and further explored the distribution value of logistics when this mode affects COVID-19. Tamke [30] developed a comprehensive mixed integer problem that considers speed-dependent energy consumption and suggests performing different flights at different speeds instead of continuously operating the UAV at maximum speed. On the contrary, choosing a model that balances the range and delivery speed of UAVs with the speed of UAVs aims to minimize the operating costs composed of truck fuel consumption costs, driver labor costs, and UAV energy costs. Mohammad et al. [31] proposed a mixed integer linear programming (MILP) model, which determines the optimal allocation of trucks and UAVs by customers, the optimal route order of trucks, and the optimal launch and reassembly positions of UAVs on the truck route. The model combines truck and UAV (UAV) operations to optimize the planning of distribution routes in a multimodal system. In the system, truck and UAV operations are synchronized, meaning that one or more UAVs are carried by a truck, and the vehicle acts as a mobile warehouse. UAVs and trucks can both deliver goods. Although the truck follows a multi-stop route, each UAV delivers goods once at a time [32,33]. The proposed optimization model minimizes the waiting time for customers in the system.

2.3. Current Research of Solving Methods for Routing Optimization Problem

There are various methods for solving optimization problems, including accurate algorithms, classic heuristic algorithms, intelligent algorithms, etc. Accurate algorithms are crucial in solving complex optimization problems by providing optimal solutions. These algorithms require rigorous mathematical models to ensure precise solutions within manageable computational complexity. Various algorithms exhibit distinct features and applicability. For instance, the branch and bound method, introduced by Stephen and Jacob [34], utilizes a tree-like structure to seek the optimal solution iteratively and is helpful in solving delivery problems with up to 65 customers. Dynamic programming algorithms, such as those employed by Li Yanfeng et al. [35] and Ouwei et al. [36], trade space for time, managing temporal complexity effectively, and providing superior solutions compared to traditional algorithms like nearest neighbor. These algorithms are essential in solving large-scale optimization challenges efficiently and effectively, providing optimal solutions even when faced with considerable complexity and uncertainty. Heuristic algorithms [37] are practical tools that can tackle large-scale problems in real-world scenarios. Although they are less precise than classical algorithms, they offer many advantages. There are two types of classic heuristic algorithms: constructive algorithms and two-stage algorithms. Constructive algorithms, such as saving algorithms and Solomon’s [38] insertion method, build solutions incrementally. Two-stage algorithms generate feasible solutions, such as Lin’s [39] 3-opt algorithm and Vigo’s [40] clustering design path algorithm, and then optimize them. One of the most prominent heuristic algorithms is the taboo search algorithm. It mimics human intelligence by iteratively refining solutions and restricting the search space to promising regions. Scholars have summarized the differences and diversity preserved in the constantly evolving VRP models.

3. Research Methodology

The distribution model of UAV and vehicles differs greatly from the problem of single agent distribution and existing research on dual agent distribution in terms of optimization objectives and constraints. Therefore, the models constructed in the past are no longer applicable to current problems. On the basis of previous research, we further clarify important factors such as constraints and objective functions, considers the constraints that UAVs can deliver in parallel in the distribution area and a single UAV can deliver multiple customers with material needs, and establishes a model with the shortest total completion time and the highest customer satisfaction.

3.1. Modeling Ideas

The modeling problem central to integrating Unmanned Aerial Vehicles (UAVs) with ground vehicles for medical supply delivery aims to address intrinsic deficiencies in healthcare delivery mechanisms.

Given a transport vehicle equipped with a set of isomorphic UAVs $\mathrm{U}=\{{\mathrm{U}}_1,{\mathrm{U}}_2,...{\mathrm{U}}_{\mathrm{k}}\}$, K is the number of UAVs. The task that needs to be executed is to distribute emergency supplies from the starting point warehouse (set as ${\mathrm{b}}_0$) to the set of customers with given material needs, represented as C = {1, 2... n}, and return to the destination warehouse (set as ${\mathrm{b}}_{\mathrm{n}+1}$) after completing the distribution task. The demand for each material demand customer is expressed in g, with vehicles traveling at an average speed of ${\mathrm{V}}^{\mathrm{T}}$ and UAVs flying at an average speed of ${\mathrm{V}}^{\mathrm{D}}$. In our problem setting, the customers of material demand points are divided into two main parts: the UAV delivery customer set and the vehicle delivery customer set. Establish a path optimization model for this type of dual subject delivery problem based on meeting corresponding constraint conditions, optimize the total completion time and customer satisfaction objectives, and seek the optimal integrated delivery path.

Figure 1. Delivery sketch map.

Figure 1. Delivery sketch map.

3.1.1. Objective Functions

Our article has two main objectives: one is to minimize delivery time and maximize customer satisfaction.

Objective 1: Minimizing the Delivery Time

The driving plan of a vehicle is actually an orderly arrangement of customers in the vehicle delivery customer set ${\mathrm{C}}^{\mathrm{T}}$, that is, the vehicle delivery route is marked as $\{{\mathrm{b}}_0,{\mathrm{b}}_1...{\mathrm{b}}_{{\mathrm{n}}^{\mathrm{r}}}\}$. To ensure that vehicles complete emergency supplies distribution and return to the warehouse in the shortest possible time, the optimal sequence for vehicle delivery to customers is determined. We introduce a variable ${\mathrm{x}}_{\mathrm{i}\mathrm{j}}$,where i, j ∈ ${\mathrm{C}}^{\mathrm{T}}$, indicating whether the vehicle is traveling from material demand point $\mathrm{i}$ to point$\mathrm{j}$. If so, then ${\mathrm{x}}_{\mathrm{i}\mathrm{j}}=1$, otherwise it is 0. The time taken for the vehicle to travel from point $\mathrm{i}$to point $\mathrm{j}$ is:

$${\mathrm{t}}_{\mathcal{i}\mathcal{j}}^{\mathrm{T}}={\displaystyle \frac{{\mathrm{L}}_{\mathcal{i}\mathcal{j}}^{\mathrm{T}}}{{\mathrm{V}}^{\mathrm{T}}}}$$

(1)

The time for vehicle delivery to customers can be expressed as:

$${\mathrm{t}}^{\mathrm{T}}=\sum _{(\mathcal{i},\mathcal{j}\in {\mathrm{C}}^{\mathrm{T}})}{\mathrm{t}}_{\mathcal{i}\mathcal{j}}^{\mathrm{T}}{\mathrm{x}}_{\mathcal{i}\mathcal{j}}$$

(2)

In this issue, if the customer’s demand for UAV delivery involves parallel delivery, the completion time of delivery is the longest time for the UAV to return to the truck among all UAV delivery paths. Due to the use of homogeneous UAVs, the payload of the UAV is set to${\mathrm{Q}}_{\mathrm{k}}^{\mathrm{D}}$ $(\mathrm{k}=\mathrm{1,2}...\mathrm{K})$. Considering the battery capacity of the UAV, the maximum distance for a single flight of the UAV is L, and the Euclidean distance between customer points for material needs is ${\mathrm{L}}_{\mathcal{i}\mathcal{\text{'}}\mathcal{j}}^{\mathrm{D}}$, and

$${\mathrm{L}}_{\mathcal{i}\mathcal{\text{'}}\mathcal{j}}^{\mathrm{D}}=\sqrt{{({\mathrm{A}}_{\mathcal{i}\mathcal{\text{'}}}-{\mathrm{A}}_{\mathcal{j}\mathcal{\text{'}}})}^2+{({\mathrm{B}}_{\mathcal{i}\mathcal{\text{'}}}-{\mathrm{B}}_{\mathcal{j}\mathcal{\text{'}}})}^2}(\mathcal{i},\mathcal{j}\in {\mathrm{C}}^{\mathrm{D}})$$

(3)

The distance from the UAV takeoff point to the material demand point is ${\mathrm{L}}_{{\mathrm{i}}^{\mathrm{*}}\mathrm{j}\mathrm{\text{'}}}({\mathrm{i}}^{\mathrm{*}}\in {\mathrm{C}}^{\mathrm{T}},{\mathrm{j}}^{\mathrm{\text{'}}}\in {\mathrm{C}}^{\mathrm{D}})$, and the nth flight route (i.e. the number of UAV flights) is defined using the path set ${\mathrm{P}}_{\mathrm{n}}$. The element ${\mathrm{P}}_{{\mathrm{n}}_{\mathrm{i}}}$ in the route represents the order of customer i being delivered by the UAV in route n. If the m-th customer of the UAV ${\mathrm{U}}_{\mathrm{k}}$ is set to ${\mathrm{k}}_{\mathrm{m}}$ (${\mathrm{k}}_{\mathrm{m}}$=0, which means that the m-th customer has not been delivered, and ${\mathrm{P}}_{{\mathrm{n}}_{{\mathrm{i}}^{\mathrm{*}}}}$ represents the launch and recovery points of the UAV, then the completion time of one operation flight for the UAV delivery group ${\mathrm{G}}_{{\mathrm{N}}_1}$ is:

$${\mathrm{t}}_{{\mathrm{G}}_{\mathrm{i}}}^{\mathrm{D}}=\mathrm{m}\mathrm{a}\mathrm{x}\left({\displaystyle \frac{{\mathrm{x}}_{\mathrm{i}\mathrm{\text{'}}\mathrm{j}\mathrm{\text{'}}}\left(\sum _{\mathrm{i}\mathrm{\text{'}}=1,\mathrm{j}\mathrm{\text{'}}={\mathrm{k}}_{\mathrm{m}}}{\mathrm{L}}_{{\mathrm{P}}_{{\mathrm{n}}_{\mathrm{i}}},{\mathrm{n}}_{\mathrm{j}}}^{\mathrm{D}}+{\mathrm{L}}_{{\mathrm{P}}_{{\mathrm{n}}_{{\mathrm{i}}^{\mathrm{*}}}},{\mathrm{n}}_{\mathrm{i}}}^{\mathrm{D}}+{\mathrm{L}}_{{\mathrm{P}}_{{\mathrm{n}}_{\mathrm{j}}},{\mathrm{n}}_{{\mathrm{i}}^{\mathrm{*}}}}^{\mathrm{D}}\right)}{{\mathrm{V}}^{\mathrm{D}}}}\right),\mathrm{n}\in \mathrm{1,2}....\mathrm{n}$$

(4)

The completion time for UAV delivery to discrete customers is:

$${\mathrm{t}}_{{\mathrm{j}}^{\mathrm{*}}}=2{\mathrm{f}}_{{\mathrm{j}}^{\mathrm{*}}}.{\displaystyle \frac{{\mathrm{L}}_{{{\mathrm{i}}^{\mathrm{*}}\mathrm{j}}^{\mathrm{*}}}}{{\mathrm{V}}^{\mathrm{D}}}}$$

(5)

So the total time to complete all grouping and discrete delivery customers is:

$${\mathrm{t}}^{\mathrm{D}}=\sum _{\mathrm{i}=1}^{\mathrm{N}}{\mathrm{t}}_{{\mathrm{G}}_{\mathrm{i}}}^{\mathrm{D}}+\sum _{{\mathrm{j}}^{\mathrm{*}}\in {\mathrm{C}}^{\mathrm{D}}\cap {\mathrm{j}}^{\mathrm{*}}\notin \mathrm{G}}{\mathrm{t}}_{{\mathrm{j}}^{\mathrm{*}}}$$

(6)

The main objective of healthcare material distribution is to deliver emergency materials to customers in demand faster. The equivalent form of the optimization objective function is as follows:

$${\mathrm{F}}_1=\mathrm{m}\mathrm{i}\mathrm{n}({\mathrm{t}}^{\mathrm{T}}+{\mathrm{t}}^{\mathrm{D}})$$

(7)

Objective 2: Maximizing Customer Satisfaction

When a public event occurs, the requirements of the delivered customers, especially their satisfaction with the delivery time, also need to be taken seriously.

The time satisfaction is represented by a number between 0 and 1. 0 represents the time when the vehicle or UAV reaches the customer’s demand point as very dissatisfied, while 1 represents the time when the vehicle or UAV reaches the customer’s demand point as very satisfied. Assuming that the latest expected time point for vehicles or UAVs to arrive at the material demand customer ${\mathrm{i}}_0$ is ${\mathrm{t}}_{{\mathrm{i}}_0}^h$, and the maximum acceptance time range for the material demand customer is $[0,{\mathrm{T}}_{{\mathrm{i}}_0}^h]$, the time satisfaction is shown in Figure 2.

The time satisfaction equation can be expressed as follows:

$${\mathrm{S}}_1=\left\{\begin{array}{c}1\\ {\displaystyle \frac{{\mathrm{t}}_{{\mathrm{i}}_0}^h-{\mathrm{t}}_{{\mathrm{i}}_0}}{{\mathrm{t}}_{{\mathrm{i}}_0}^h-{\mathrm{T}}_{{\mathrm{i}}_0}^h}}\\ 0\end{array},{{\mathrm{t}}_{{\mathrm{i}}_0}^h\le {\mathrm{t}}_{{\mathrm{i}}_0}\le \mathrm{T}}_{{\mathrm{i}}_0}^h\right.$$

(8)

From the perspective of distribution material quantity, satisfaction is positively correlated with the satisfaction rate of distribution material quantity. Using 0 to indicate that the number of materials delivered by UAVs or vehicles to customers in need completely does not meet the demand, and the demand customers are very dissatisfied. Using 1 to indicate that the number of materials delivered to customers in need completely meets the demand, and the demand customers are very satisfied. The satisfaction chart regarding the satisfaction rate of material delivery volume is as follows Figure 3:

The relationship expression is as follows:

$${\mathrm{S}}_2=\left\{\begin{array}{c}0\\ {\displaystyle \frac{{\mathrm{u}}_{{\mathrm{i}}_0}-{\mathrm{u}}_h}{{\mathrm{u}}_{\mathrm{l}}-{\mathrm{u}}_{{\mathrm{i}}_0}}}\\ 1\end{array},{\mathrm{u}}_{\mathrm{l}}\le {\mathrm{u}}_{{\mathrm{i}}_0}\le {\mathrm{u}}_h\right.$$

(9)

The optimal objective function for customer satisfaction obtained by combining the time satisfaction sector and the quantity satisfaction is:

$${\mathrm{F}}_2=\mathrm{m}\mathrm{a}\mathrm{x}(\mathsf{\alpha }{\mathrm{S}}_1+\mathsf{\beta }{\mathrm{S}}_2)$$

(10)

$${\mathrm{s}.\mathrm{t} \mathrm{g}}_{\mathrm{i}}^{\mathrm{D}}\mathrm{o}\mathrm{r}{\mathrm{g}}_{\mathrm{i}}^{\mathrm{T}}\le {\mathrm{g}}_{\mathrm{i}}$$

(11)

$${{\mathrm{t}}_{{\mathrm{i}}_0}\le \mathrm{T}}_{{\mathrm{i}}_0}^h$$

(12)

The demand and time are equally important, here α Related to β Taking the mean of 0.5, the constraint (27) indicates that the number of materials delivered by vehicles or UAVs to the material demand point cannot exceed the customer’s demand, and the constraint (28) indicates that the time for delivery to the material demand point cannot exceed the maximum tolerance time.

3.1.2. Constraints

The outlined constraints are in place to ensure the material supply process is efficient and effective in a integrated delivery system between vehicles and UAVs. These constraints dictate the behavior of the vehicles and UAVs, optimizing the delivery process and enhancing overall system efficiency.

$$\sum _{\mathrm{i}\in {\mathrm{C}}^{\mathrm{T}}}{\mathrm{x}}_{0\mathrm{i}}=\sum _{\mathrm{i}\in {\mathrm{C}}^{\mathrm{T}}}{\mathrm{x}}_{\mathrm{i}0}=1$$

(13)

$$\sum _{\mathrm{i}\in {\mathrm{C}}^{\mathrm{T}}\cup \left\{0\right\},\mathrm{j}\in {\mathrm{C}}^{\mathrm{T}}\cup \{\mathrm{n}+1\}}{\mathrm{x}}_{\mathrm{i}\mathrm{j}}=\sum _{\mathrm{i}\in {\mathrm{C}}^{\mathrm{T}}\cup \left\{0\right\},\mathrm{j}\in {\mathrm{C}}^{\mathrm{T}}\cup \{\mathrm{n}+1\}}{\mathrm{x}}_{\mathrm{i}\mathrm{j}}$$

(14)

$$\sum _{\mathrm{i}\in {\mathrm{C}}^{\mathrm{T}}\cup \left\{0\right\},\mathrm{j}\mathrm{\text{'}}\in {\mathrm{C}}^{\mathrm{D}}}{\mathrm{x}}_{\mathrm{i}\mathrm{j}\mathrm{\text{'}}}=0$$

(15)

$$\sum _{\mathrm{i}\in {\mathrm{C}}^{\mathrm{T}}}{\mathrm{g}}_{\mathrm{i}}^{\mathrm{T}}+\sum _{\mathrm{i}\in {\mathrm{C}}^{\mathrm{D}}}{\mathrm{g}}_{\mathrm{i}}^{\mathrm{D}}\le {\mathrm{Q}}^{\mathrm{T}}$$

(16)

$$\sum _{\mathrm{i},\mathrm{j}\in {\mathrm{C}}^{\mathrm{T}}}{\mathrm{x}}_{\mathrm{i}\mathrm{j}}+\sum _{\mathrm{i}\mathrm{\text{'}},\mathrm{j}\mathrm{\text{'}}\in {\mathrm{C}}^{\mathrm{D}}}{\mathrm{x}}_{\mathrm{i}\mathrm{\text{'}}\mathrm{j}\mathrm{\text{'}}}=1$$

(17)

$$\sum _{\mathrm{i}\mathrm{\text{'}},\mathrm{j}\mathrm{\text{'}}\in {\mathrm{C}}^{\mathrm{D}}}{\mathrm{x}}_{\mathrm{i}\mathrm{\text{'}}\mathrm{j}\mathrm{\text{'}}}=\sum _{\mathrm{i}\mathrm{\text{'}},\mathrm{j}\mathrm{\text{'}}\in {\mathrm{C}}^{\mathrm{D}}}{\mathrm{x}}_{\mathrm{i}\mathrm{\text{'}}\mathrm{j}\mathrm{\text{'}}}$$

(18)

$$\sum _{\mathrm{i}\mathrm{\text{'}}=1,\mathrm{j}\mathrm{\text{'}}={\mathrm{n}}_{\mathrm{m}},\mathrm{n}\in {\mathrm{P}}_{\mathrm{n}}^{\mathrm{D}}}{\mathrm{L}}_{{\mathrm{P}}_{{\mathrm{n}}_{\mathrm{i}}},{\mathrm{n}}_{\mathrm{j}}}^{\mathrm{D}}+{\mathrm{L}}_{{\mathrm{P}}_{{\mathrm{n}}_{{\mathrm{i}}^{\mathrm{*}}}},{\mathrm{n}}_{\mathrm{i}}}^{\mathrm{D}}+{\mathrm{L}}_{{\mathrm{P}}_{{\mathrm{n}}_{\mathrm{j}}},{\mathrm{n}}_{{\mathrm{i}}^{\mathrm{*}}}}^{\mathrm{D}}<\mathrm{L}$$

(19)

$$h_{{\mathrm{i}}^{\mathrm{*}}\mathrm{j}}^{\mathrm{p}}=1,{\mathrm{i}}^{\mathrm{*}}\in {\mathrm{C}}^{\mathrm{T}}\cap {\mathrm{i}}^{\mathrm{*}}\ne 0,\mathrm{j}\mathrm{\text{'}}\in {\mathrm{C}}^{\mathrm{D}}$$

(20)

$$h_{{\mathrm{i}}^{\mathrm{*}}\mathrm{j}}^{\mathrm{p}}+{\mathrm{x}}_{\mathrm{j}\mathrm{\text{'}}\mathrm{k}\mathrm{\text{'}}}\ge 1,{\mathrm{i}}^{\mathrm{*}}\in {\mathrm{C}}^{\mathrm{T}},\mathrm{j}\mathrm{\text{'}}\in {\mathrm{C}}^{\mathrm{D}},\mathrm{k}\mathrm{\text{'}}\in {\mathrm{C}}^{\mathrm{D}}$$

(21)

$$\sum _{\mathrm{j}\mathrm{\text{'}}\in {\mathrm{C}}^{\mathrm{D}}}{h}_{{\mathrm{i}}^{\mathrm{*}}\mathrm{j}}^{\mathrm{p}}=1,\forall {\mathrm{i}}^{\mathrm{*}}\in {\mathrm{C}}^{\mathrm{T}}$$

(22)

$$\sum _{\mathrm{i}\in {\mathrm{C}}^{\mathrm{T}}}{\mathrm{x}}_{\mathrm{i}\mathrm{j}}=h_{{\mathrm{i}}^{\mathrm{*}}\mathrm{j}}^{\mathrm{p}},{\mathrm{j}=\mathrm{j}}^{\mathrm{*}}\in {\mathrm{C}}^{\mathrm{T}}$$

(23)

$$\sum _{\mathrm{i}\mathrm{\text{'}}\in {\mathrm{C}}^{\mathrm{D}}}{\mathrm{g}}_{\mathrm{i}\mathrm{\text{'}}}^{\mathrm{D}}\le {\mathrm{Q}}^{\mathrm{D}}$$

(24)

$${\mathsf{\alpha }}_{\mathrm{i}}-{\mathsf{\alpha }}_{\mathrm{j}}+1\le (\mathrm{n}+2)+2(1-{\mathrm{x}}_{\mathrm{i}\mathrm{j}})\forall \mathcal{i}\in \{0....\mathrm{n}\},\forall \mathrm{j}\in \{1....\mathrm{n}+1\},1\le {\mathsf{\alpha }}_{\mathrm{i}}\le \mathrm{n}+2$$

(25)

$${\mathsf{\alpha }}_{\mathrm{i}}-{\mathsf{\alpha }}_{\mathrm{j}}+1\le \mathrm{n}(1-{\mathrm{x}}_{\mathrm{i}\mathrm{\text{'}}\mathrm{j}\mathrm{\text{'}}})\forall \mathrm{i}\mathrm{\text{'}},\mathrm{j}\mathrm{\text{'}}\in {\mathrm{C}}^{\mathrm{D}}1\le {\mathsf{\alpha }}_{\mathrm{i}}\le \mathrm{n}$$

(26)

$${\mathrm{w}}_{\mathrm{i}\mathrm{\text{'}}{\mathrm{i}}^{\mathrm{*}}\mathrm{j}\mathrm{\text{'}}}^{\mathrm{u}}=0$$

(27)

$${\mathrm{t}}_{\mathrm{j}}^{\mathrm{T}}={\mathrm{t}}_{\mathrm{i}}^{\mathrm{T}}+{\mathrm{t}}_{{\mathrm{i}}^{\mathrm{*}}\mathrm{j}}^{\mathrm{T}}+{\mathrm{f}}_{{\mathrm{j}}^{\mathrm{*}}}{\mathrm{t}}_{{\mathrm{j}}^{\mathrm{*}}}^{\mathrm{D}}+h_{{\mathrm{i}}^{\mathrm{*}}\mathrm{j}}^{\mathrm{p}},{\mathrm{t}}_{{\mathrm{G}}_{\mathrm{i}}}^{\mathrm{D}},{\mathrm{i}=\mathrm{i}}^{\mathrm{*}}\in {\mathrm{C}}^{\mathrm{T}}$$

(28)

$${\mathrm{t}}_{\mathrm{j}\mathrm{\text{'}}}^{\mathrm{D}}={\mathrm{t}}_{\mathrm{i}\mathrm{\text{'}}}^{\mathrm{D}}+{\mathrm{t}}_{\mathrm{i}\mathrm{\text{'}}\mathrm{j}\mathrm{\text{'}}}^{\mathrm{D}}$$

(29)

4. Results and Discussions

Optimization problems have always been characterized by high complexity and difficulty in solving. With the in-depth research of algorithms, intelligent heuristic algorithms have gradually become a new method for solving optimization problems. This section first preprocesses the clustering of UAV delivery sets, and then designs artificial bee colony algorithms to optimize vehicle delivery paths and UAV delivery paths respectively.

4.1. Data Preprocessing for UAV Delivery to Customers

Due to the characteristics of large quantities and clustering when public events occur, this article adopts the DBSCAN algorithm to cluster the unmanned aerial vehicle (UAV) delivery material demand customers that cannot be reached by vehicles according to the density. DBSCAN (Density Based Spatial Clustering of Applications with Noise) is a density based data clustering method proposed by Martin Ester et al. [54] in 1996. This algorithm has certain advantages over many clustering algorithms:

It can divide high-density areas into clusters, and can also divide clusters of different shapes in noisy data.

Compared with the K-MEANS clustering method, there is no need to input the number of clusters to be divided in advance.

You can input parameters to filter noise when needed.

After clustering, the material demand customers are divided into different groups and individual discrete customers, and each discrete customer is treated as a separate category. At the same time, the closest vehicle delivery customer point to the clustered customer is identified as the density corresponding customer, which is used as the location for the launch and recovery of the vehicle mounted UAV.

4.1.1. UAV Customer Clustering Processing Based on DBSCAN Algorithm

The specific steps of DBSCAN algorithm for clustering UAV delivery customers are as follows:

Read coordinate data of material demand customer points from the dataset

Input the required parameters in the DBSCAN algorithm

Output the results of clustering and the number of clusters

Determine the coordinates of the closest distance between vehicle delivery customers and customers in each category after clustering

Determine the position coordinate point of the vehicle as the location for UAV launch and retrieval

The DBSCAN algorithm flow for customer clustering processing in UAV delivery is shown in Figure 4.

4.1.2. Dual Objective Processing in the Model

The model constructed in this article includes two objectives: (1) the shortest total delivery time and (2) the highest customer satisfaction. Based on the high complexity of the dual objective function in the optimization process, the dual objective function model is transformed into a single objective function optimization model for solution. Given that the approach and operation of adding penalty factors are both easy, this method is adopted in the article to handle the dual objectives in the model. The specific processing steps are as follows: Step 1: Since the highest customer satisfaction is a maximum problem, directly adding a penalty factor will not match the total delivery time of the minimum problem and does not conform to logical reasoning. Therefore, taking a negative customer satisfaction value and adding a time penalty factor will result. Then add it linearly to the total delivery time.

Step 2: Taking a negative value for customer satisfaction $\left[\mathrm{0,1}\right]$ may lead to an imbalance between the two objective functions, resulting in a negative total objective function value. Therefore, it is necessary to select an appropriate time penalty factor $\mathsf{\gamma }\in \left[\mathrm{0,1}\right]$. In the text, the average satisfaction value is set at 0.5. When the satisfaction is greater than $0.5,\mathsf{\gamma }$ take 0.25; When satisfaction is less than $0.5,\mathsf{\gamma }$ take 0.75. Make both objectives achieve good results in optimization. However, the optimization direction of the dual objective function is different, and the optimization of objective F1 will inevitably affect the optimization of objective F2. In the end, a satisfactory solution can only be obtained according to actual requirements.

The final optimization objective function after processing is:

$\mathrm{m}\mathrm{i}\mathrm{n}\mathrm{F}={\mathrm{F}}_1-\mathsf{\gamma }{\mathrm{F}}_2$ (Where F is the transformed objective function)

The rate of return is positively correlated with the probability of selecting a food source, so the probability of bees being recruited to which food source is also positively correlated with the rate of return on the food source 55,56]. The process of artificial bee colony algorithm is shown in Figure 5.

The ABC algorithm has the characteristics of fewer parameters, simplicity, and easy implementation. In the past decade, ABC has been successfully applied to solve TSP problems, and its powerful optimization capabilities and outstanding search advantages in TSP and VRP have been tested. TSP-MD is an extension of TSP, and the close relationship between TSP and TSP-MD lays the foundation for handling TSP-MD.

4.1.3. Artificial Bee Colony Algorithm Initialization and Neighborhood Search Strategy

The UAV Road section is the path generated by dividing the number of UAVs into individual units in a group. The optimal solution of the vehicle and the optimal solution of the UAV are combined using cross chain coding to form a joint delivery path, as shown in Figure 6:

When the artificial bee colony algorithm executes a process, initialization is the first crucial step, which includes the population size $\mathrm{N}\mathrm{P}$, the number of iterations required, and the maximum search limit. The initial solution is generated through a random method, which is the path sequence $\mathrm{X}\mathrm{i}(\mathrm{i}=\mathrm{1,2},3,...,\mathrm{S}\mathrm{N})$ from the starting warehouse to the visiting customer and then back to the warehouse. SN is the number of food sources that need to be delivered to the customer. In the search process, neighborhood search strategies are divided into three methods: exchange, insertion, and reverse order, as shown in Figure 7.

4.1.4 Decoding Strategy for UAV Path

The initial solution is constructed using a UAV insertion path strategy. If $\mathrm{N}$ UAVs simultaneously deliver $\mathrm{M}$ customers, the decoding diagram of element $\mathrm{N}+\mathrm{M}-1$ in the obtained solution is shown in Figure 8:

4.1.5. Adaptive Probability Design for Following Bees

This section describes the design of adaptive probability for following bees in artificial bee colony algorithms. The common feature of heuristic algorithms is that they are selected according to a certain probability, and how the selected probability is allocated has a significant impact on the optimization efficiency of the algorithm. Optimizing randomness is one of the characteristics of artificial bee colony algorithms, and the follower bees in the algorithm heavily rely on probability to choose when searching for the optimal process. In the ABC algorithm, the follower bee adopts a roulette wheel method to select the leader bee. Its characteristic is that as the fitness value increases, the probability of being selected increases, which is the basic probability selection, as shown in Figure 9. This method will make bees quickly gather towards food sources with high nectar value, making it difficult to ensure the diversity of food sources. To some extent, it will abandon many potential high-value food sources and prematurely fall into local optima. In order to improve the optimization effect of the algorithm and increase the diversity of food sources, an adaptive probability selection strategy is proposed. The strategy based on adaptive probability is based on the fitness probability of the population, and the selection probability is allocated according to changes in the number of iterations, in order to achieve the goal of selecting the probability of food sources even when the fitness value is small, as shown in Figure 10. The design of the adaptive process requires calculating the fitness probability of the food source, and the probability of the food source being selected is calculated by formula 4.1:

$$\mathcal{p}={\left|\mathcal{p}\mathrm{r}\mathrm{o}\mathrm{b}-{\displaystyle \frac{1}{\mathcal{i}\mathrm{t}\mathrm{e}\mathrm{r}}}\right|}^2$$

(4.1)

Among them, $\mathcal{p}\mathrm{r}\mathrm{o}\mathrm{b}$ is the probability of food source fitness, and $\mathcal{i}\mathrm{t}\mathrm{e}\mathrm{r}$ is the number of iterations.

By comparing the results of basic probability selection in Figure 9 with the adaptive probability selection in Figure 10, it can be seen that the adaptive probability can exhibit a parabolic shape, allowing for a high chance of being selected even when the initial fitness value is small. This not only improves the probability of food sources being selected when the initial fitness value is small before improvement, but also takes into account the probability of being selected when the fitness value is large, achieving the goal of diverse populations.

Based on the above design process, the detailed steps for solving the model using artificial bee colony algorithm are obtained:

Order 1: reads the coordinate point information of the delivery demand customer. The first allocation is the collection of demand customers who can receive delivery services from vehicles and cannot receive delivery from vehicles after departing from the distribution center. The number of demand customers is n, code each material demand customer, and then obtain the expected delivery time, demand, and vehicle and UAV load characteristics of each material demand customer.

Order 2: calculates the distance between each customer and forms a distance matrix.

Order 3: This includes: population size, number of leading bees, number of following bees, maximum number of leading bee searches, number of reconnaissance bees, number of iterations, number of customers, and number of UAVs.

Order 4: initializes the population based on the method used in the algorithm design to generate initial solutions, and calculates the fitness values of each solution based on the objective function and fitness function.

Order 5: runs the algorithm and begins the iteration process, repeating Order 6 to Order 12.

Order 6: leads the bee phase to execute all solutions once and continue searching for new solutions within the neighborhood of the solutions. And update the solution with a higher fitness value than the original solution based on the principle of greed.

Order 7: If the state of the solution in Order 7 has not been updated, the number of local searches corresponding to the solution is limited to limit+1. If the solution is updated, the corresponding search frequency variable limit=0.

Order 8: obtains the basic fitness probability based on the fitness values of all solutions, and then calculates the adaptive probability values of each solution being selected by the following bee.

Order 9: follows the bee to select all solutions according to the selected probability value calculated in Order 8. It continues to search for new solutions locally near the current solution, evaluates the fitness of the found new solutions, and still follows the greedy principle of selecting the optimal solution from the new solutions while saving the original solution. If the current solution has not been updated, the local search frequency variable corresponding to the current solution is also limited = limit + 1. If the solution is updated, limit = 0.

Order 10: If there is no updated solution within the Limit, discard the solution and switch to the reconnaissance bee stage to find a regenerated solution to replace it.

Order 11: records the current optimal solution.

Order 12: determines whether the global maximum cycle of iterations has been reached. If it has, it indicates the end of the algorithm. Otherwise, it goes to Order 6 and Order 11, where the recorded solution is the global optimal solution.

4.2. A Demonstrative Case

The ABC algorithm program in this chapter is written using Matlab 2016b software, and the PC parameter configuration for running the algorithm is Intel (R) Core i5 CPU T6600@2.23GHz 8GRAM/Windows 10 operating system. The basic instance data for validating the algorithm is sourced from the Solomon case study and literature dataset, and all results are obtained using the same configuration parameters. The various parameters of vehicles and UAVs, as well as the basic parameter settings of the artificial bee colony algorithm, are shown in Table 1. Using a 60-point case study as the experimental data for the customer set of emergency material distribution demand, the solution is studied, as shown in Table 2. Cluster the customer set for UAV delivery, as shown in Figure 11. The UAV delivery customers are divided into 3 categories and 2 independent delivery points. The maximum number of iterations and Limit values within the algorithm are set to 1500 and 100, respectively. The artificial bee colony algorithm designed in this paper and commonly used genetic algorithms are compared to solve for the optimal delivery paths of vehicles and UAVs, as shown in Figure 4.9-4.20. The paths of each UAV in Figure 12, Figure 13, Figure 14 and Figure 15 are represented in different colors.

4.3. Discussions

The comparison between artificial bee colony algorithm and genetic algorithm for optimizing vehicle and UAV paths shows that artificial bee colony algorithm reduces algorithm running time by 41.7% and total delivery time by 62.4%. Verified the effectiveness of the artificial bee colony algorithm. From Figures 4.18, 4.19, and 4.20, it can be seen that the convergence of the artificial bee colony algorithm in solving the examples in this paper is better than that of the genetic algorithm, and the quality of solving the examples in this paper is better.

5. Improved Artificial Bee Colony

In basic artificial bee colony algorithms, optimization is carried out based on the foraging behavior of leading bees, following bees, and scouting bees [63,64]. In the early stage, the local optimization effect is obvious and the convergence speed is fast, but in the later stage, it is easy to fall into local optima and it is difficult to obtain a global optimal solution. Therefore, based on the artificial bee colony algorithm, a simulated annealing strategy is introduced into the algorithm.

(1) Simulated annealing strategy design

The percentage difference between the new route and the current memory route:

$\mathrm{D}\mathrm{e}\mathrm{l}\mathrm{t}\mathrm{a}=\left(\mathrm{F}\mathrm{i}\mathrm{t}\right(\mathrm{n}\mathrm{e}\mathrm{w})-\mathrm{F}\mathrm{i}\mathrm{t}(\mathrm{c}\mathrm{u}\mathrm{r}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{t}\left)\right)/\mathrm{F}\mathrm{i}\mathrm{t}\left(\mathrm{c}\mathrm{u}\mathrm{r}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{t}\right)$Probability calculation:$\mathrm{P}=\mathrm{e}\mathrm{x}\mathrm{p}(-\mathrm{d}\mathrm{e}\mathrm{l}\mathrm{t}\mathrm{a}1/\mathrm{T})$Temperature calculation: $\mathrm{T}=\mathsf{\eta }\mathrm{T}0$(2) This section focuses on improving the local search process of the basic artificial bee colony algorithm. When following the bee to select a neighborhood, the preservation of the solution should be combined with the Metropolis rule of simulated annealing. When $\mathrm{n}\mathrm{e}\mathrm{w}-\mathrm{c}\mathrm{u}\mathrm{r}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{t}<0$, the current solution is updated. When $\mathrm{n}\mathrm{e}\mathrm{w}-\mathrm{c}\mathrm{u}\mathrm{r}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{t}>0$, concepts such as “probability” and “cooling temperature” are proposed to determine acceptable solutions. The algorithm improvement is shown in Figure 16:

(3) Detailed steps and process for improving the algorithm:

Process 1: Initial parameter setting: Maximum number of iterations after reaching $\mathrm{m}\mathrm{a}\mathrm{x}\mathrm{i}\mathrm{t}\mathrm{e}\mathrm{r}$, annealing coefficient$\mathsf{\eta }$， Iteration to the current $\mathrm{C}\mathrm{u}\mathrm{r}\mathrm{r}\mathrm{i}\mathrm{t}\mathrm{e}\mathrm{r}$, set the initial temperature value $\mathrm{T}=\mathrm{T}0$.

Process 2: Lead the bee to perform a random neighbourhood search operation on the current honey source location $\mathrm{X}$, obtain the new honey source $\mathrm{X}\mathrm{\text{'}}$, and calculate the objective function value $\mathrm{F}\left(\mathrm{X}\mathrm{\text{'}}\right)$ for that location. Accept the position by judging the size of the function value. If $\mathrm{F}\left(\mathrm{X}\mathrm{\text{'}}\right)$ is less than $\mathrm{F}\left(\mathrm{X}\right)$, the new honey source position $\mathrm{X}\mathrm{\text{'}}$ is adopted; Otherwise, abandon the honey source at that location.

Process 3: Follow the bee to select the hired bee through adaptive probability and search in the vicinity of the honey source where the selected hired bee is located. When a new honey source $\mathrm{X}\mathrm{\text{'}}\mathrm{\text{'}}$ is found, calculate the objective function value $\mathrm{F}\left(\mathrm{X}\mathrm{\text{'}}\mathrm{\text{'}}\right)$ for that location. Accept the position by judging the size of the function value. If $\mathrm{F}\left(\mathrm{X}\mathrm{\text{'}}\mathrm{\text{'}}\right)$ is less than $\mathrm{F}\left(\mathrm{X}\right)$, the new honey source position $\mathrm{X}\mathrm{\text{'}}\mathrm{\text{'}}$ is adopted; Otherwise, according to the designed annealing strategy, the probability $\mathrm{P}$ will accept the honey source at the new location.

Process 4: Reconnaissance bees follow the original criteria of the algorithm to determine and search for the honey source.

Process 5: Update $\mathrm{T}$ according to the temperature change method with attenuation coefficient, and calculate using the temperature calculation formula.

Process 6: Repeat Process 2, Process 3, and Process 4 until the maximum number of iterations set. In Process 3, following the simulated annealing strategy embedded with bees improves the execution efficiency of the artificial basic bee colony algorithm on the one hand; At the same time, by accepting probability $\mathrm{P}$, it is easy to know that in the early stage, if the temperature $\mathrm{T}$ is set to be high, the algorithm will have a higher acceptance probability for slightly worse solutions, making it easier for the algorithm to jump out of local optima. In this case, the level of development ability of the bee colony is very high. At the later stage of the algorithm, as the temperature drops to a relatively low value, the probability of accepting non optimal solutions gradually decreases, which better ensures the local search of the bee colony at low temperatures and can better find the optimal solution.

5.1. Results Analysis

This section uses the same example and uses the improved bee colony algorithm to optimize the paths for vehicle and UAV delivery to customers, as shown in Figures 17, 19, 22, and 24. The indicators for solving vehicle paths using the improved bee colony algorithm and the classical artificial bee colony algorithm are compared, as shown in Table 4. At the same time, provide a comparison chart of the iterative optimization curves of the algorithm, as shown in Figures 18, 19, 21, 23, and 24.

Figure 17. Improved Bee Colony Algorithm for Path Optimization

Figure 17. Improved Bee Colony Algorithm for Path Optimization

Figure 18. Comparison of iterative optimization curves for solving vehicle paths

Figure 18. Comparison of iterative optimization curves for solving vehicle paths

Figure 19. Solving path optimization for the class1 of UAV delivery customers using the improved the bee colony algorithm

Figure 19. Solving path optimization for the class1 of UAV delivery customers using the improved the bee colony algorithm

Figure 20. Comparison chart of iterative optimization curves between improved bee colony algorithm and artificial bee colony algorithm for solving the class 1 of UAV customer

Figure 20. Comparison chart of iterative optimization curves between improved bee colony algorithm and artificial bee colony algorithm for solving the class 1 of UAV customer

Figure 21. Improved bee colony algorithm to solve the path optimization for the class 2 of unmanned aerial vehicle delivery customers

Figure 21. Improved bee colony algorithm to solve the path optimization for the class 2 of unmanned aerial vehicle delivery customers

Figure 22. Comparison of iterative optimization curves between improved bee colony algorithm and artificial bee colony algorithm for solving the class 2 of UAV customers

Figure 22. Comparison of iterative optimization curves between improved bee colony algorithm and artificial bee colony algorithm for solving the class 2 of UAV customers

Figure 23. Improved bee colony algorithm to solve the path optimization for the class 3 of unmanned aerial vehicle delivery customers

Figure 23. Improved bee colony algorithm to solve the path optimization for the class 3 of unmanned aerial vehicle delivery customers

Figure 24. Comparison chart of iterative optimization curves between improved bee colony algorithm and artificial bee colony algorithm for solving the UAV customers of class 3

Figure 24. Comparison chart of iterative optimization curves between improved bee colony algorithm and artificial bee colony algorithm for solving the UAV customers of class 3

Based on the above figures, it is evident that as the number of iterations increases, the flight time decreases gradually. This result proves that the enhanced algorithm proposed in this article effectively optimizes the path and reduces the time for vehicle-UAV integrated delivery. Moreover, our model is cost-effective, making it suitable for real-time applications.

The results of improving the bee colony algorithm compared to artificial bee colony algorithm for solving vehicles and UAVs are shown in Table 5:

By comparing and analyzing the results in Table 5.2, the improved bee colony algorithm using simulated annealing strategy reduces the number of optimization iterations by 47% compared to the artificial bee colony algorithm, and also reduces the running time by 45%. This indicates that the improved bee colony algorithm can obtain the optimal solution with fewer iterations and shorter running time. The total path time obtained by solving also decreased by 14.7%, indicating that the improved bee colony algorithm can effectively jump out of the local optimal solution in the mid-term, reflecting the superiority of the improved bee colony algorithm.

6. Conclusions and Discussions

6.1. Conclusions

This study established a mathematical model for distributing emergency supplies using vehicles and UAVs and developed a solution model based on the Artificial Bee Colony (ABC) algorithm. An improved ABC algorithm, incorporating a simulated annealing strategy, was introduced to tackle large-scale optimization problems effectively. The research explored various delivery methods and operational modes for the integrated distribution of Vehicle-UAV, focusing on optimizing TSP (Traveling Salesman Problem) and VRP (Vehicle Routing Problem) paths. The goal was to establish a robust path model and design algorithms to solve it, ensuring time-optimal delivery paths that satisfy customer demands.The effectiveness of the classic ABC algorithm was validated against genetic algorithms, demonstrating superior performance in terms of average running time, iteration optimization times, and solution quality. Despite the ABC algorithm's fast early convergence, it showed limitations in local search capabilities during later stages. To address this, a simulated annealing strategy was introduced, leading to an improved bee colony algorithm. In tests with 60 customer examples, the improved algorithm achieved a 14.7% reduction in total path time, effectively escaping local optima during later stages. Compared to the classic ABC algorithm, it reduced optimization iterations by 47% and algorithm running time by 45%, highlighting the superiority of the improved bee colony algorithm.

6.2. Future Work

The article has contributed to relaxing the constraints of UAV delivery and has proposed a Vehicle-UAV integrated delivery model for the emergency supplies during crisis. An artificial bee colony algorithm was developed to solve the model, demonstrating fast convergence speed in the early stage and verifying its feasibility. However, due to the tendency to fall into local optima in the middle and later stages, an improved bee colony algorithm was proposed, aiming to enhance the model-solving methods. Nevertheless, there are still several other aspects that need to be addressed in future research endeavors. New material demands for emergency supplies might emerge during the distribution process. In response to the diverse scenarios, dynamic distribution material plans and dynamic coordination distribution plans involving multiple emergency material warehouses should be established. Additionally, the docking operation of UAVs should be made more flexible and active. Enabling the UAVs to cease loading vehicles at any time in the most time-saving situation and enabling them to park other nearby delivery vehicles as needed would contribute to achieving timely material delivery and optimizing flight delivery time. In the future, further research can explore the use of the artificial bee colony algorithm to develop initial solutions and subsequently select more suitable hybrid algorithms to optimize the path problem of UAVs, with the goal of achieving the optimal solution more efficiently.