Drone-Based Emergent Distribution of Packages to an Island from a Land Base

Abstract

An island logistics system is vulnerable in emergency conditions and even isolated from land logistics. Drone-based distribution is an emerging solution investigated in this study to transport packages from a land base to the islands. Considering the drone costs, drone landing platforms in islands, and incorporation into the island ground distribution system, this study categorizes the direct, point-to-point, and cyclic bi-stage distribution modes: in the direct mode, the packages are distributed from the drone base station to the customers directly by drones; in the point-to-point mode, the packages are transported to the drone landing platform and then distributed to the customers independently; in the cyclic mode, the packages are distributed from a drone landing platform by a closed route. The modes are formulated, and evaluation metrics and solution methods are developed. In the experiments based on an island case, the models and solution methods are demonstrated, compared, and analyzed. The cyclic bi-stage distribution mode can improve drone flying distance by 50%, and an iterative heuristic algorithm can further improve drone flying distance by 27.8%, and the ground costs by 3.16%, average for the settings of twenty to sixty customers and two to four drone landing platforms. Based on the modeling and experimental studies, managerial implications and possible extensions are discussed.

1. Introduction

With the increasing importance and frequent economic activities in island regions, the planning and construction of an island’s logistics system have entered a period of rapid development. In the face of challenges such as bad weather and road disruption, how to transport emergency medical products and urgent living materials timely and safely to the island customers, and how to better control the total cost of the system while improving efficiency and safety. In recent years, the rapid development of drone and artificial intelligence technologies has made drones widely tested and even used in various fields [1]. The advantages of drones, such as small body, fast speed, and freedom from terrain conditions, make it possible to replace traditional logistics methods to complete tasks more efficiently and effectively in some fields, e.g., emergency medical medicine delivery in disaster areas, agricultural pesticide spraying, field monitoring, power line patrol inspection, pastoral logistics distribution, etc. Many applications show that the research of drones in relevant fields is of great practical significance and commercial value [2].

In emergency logistics, the environment is complex and the demands are unpredictable, resulting in the mismatching among demands, supplies, and predictions. As an emerging distribution method, using drones can shorten the delivery time and deliver to places that cannot be reached by traditional transportation methods such as islands, and isolated or remote customers [1]. Drone-based distribution system incurs some advantages in emergency management: It can overcome various obstacles, quickly respond to emergency demands, and achieve accurate and efficient security; it has low cost and strong emergency capability, and can realize uninterrupted operations; it can adapt to challenging tasks, and achieve zero casualty rate in the distribution process. Distribution is a crucial link in emergency logistics and an important link that affects the response speed of the entire emergency system. Medical transport using drones can be used in urgent situations, where the main variable that has an impact on the success of life and health saving is the breaking of barriers to reaching difficult-to-reach places. In the context of the spread of the SARS-CoV-2 virus, drones may be used to provide diagnostic screening tests [3], medicinal products, and septic materials, transport of samples [4] of biological material, as well as information campaigns on how to deal with an epidemic, quarantine, or isolation at home [5]. With the rapid development of drone-related technologies, the applications of drones have been greatly improved, and at the same time, more requirements and implementations have emerged. The applications of 5G, the Internet of Things, GPS, artificial intelligence, and other technologies, as well as the deployment of communication infrastructures, have brought broad application prospects for drones. With the exploration of some leading retailing and distribution companies, e.g., amazon.com, JD.com, and sf-express.com, drone-based logistics has been promoted and applied in practice, gradually entering people’s vision. Especially, using drones can better solve the uncertainty and timeliness of emergency logistics with its advantages.

Islands usually lack relatively developed road networks, while their transportation means are limited to roads [6]. Moreover, some islands are far away from the mainland and are scattered, e.g., the eastern islands of Zhoushan in China, and the waterborne transportation routes between the islands and reefs are complicated. It is difficult for the islands to achieve their self-sufficient production and living needs. Due to the particularity and complexity of the geographical locations, an island is generally a closed system. When dealing with all kinds of disasters, an island has a prominent ecological and economic vulnerability, showing the characteristics of vulnerability, low self-regulation ability, and prominent disaster chain phenomenon. When the roads are damaged and impassable, the traditional truck and manual delivery methods are inefficient and difficult to be recovered timely, and the traditional logistics services are inaccessible. Rapid technological developments in autonomous unmanned aerial vehicles (UAV or drones) and evolving legislation may soon open the way for their large-scale implementation in the last-mile delivery of products [7]. Drone-based delivery can play the advantage of ignoring the terrain, shortening the delivery time, and improving efficiency, and can be effectively used for post-disaster relief and emergency distribution. The drone-related technologies have also been made achievements and gradually improved to ensure delivery capabilities and efficiency [8]. The application of drones will solve the bottleneck problems of high timeliness, high security, and high accuracy in the current island emergency logistics distribution.

This paper mainly studied the application of drones in emergency distribution in islands. Although drones have significant advantages in logistics distribution, there are still many limitations, e.g., limited load capabilities, energy and charging problems, landing safety, and landing platforms. The drone-based island emergency logistics system studied in this paper consists of three echelons of components: the onshore drone base station, the island mediate distribution depots, and the island terminal customers. Since the landing of drones often requires a landing platform to ensure landing safety, and each landing platform has additional costs, this paper considers two distribution situations: direct and bi-stage distribution. In the direct mode, drones directly distribute packages to the customers from the offshore drone base station. In the bi-stage mode, drones distribute the packages from the drone base station to the island depots with landing platforms first, and then island riders distribute the packages to the customers by ground transportation. Riders are persons taking motorbikes or mopeds for last-mile distribution. Riders are generally familiar with the roads and traffic conditions. The bi-stage mode further consists of two sub-modes, bi-stage point-to-point, and cyclic distribution modes, where the cyclic mode optimizes the ground distribution routes compared to the point-to-point mode by using models and solution algorithms for the general traveling salesman problem (TSP). The three modes and their evaluation metrics are formulated, and solution methods are developed.

The rest sections are organized as follows. Section 2 reviews the related studies on drone applications in emergency management and drone routing problems. Moreover, the incremental contributions to literature are elucidated. Then, we investigate the drone-based island emergency distribution problem, including the three distribution modes and evaluation metrics. In Section 4, the three modes are formulated, and Section 5 develops an iterative heuristic algorithm to solve the cyclic bi-stage distribution problem. We conducted a series of numerical experiments to examine the proposed models and algorithms in Section 6. Finally, we conclude the study in Section 7.

2. Related Studies

2.1. Drones in Emergency Management

With the continuous impact of climate change and human activities, a series of major risk events with sudden, random, and uncontrollable characteristics have erupted around the world. The new characteristics of chain, cluster, and multi-hazard superposition of disasters have brought huge challenges to the emergency management work of various governments, resulting in incalculable economic and social losses [9].

The key to emergency relief is how to deploy emergency supplies quickly and efficiently to the disaster-affected areas in the case of bad environment and poor road network conditions, minimize the disaster, and ensure the safety of life and property of the people in the disaster area. Drones, as one of the key technologies, have the advantages of high mobility, flexible deployment, and good line-of-sight links. In recent years, with the development of science and technology, drones have been widely applied in emergency scenarios [10].

Table 1. Pioneering studies on emergency management using drones.

Study	Research Problem	Method	Scenario
[11]	- A multi-period facility location. - Uncertainties in failure.- Maximize coverage with reliability constraints.	MILP	Emergency
[21]	- Air-ground cooperative emergency communication network based on drones. - Trajectory optimization of emergency drone base stations.	ML	Network
[15]	- Synchronized truck–drone delivery. - Priority of different disaster-stricken sites.	IP+DP	Disaster
[8]	- Individual attitudes, perceptions, and intentions for drone adoption. - Survey data of 146 mountain rescuers. - Intention to use drones is driven by performance gains and facilitating conditions.	RA	Disaster
[16]	- Consider disaster or enemy in emergency management. - VRP with trucks, drones, and random attacks.	ALNS	Disaster
[13]	- Use existing drone infrastructure. - Establish a drone-enabled backup transport system. - Policy on public-public and public–private partnerships.	A	Emergency
[26]	- Recovery of post-disaster wireless communication. - Cooperating drones for downlink transmission over rescue vehicles on the ground.	H	Network
[17]	- Cope drone facilities effectively. - Schedule drones in humanitarian logistics. - Uncertainties of drone operating conditions.	MH	Disaster
[14]	- Drone-based wireless network. - A set-covering problem with a realistic coverage radius.	MH	Network
[22]	- Identify areas with a high density of low mobility or stationary users. - Optimize drone base stations and user assignments.	ML+MOP	Network
[18]	- Locate depots considering inventories, service regions, and stochastic demands. - Minimize the overall system cost.	H	Disaster
[19]	- Routing drones through sampling locations. - Consider locations’ priorities and spatial correlations.	ALNS	Disaster
[23]	- Design drone trajectory and resource allocation to maximize average throughputs considering co-channel interference and completion time.	BP+PSO	Network
[14]	- A location–allocation problem with coverage distance and capacity of drones. - consider demand uncertainties.	MH	Emergency
[24]	- Optimize delays for multi-drone emergency networks. - Consider users’ priorities and communication delays.	H	Network
[25]	- Integrate drones and ground mobile devices. - Optimize transitions from the control center to the drones.	H	Network
[20]	- Locate drone landing sites in complicated environments. - Consider terrain uncertainty, safety, fuel consumption, and path planning.	H	Disaster
[12]	- Schedule bloodmobiles and drones to collect blood from donors. - Consider uncertain blood demands and donors.	MILP	Emergency

Note: A = mathematical analytics; ALNS = adaptive large neighborhood search; BP = bilevel program; DP = dynamic programming; dynamic programming; H = general heuristics; IP = integer program; MH = math-heuristics MILP = mixed-integer linear program; ML = machine learning including deep learning and reinforcement learning; MOP = multi-objective optimization; PSO = particle swarm optimization; RA = regression analysis.

summarizes studies on drone-based emergency management from three aspects: the problem features, methods, and application scenarios. Possibly due to the complexities of the studied problems, most studies in Table 1 did not formulate regular models, e.g., IP and MILP [11,12]. As a result, various heuristics and analytical methods were developed. Summarizing the reviewed papers, we can identify and classify the scenarios into three: emergency [11,12,13,14], disaster [15,16,17,18,19,20], and network [14,21,22,23,24,25]. The disaster scenario mainly concerns post-disaster delivery and various supply and logistics aspects. The network scenario includes two aspects of communication networks, network recovery after disasters or destruction by using drones temporally and movable drone-based communication networks for military or special objectives. The emergency scenario represents other scenarios, mainly emergency medical package delivery in Table 1. We can also identify two key features from the table, uncertainty, and synchronization, which may make the solution methods challenging.

2.2. Drone-Based Routing Problems

Drones have been applied in various scenarios, depending on drone performances, including hovering stabilities, loading capacities, and endurance mileages. They are primarily determined by technological innovations in drone electricity, machinery, and control systems of drones.

Considering a fleet of drones, the performances, especially efficiency-related measures, mainly depend on scheduling, routing, coordinating, and management strategies. In this study, we are mainly concerned with transportation and logistics-related scenarios. To satisfy the logistics missions, the drones should be affected organized, scheduled, and managed. Especially, to accomplish many tasks, the drones should coordinate with other devices and facilities to strengthen the drones’ strength and overcome the shortages. In many studies, drones will cooperate with trucks, where trucks undertake branch transportation and act as drone landing platforms and charging stations.

In

Table 2. Pioneering studies on drone-based routing problems.

Study	Problem Features	Devices #	Model *	Algorithm **
[29]	- Extend the crossing postman problem to visit different shapes of dimensional elements. - Drones visit geographical elements to deliver goods or services.	D	MILP	MH
[27]	- Drones monitor a set of areas with different accuracy requirements. - A VRP considering flying heights, which impacts the accuracy level and service time.	D	IP	Tabu
[30]	- A multi-depot VRP with separation distance constraints.	D	MILP	MH+Tabu
[43]	- Routing drones in patrolling missions.	D	SP	H
[31]	- A drone travels with a truck, takes off from its stop to serve customers, and lands at a service hub to travel with another truck if the flying range and loading capacity limitations are satisfied.	D+T	MILP	MH
[41]	- Consider a tandem between a truck and $k$ drones. - Each drone loads one or more packages to deliver to customers. - Each drone may return to the truck to swap/recharge batteries.	D+T	-	H
[32]	- Use drones to gather information from targets in military missions. - Drones hold limited sensors considering weights and flying time. - Assign sensors to drones and maximize intelligence gain.	D	MILP	MH
[33]	- Consider the mothership and drone routing problem. - A drone launches from the mothership, and then return to refuel. - The drone has a limited range of time units.	D+S	MILP	MH
[28]	- Consider the drone routing problem with recharging stops.	D	IP	H+MH
[34]	- Coordinate one mothership with one drone. - Minimize the overall weighted distance traveled while satisfying percentages of visits to targets.	D+S	MILP	MH
[44]	- A large drone carries multiple small drones to distribution regions; - Large drone does not directly deliver parcels, but rather launches small drones to deliver parcels. - Allows each small drone to deliver multiple parcels in a flight. - Considering energy consumption.	D+D	DP	VNS
[39]	- Synchronize drones and delivery trucks. - Trucks can work as mobile launching and retrieval sites. - A VRP with Time Windows and Synchronized Drones.	D+T	MOP	ACO+NSGA-II
[35]	- Consider multi-trip, recharge, and energy savings operations. - Drones’ energy consumption is modeled as a nonlinear function of payload and travel distance.	D	MILP	MH
[36]	- A drone travels to multiple locations to perform specified observation tasks and rendezvous periodically with the truck to swap its batteries.	D+T	MILP	NS
[37]	- Drones launch from or land on a moving truck without stopping. - Synchronize trucks and drones on arcs. - Trucks dispatch and retrieve drones at suitable moving locations on arcs of truck routes.	D+T	MILP	ALNS
[38]	- A truck and drone routing problem. - A drone can be launched from the truck at an intermediate depot once (single-trip drone) or several times (multiple-trip drone).	D+T	MILP	GRASP
[42]	- Consider a synchronized truck–drone operation by allowing multiple drones to fly from a truck, serve one or multiple customers, and return to the same truck for a battery swap and package retrieval.	D+T	MIP	LNS
[26]	- Introduce a new truck–drone tandem that allows the truck to stop at non-customer locations.	D+T	MILP	SA+VNS
[40]	- A truck collaborates with drones to perform parcel deliveries, and each customer can be served earlier and later than the required time with a given tolerance. - Optimize total distribution cost and overall customer satisfaction.	D+T	MOP	H

Note: #: D = drone, T = truck, S = mothership; *: DP = dynamic programming, IP = integer program, MILP = mixed-integer linear program, MIP = mixed-integer program, MOP = multi-objective optimization problem, SP = stochastic program; **: ACO = Ant Colony Optimization algorithm, ALNS = variable neighborhood descent algorithm, GRASP = greedy randomized adaptive search procedure, H = Heuristics, LNS = Large Neighborhood Search, MH = Math-heuristics, NS = Neighborhood search algorithm, NSGA-II = Non-dominated Sorting Genetic Algorithm, SA = simulated annealing, Tabu = Tabu search algorithm, VNS = variable neighborhood search.

, we summarize some pioneering studies in four columns. The drone routing problems are featured in various aspects, mainly including application scenarios, new features comparing existing popular routing problems, and methodological challenges. These features may further be reflected in the other three columns in the table. First, in logistics, constrained by load weights and endurance distances, drones generally coordinate with trucks or other kinds of motherships. Second, the drone routing problems are generally formatted as IP [27,28] and MILP [26,29,30,31,32,33,34,35,36,37,38] models. Considering various performance metrics, e.g., cost, distance, energy consumption, and consumer satisfaction, the drone routing problems can be formulated and analyzed by multi-objective optimization models [39,40]. Fourth, the drone routing problems generally couple drones and other devices, e.g., trucks [26,31,36,37,38,39,40,41,42] and motherships [33,34], which makes the models challenging in algorithm development. Various algorithms, including meta- and math-heuristics [29,30,31,32,33,34,35], have been studied considering the problem and model features.

2.3. Incremental Contributions to Literature

Based on the background study in Section 1, and the two streams of studies on drones in Section 2.1 and Section 2.2, this study contributes to the literature in the following three aspects. First, an island is a special scenario suitable for drone-based applications and the island logistics systems are practically vulnerable, while drones can help build connections with the offshore logistics systems. Second, as studied in Section 2.1 and Section 2.2, drones in logistics are generally applied to “last-mile” delivery, while this study takes drones to serve the first-stage distribution from the offshore base station to the islands. Third, comparing the studies in Table 1 and Table 2, this study developed three drone-based distribution modes considering the balances between drone-based and ground distributions that are medicated by drone landing platforms. The above distinct features activated this study.

3. Problem Statement

3.1. The Problem

As depicted in Figure 1, eight customers on an island are distributed by the nearby land freight terminal. First, the packages are distributed from the land terminal to the island terminal by cargo ships, generally motorboats, for an emergency. Then, the packages are distributed to the island customers by riders or other means. However, when it is foggy, windy, or there are other climate conditions not suitable for berthing and shipping, the emergency demands from the customers cannot be met. Due to the COVID-19 pandemic, vaccines, drugs, and medicine packages are urgent. Considering such conditions, it is beneficial to apply drones for island distribution.

In this study, we consider a drone base station at the original freight terminal, whose location is denoted by $O$. The coordinate is denoted by $\left(X^o,Y^o\right)$. A set of customers are denoted by $I$, indexed by $i,j$ generally. The distance between the customers $i,j$ are denoted by $C_{ij}$, and the straight distance from the station to the customer $i$ is $C_i^o$.

3.2. Three Drone-Based Distribution Modes

Drones can be used in the scenario described in Section 3.1. In the following, three solutions are developed considering drones in the island distribution scenario.

First, in the “direct distribution” mode, the drones take and send the packages to the customers directly, as depicted in Figure 2a. Here, the scenario of island distribution is used as the same described in Figure 1. As seen from these figures, the freight terminal in Figure 1 is replaced by a drone base station in Figure 2. The transportation means from the offshore to the island use drones other than ships. In this study, we consider the simplest direct distribution mode, where a drone serves a customer one time. In such a mode, the island distribution system can use drones of lower costs comparatively. However, each customer must install a drone landing platform to accept the packages from the drones.

The drone landing platform is similar to an express cabinet with a platform on top of it to accept packages dropped by drones, as depicted in Figure 3. In the vertical view of the platform, it consists of an automatic door that will be opened when accepting and sensing a drone to drop packages onto it. There are graphics on the platform, and there are sensors attached to it. They can help drones to identify the platform (DLP) and locate the door to open for accepting the packages.

A DLP can be shared within a residential quarter or with some close neighbors. It also can be installed as an express cabinet at the same time, while it should be suitable for drones to locate and deliver packages. As a result, it will be costly and will occupy space.

Second, as depicted in Figure 2b, an island distribution system can consist of two stages. In the first stage, the packages are distributed from the original terminal to some island DLPs by drones. In the second stage, the riders distributed the packages from the DLPs to the customers one by one, or the customers can go to the DLP and accept the packages from the cabinets. Because the packages are distributed from the DLPs to the customers one by one, we entitle this solution as a “Point-to-point bi-stage solution”.

Based on the second mode—“Point-to-point bi-stage distribution”, a specialized rider can be managed to distribute the packages to the customers in a cycle—visiting the customers sequentially for a DLP and finally returning to the DLP. This problem is famous as the TSP. This solution is beneficial when the batch of distribution consists of several packages to different customers.

Three notations, $dir, ptp, cyc$, are used to represent the three drone-based distribution modes.

3.3. Distribution Mode Evaluation Metrics

As studied in Section 3.1 and Section 3.2, three solutions incur different strategies for using drones, riders, and DLPs. To evaluate these solutions, four metrics are developed as follows.

(1) Drone cost ($z^{dd}$). In the “direct distribution” solution, the distribution system can use drones with smaller load capacities because a drone is used to distribute packages to a single customer. In the “Point-to-point” or “Cyclic” “bi-stage distribution” solutions, a drone is used to distribute packages to a group of customers. As a result, drones with larger load capacities should be considered. It is a capacity decision-making problem and will not be studied mathematically or experimentally in this study. The comparison study contributes to the knowledge of different solutions.

(2) Drone flying distance ($z^{df}$). In the three solutions studied above, drones are used to transport packages from the land to the island. The drone flying cost is primarily determined by the drones’ flying distances and packages’ weights, while affected by the drones’ original and depreciation costs, or the rental protocols. In the following study, it is presumed that the drone flying cost depends on the flying distances and packages.

(3) Drone landing platform ($z^{DLP}$). DLP is an important and costly facility in a drone-based distribution system. Its cost depends on its size and capacity, installed systems, and even the space occupied by it. Its capacity depends on the delivered packages and the frequencies, and the riders’ schedules.

(4) Ground distribution ($z^{gd}$). In the two “bi-stage distribution” solutions, riders will distribute the packages from the DLPs to the customers. The “cyclic” solution will incur fewer traveling distances than the “point-to-point” solution. The ground distribution cost consists of fixed and variable costs. Generally, the system should pay fixed basic salaries to the riders and fixed costs of distribution facilities and devices, while the traveling time and distances construct variable costs. In the following study, the variable cost of traveling distances is considered.

Table 3. Comparisons among three drone-based distribution solutions.

Distribution Mode	Drone Cost ( ${\mathit{z}}^{\mathit{d}\mathit{d}}$)	Fly Distance ( ${\mathit{z}}^{\mathit{d}\mathit{f}}$)	DLPs ( ${\mathit{z}}^{\mathit{D}\mathit{L}\mathit{P}}$)	Ground Distribution ( ${\mathit{z}}^{\mathit{g}\mathit{d}}$)
Direct ($dir$)	Low	High	Many	None
Point-to-point ($ptp$)	High	Low	A few	High
Cyclic ($cyc$)	High	Low	A few	Low

summarizes the comparisons in four metrics among three distribution solutions. Here, the following words are used to express the comparative relations: Low, High, Many, and A few. They indicate the differences among the modes and do not indicate concrete value bounds.

4. Formulation

In the following, the models of the three drone-based island distribution modes are formulated, as well as the four metrics for evaluating the modes.

(1)	Names
$dir$	The direct distribution mode
$ptp$	The point-to-point bi-stage distribution mode
$cyc$	The cyclic bi-stage distribution mode
$O$	The identity of the drone base station on the land side
(2)	Set
$I$	A set of customers, indexed by $i,j$
(3)	Data
$C_i^o$	The distance from the drone base station to the customer $i$
$C_{ij}$	The distance between the customers, $i,j$
$D_i$	The number of packages distributed to the customer $i$
(4)	Variable
$z^{dd}$	The cost of buying drones
$z^{df}$	The drone flying distances are weighted by the number of packages
$z^{DLP}$	The number of DLPs
$z^{gd}$	The riders’ traveling distances are weighted by the number of packages

4.1. Direct Distribution

As studied in Section 3, in the direct distribution mode, the packages are distributed from the drone base station to the customers directly. $\left[\mathrm{M}\left(dir\right)\right]$ computes three metrics, namely, $z^{df},z^{DLP},z^{gd}$, as defined in (1). It is denoted as a model, although it is not a formal mathematical program.

$$\left[\mathrm{M}\left(dir\right)\right]\left\{\begin{array}{c}{z}^{df}\left(dir\right)={\displaystyle \sum }_i{C}_i^o{D}_i\\ z^{DLP}\left(dir\right)=\left|I\right|\\ z^{gd}\left(dir\right)=0\end{array}\right.$$

(1)

Further, $\left[\mathrm{M}\left(dir\right)\right]$ can be represented by a function with inputs and outputs as (2).

$$\left(z^{df},z^{DLP},z^{gd}\right)\leftarrow f^{\left[\mathrm{M}\left(dir\right)\right]}\left(I,O,C^o,D\right).$$

(2)

4.2. Point-to-Point Bi-Stage Distribution

Aside from the notations defined in

Table 4. A dataset $I20P3$ used in the demonstration.

$$\mathit{I}$$	Latitude	Longitude	$${\mathit{D}}_{\mathit{i}}$$	$${\mathit{C}}_{\mathit{i}}^{\mathit{o}}$$	$$\mathit{I}$$	Latitude	Longitude	$${\mathit{D}}_{\mathit{i}}$$	$${\mathit{C}}_{\mathit{i}}^{\mathit{o}}$$
1	122.3933	29.99097	16	4.034	11	122.4026	30.02376	2	7.762
2	122.3977	30.01769	8	6.976	12	122.4048	30.02670	2	8.14
3	122.4042	30.02551	4	7.995	13	122.4191	30.00947	18	7.055
4	122.3951	29.99045	20	4.051	14	122.3976	30.00798	2	5.946
5	122.3971	30.01360	14	6.524	15	122.3974	30.00014	11	5.125
6	122.3911	30.01322	10	6.335	16	122.3961	30.01555	1	6.705
7	122.3981	30.00404	5	5.554	17	122.4232	30.00806	20	7.162
8	122.3972	30.00713	16	5.845	18	122.4023	30.02311	6	7.686
9	122.397	29.99425	8	4.512	19	122.3976	30.00798	1	5.946
10	122.3956	30.01369	14	6.492	20	122.3994	30.00887	10	6.099

, the binary variable $x_{ij}$ is introduced to represent the assignment of the customer $j$ to a candidate DLP at the customer $i$. When $x_{ii}=1$, a DLP is installed at the customer $i$. Then, a new variable $y_i$ is used to represent the packages distributed to the DLP installed at $i$. Then, a multi-objective model is devised as $\left[\mathrm{M}\left(ptp\right)\right]$.

$$\left[\mathrm{M}\left(ptp\right)\right]\min \left(z^{df},z^{DLP},z^{gd}\right)$$

where:

$$z^{df}={\displaystyle \sum }_i{C}_i^o{y}_i$$

(3)

$$z^{DLP}={\displaystyle \sum }_i{x}_{ii}$$

(4)

$$z^{gd}={\displaystyle \sum }_{ij}{x}_{ij}{D}_j{C}_{ij}$$

(5)

Subject to

$${\displaystyle \sum }_i{x}_{ij}=1,\forall j$$

(6)

$$x_{ij}\le x_{ii},\forall i,j$$

(7)

$$y_i={\displaystyle \sum }_j{x}_{ij}{D}_j,\forall i$$

(8)

$$x_{ij}\in \left\{0,1\right\},y_i\ge 0,\forall i,j$$

(9)

In Constraint (6), all customers should be serviced. In Constraint (7), only DLP can be assigned to a customer. In Constraint (8), the accepted packages at a DLP are computed. Finally, Constraint (9) defines the variable integrities and domains.

When the ground distribution cost should be optimized with the highest priority under a determined number of DLPs, the model can be transferred into a mixed-integer linear program (MILP), as formulated by $[\mathrm{M}(ptp,{\overline{f}}^{DLP})]$.

$$\left[\mathrm{M}\left(ptp,{\overline{z}}^{DLP}\right)\right]\min z^{gd}$$

Subject to

$$Constraints \left(\left(2\right)-\left(7\right)\right)$$

$$z^{DLP}={\overline{z}}^{DLP}$$

(10)

$\left[\mathrm{M}\left(ptp,{\overline{z}}^{DLP}\right)\right]$ is a single-objective linear program with a parameter ${\overline{z}}^{DLP}$, which is the number of DLPs. After solving the model, obtain $x,y$ and then $z^{df}\left(ptp\right)$ can be computed by (3). In the following sections, $\left[\mathrm{M}\left(ptp,{\overline{z}}^{DLP}\right)\right]$ can be used in the following form,

$$\left(z^{df},z^{DLP},z^{gd},x,y\right)\leftarrow f^{\left[\mathrm{M}\left(ptp,{\overline{z}}^{DLP}\right)\right]}\left(I,O,C^o,C_{ij},D;{\overline{z}}^{DLP}\right).$$

4.3. Cyclic Bi-Stage Distribution

In Section 4.2, the ${\overline{z}}^{DLP}$ DLPs are chosen by the model $\left[\mathrm{M}\left(ptp,{\overline{z}}^{DLP}\right)\right]$, as well as the assignments ($x$) of customers to the DLPs. So, the DLPs of the solution can be denoted by $I^{DLP}$; for each $d\in I^{DLP}$, the assigned customers are $I_d^{DLP}\subseteq I$. The above process can be interpreted as follows,

$$I^{DLP},\left\{I_d^{DLP}|\forall d\in I^{DLP}\right\}\leftarrow x.$$

In the following, a TSP model is used to find a sequence of the customers in $N$ with minimal traveling distance. Here, $N$ is indexed by $i,j$ default.

Two decision variables are introduced here. First, a binary variable $x_{ij}^c$ is used to represent sequential visits to the customers. If $x_{ij}^c=1$, a rider visits $i$ first and then adjacently visits $j$; otherwise, 0. Another positive variable $u_i$ is used to represent the visiting order of the customer $i$, which is also used to avoid visiting loops in a TSP solution.

$$\left[\mathrm{M}\left(cyc,N\right)\right]\min z^{gd}$$

where

$$z^{gd}={\displaystyle \sum }_{ij}{D}_i{C}_{ij}{x}_{ij}$$

(11)

Subject to

$${{\displaystyle \sum }}_{j\in N,i\ne j}{x}_{ij}^c=1, \forall i$$

(12)

$${{\displaystyle \sum }}_{i\in N,i\ne j}{x}_{ij}^c=1, \forall j$$

(13)

$$u_i-u_j+\left|N\right|·x_{ij}^c\le \left|N\right|-1, \forall i\in N,j\in N\backslash 0,i\ne j$$

(14)

$$1\le u_i\le \left|N\right|-1,\forall i\in N\backslash 0$$

(15)

$$x_{ij}^c\in \left\{0,1\right\}, \forall i,j\in N$$

(16)

In the model $\left[\mathrm{M}\left(cyc,N\right)\right]$, Constraint (11) further defines the objective to minimize the total traveling cost. For each customer, Constraints (12) and (13) ensure that each customer is visited once and only once. Constraint (14) is used to avoid loops in the final sequence determined by $x$ and $u$. Finally, Constraints (15) and (16) define the domain of $u$ and the integrity of $x$. The solving of $\left[\mathrm{M}\left(cyc,N\right)\right]$ can be expressed by $f^{\left[\mathrm{M}\left(cyc,N\right)\right]}$, as follows,

$$\left(z^{gd}\left(cyc,N\right),x^c,u\right)\leftarrow f^{\left[\mathrm{M}\left(cyc,N\right)\right]}\left(I,C_{\cdot }^o,C_{ij},D_{\cdot }\right).$$

Then, using $f^{\left[\mathrm{M}\left(cyc,N\right)\right]}$, by iterating the DLP $d$ and its corresponding serviced customers $I_d^{DLP}$, the total ground distribution distance can be computed by (17), while the $z^{df}\left(cyc\right),z^{DLP}\left(cyc\right)$ are taken from $z^{df}\left(ptp\right),z^{DLP}\left(ptp\right)$, not changing, as formulated in (18).

$$z^{gd}\left(cyc\right)\leftarrow {\displaystyle \sum }_{d\in I^{DLR}}{f}^{\left[\mathrm{M}\left(cyc,I_d^{DLP}\right)\right]}\left(I,C^o,C_{ij},D;I_d^{DLP}\right)$$

(17)

$$\left(z^{df}\left(cyc\right),z^{DLP}\left(cyc\right)\right)\leftarrow \left(z^{df}\left(ptp\right),z^{DLP}\left(ptp\right)\right)$$

(18)

The process above uses a bi-stage method to compute the “point-to-point” solution and then optimizes it to obtain a “cyclic” distribution solution. The whole process is denoted by $f^{cyc}$ as follows,

$$\left(z^{df},z^{DLP},z^{gd},x,y\right)\leftarrow f^{cyc}\left(I,O,C^o,C_{ij},D;{\overline{z}}^{DLP}\right).$$

The process $f^{cyc}$ can be improved by iterative heuristics, which is studied in Section 5.

5. Iterative Heuristics for Improving Cyclic Bi-Stage Distribution

The “$cyc$” mode aims at optimizing the DLPs and ground distribution at the same time, while $f^{cyc}$ (Section 4.3) optimizes them in two stages. Therefore, $f^{cyc}$ can be improved by iterating these two stages, as developed in Algorithm 1.

Algorithm 1	Iterative Heuristics for Improving Cyclic Mode (IH)
Input	$I,O,C^o,C,D$: the known data
	${\overline{z}}^{DLP}$: the number of DLPs to be installed.
Output	$z^{df},z^{gd},x,y,I^{DLP},I_{+}^{DLP}$: the metrics and bi-stage distribution network.
Variable	${\check{z}}^{df}=+\infty$, the initial drone flying distance;
	${\check{z}}^{gd}=+\infty$, the initial ground distribution distances.
Steps
Step 1	Choose DLPs in the point-to-point bi-stage distribution mode.
	$\left(z^{df}\left(ptp\right),z^{DLP}\left(ptp\right),z^{gd}\left(ptp\right),x,y\right)\leftarrow f^{\left[\mathrm{M}\left(ptp,{\overline{z}}^{DLP}\right)\right]}\left(I,O,C^o,C_{ij},D\right)$.
Step 2	Prepare cyclic distribution data
	Generate $I^{DLP}$, $I_{+}^{DLP}=\left\{I_d^{DLP}\subseteq I|d\in I^{DLP}\right\}$ from $x$.
Step 3	Optimize distribution routes in the cyclic bi-stage mode
Step 3.1	Compute $z^{gd}\left(cyc\right)$ by solving $\left[\mathrm{M}\left(cyc,I_d^{DLP}\right)\right]$ for each $d\in I^{DLP}$
	$z^{gd}\left(cyc\right)\leftarrow {\sum }_{d\in I^{DLR}}{f}^{\left[\mathrm{M}\left(cyc,I_d^{DLP}\right)\right]}\left(I,C^o,C_{ij},D\right)$.
Step 3.2	Use $z^{df}\left(ptp\right),z^{DLP}\left(ptp\right)$ as they are not changed.
	$\left(z^{df}\left(cyc\right),z^{DLP}\left(cyc\right)\right)\leftarrow \left(z^{df}\left(ptp\right),z^{DLR}\left(ptp\right)\right)$.
Step 3.3	Merge the results of $f^{cyc}$
	$\left(z^{df},z^{DLP},z^{gd},x,y\right)\leftarrow f^{cyc}\left(I,O,C^o,C_{ij},D\right)$.
Step 4	Terminal criteria
	If $z^{df}<{\check{z}}^{df} and z^{gd}<{\check{z}}^{gd}$ Then
	${\check{z}}^{df},{\check{z}}^{gd}\leftarrow z^{df}, z^{gd}$
	Else
	Go to Step 8
Step 5	Update the DLPs by using the closed customer in each cyclic route
	For $d\in I^{DLP}$
	$d^{min}=\arg \underset{i\in I_d^{DLR}}{\min }{C}_i^o$
	Replace $d$ in $I^{DLR}$ by $d^{min}$
Step 6	Update the assignments of customers to DLPs
	Set $I_d^{DLP}=\varnothing$ for each $d\in I^{DLP}$
	For each $i\in I$,
	$d=\arg \underset{j\in I^{DLR}}{\min }{C}_{ij}$.
	insert $i$ into $I_d^{DLR}$
Step 7	Go to Step 3
Step 8	Return $z^{df},z^{gd},x,y,I^{DLP},I_{+}^{DLP}$

Algorithm 1 will iteratively improve $\left(z^{df},z^{gd}\right)$. In Step 1, an initial solution is obtained by $\mathrm{M}\left(ptp,{\overline{z}}^{DLP}\right)$. Step 3 provides a holistic flow of solving $f^{cyc}$. Step 4 sets the termination criteria so that $\left(z^{df},z^{gd}\right)$ should all be improved. Step 5 optimizes the drone flying distances by choosing the closest customer to the drone base station as the DLP for each “cycle”. Step 6 updates the assignments of customers to each DLP, and goes to the iteration loop in Step 7.

The algorithm is denoted by $f^{\left[\mathrm{IH}\right]}$ as follows,

$$z^{df},z^{gd},x,y,I^{DLP},I_{+}^{DLP}\leftarrow f^{\left[\mathrm{IH}\right]}\left(I,O,C^o,C,D;{\overline{z}}^{DLP}\right).$$

6. Numerical Experiments

6.1. Parameter Estimation

As studied in Section 3, a China island (Mount Putuo) is chosen as the scenario of the study, where the drone base station is set in the Zhujiajian passenger and cargo transport terminal $\left(O\right)$. Mount Putuo is in Putuo District, Zhoushan City, Zhejiang Province, the southern edge of Hangzhou Bay, and the eastern sea area of Zhoushan Islands. It is a famous tourist resort in China and one of the four famous Buddhist mountains in China. From the real estate, important facilities, temples, and scenic spots, 287 locations are identified as customers $\left(I\right)$ in this study. The customers’ demands $\left(D\right)$ are estimated by considering the residents and tourists. The latitudes and longitudes of these customers are generated by using a map service (www.amap.com), which provides a callable Application Program Interface (API) to obtain the latitude and longitude from a given address after we obtain the addresses of the 287 communities. The distances from the drone base station to the customers are computed using geographical spherical straight distances $\left(C^o\right)$ directly, while the distances among the customers are 1.4 times straight distances $\left(C\right)$, considering the rugged mountain roads. The number of DLPs $\left({\overline{z}}^{DLP}\right)$ is set to 3 generally, which can be adjusted in the experiments.

6.2. Dataset Generation

The datasets are generated based on the parameters estimated above. The computing complexity is mainly determined by the number of customers and the DLPs. As a result, we use $InPm$ to represent the datasets, where $n$ is the number of customers, $n=\left|I\right|$, and $m$ is the number of DLPs. Here $n<287$ and $m={\overline{z}}^{DLP}$. Generally, $I20P3$ datasets are used for demonstration, while other $\left(n,m\right)$ configurations are used to study the problem features and algorithm performances.

6.3. Experiments and Results

We conducted three groups of experiments to study the three drone-based distribution modes and their solution methods in the following subsections. In the experiments, a personal computer with the CPU, Intel(R) Xeon(R) CPU E3-1535M v6 @ 3.10GHz 3.10 GHz, and 64 GB RAM. All the solution algorithms are implemented by Python 3.7 and the models, $\left[M\left(ptp\right)\right], \left[M\left(cyc,N\right)\right]$, are solved by Cplex 12.9 (https://www.ibm.com/products/ilog-cplex-optimization-studio, accessed on 19 March 2023).

6.3.1. Demonstration of Drone-Based Distribution Modes and Solution Methods

In Section 3, the three drone-based distribution modes ($dir,ptp,cyc$) are described conceptually, and further, the direct and point-to-point bi-stage distribution modes are formulated (as $\left[\mathrm{M}\left(dir\right)\right]$, $\left[\mathrm{M}\left(ptp,{\overline{z}}^{DLP}\right)\right]$) in Section 4. Considering the complexity of solving the cyclic bi-stage distribution model ($cyc$), an iterative heuristic algorithm ($f^{\left[\mathrm{IH}\right]}$) in Section 5 is developed to improve the solution method $f^{cyc}$ devised in Section 4. In Table 4, a dataset with 20 customers and setting two DLPs is generated to demonstrate these modes, models, and algorithms.

The results corresponding to the four solution methods are depicted in Figure 4. Intuitively, the bi-stage modes (Figure 4b–d) reduce the drones’ flying distances of the direct distribution mode (Figure 4a) apparently; the cyclic mode (Figure 4c,d) is competitive in reducing ground distribution distance compared to the point-to-point mode (Figure 4b); and the iterative heuristic algorithm can improve the results by changing the DLPs and the ground distribution routes (Figure 4c,d).

Table 5. Evaluating the solutions of four solution methods.

Solution Method	Drone Fly Distance ( ${\mathit{z}}^{\mathit{d}\mathit{f}}$/km)	Ground Distribution (${\mathit{z}}^{\mathit{g}\mathit{d}}$/km)	DLPs ( ${\overline{\mathit{z}}}^{\mathit{D}\mathit{L}\mathit{P}}$)
$\left[\mathrm{M}\left(dir\right)\right]$	1126.370	0.00	20
$\left[\mathrm{M}\left(ptp,{\overline{z}}^{DLP}=3\right)\right]$	1114.741	113.59	3
$f^{cyc}\left({\overline{z}}^{DLP}=3\right)$	1114.741	103.28	3
$f^{\left[\mathrm{IH}\right]}\left({\overline{z}}^{DLP}=3\right)$	1015.856	138.57	3

presents three criteria for the four solution methods. In the direct distribution mode, all customers will install DLPs, while three DLPs are set in the bi-stage modes. As a result, the “$dir$” mode incurs no ground distribution. The bi-stage mode can reduce the drone flying distance by $1.03\%$ ($\left(1126.370-1114.741\right)/1126.370$). The iterative heuristics can reduce $9.71\%$ ($\left(1126.370-1015.856\right)/1126.370$). The “$cyc$” mode can reduce $9.08\%$ ($\left(113.59-103.28\right)/113.59$) compared to the “$p2p$” mode.

6.3.2. Impacts of DLPs on Bi-Stage Distribution Modes

The DLPs are important facilities demanding great investment in drone-based bi-stage distribution modes. In

Table 6. Comparative results of solving $\left[\mathrm{M}\left(ptp\right)\right]$ and $f^{cyc}$.

DLPs	I	$${\mathit{z}}^{\mathit{d}\mathit{f}}\left(\mathit{p}\mathit{t}\mathit{p}\right)$$	$${\mathit{z}}^{\mathit{g}\mathit{d}}\left(\mathit{p}\mathit{t}\mathit{p}\right)$$	$$\mathbf{CT}\left(\mathit{p}\mathit{t}\mathit{p}\right)$$	$${\mathit{z}}^{\mathit{d}\mathit{f}}\left(\mathit{c}\mathit{y}\mathit{c}\right)$$	$${\mathit{z}}^{\mathit{g}\mathit{d}}\left(\mathit{c}\mathit{y}\mathit{c}\right)$$	$$\mathbf{CT}\left(\mathit{c}\mathit{y}\mathit{c}\right)$$	$${\mathit{z}}^{\mathit{g}\mathit{d}*}$$	$$\mathit{C}{\mathit{T}}^{*}$$
		km	km	s	km	km	s	%	%
2	20	746.26	331.20	0.09	746.26	147.91	0.11	55	22
2	25	959.07	439.08	0.09	959.07	176.55	0.16	60	78
2	30	1092.65	469.60	0.07	1092.65	169.62	0.25	64	257
2	35	1311.25	516.32	0.09	1311.25	187.60	0.42	64	367
2	40	1401.83	545.48	0.12	1401.83	165.13	0.73	70	508
2	45	1653.96	679.77	0.13	1653.96	213.78	1.52	69	1069
2	50	1863.37	734.60	0.27	1863.37	207.94	2.00	72	641
2	55	2003.49	882.24	0.30	2003.49	223.64	1.99	75	563
2	60	2124.44	985.03	0.35	2124.44	231.49	3.61	76	931
3	20	693.33	206.18	0.04	693.33	135.89	0.13	34	225
3	25	1012.12	283.67	0.05	1012.12	194.25	0.15	32	200
3	30	1174.21	307.95	0.07	1174.21	186.16	0.33	40	371
3	35	1406.59	362.19	0.10	1406.59	204.13	0.37	44	270
3	40	1473.53	388.87	0.12	1473.53	162.08	0.47	58	292
3	45	1634.10	471.74	0.19	1634.10	203.46	2.29	57	1105
3	50	1804.43	512.30	0.22	1804.43	201.05	1.59	61	623
3	55	1946.67	601.58	0.30	1946.67	219.70	1.59	63	430
3	60	2063.12	646.33	0.34	2063.12	235.45	2.13	64	526
4	20	744.82	129.95	0.04	744.82	115.83	0.11	11	175
4	25	988.11	174.15	0.05	988.11	158.59	0.19	9	280
4	30	1150.21	198.43	0.07	1150.21	150.49	0.33	24	371
4	35	1372.56	241.42	0.10	1372.56	177.88	0.44	26	340
4	40	1436.39	260.79	0.11	1436.39	155.38	0.52	40	373
4	45	1659.52	348.32	0.19	1659.52	200.87	2.02	42	963
4	50	1818.21	368.00	0.28	1818.21	202.92	1.50	45	436
4	55	1977.30	448.16	0.23	1977.30	223.58	1.38	50	500
4	60	2098.81	477.09	0.32	2098.81	229.31	1.71	52	434
Avg.		1467.05	444.83	0.16	1467.05	188.17	1.04	50	457

, the datasets $InPm$ ($n=20,25,\cdots ,60$, and $m=2,3,4$) are tested by solving the point-to-point and cyclic bi-stage distribution modes. In the table, $z^{df}cyc$ is equal to $z^{df}ptp$, because the “$cyc$” mode only improves the ground distribution of “$ptp$”. The ground distribution distances are reduced (defined as $z^{gd*}$), and the computing time is increased (defined as $C{T}^{*}$) from “$ptp$” to “$cyc$”.

$$z^{gd*}=\dot{\frac{z^{gd}ptp-z^{gd}cyc}{z^{gd}ptp}}100/\%$$

$$C{T}^{*}=\dot{\frac{CTcyc-CTptp}{CTptp}}100/\%$$

Averagely, the ground distribution distance can be reduced by 50%, and the computing time will increase by $457\%$, as presented in the last line in Table 6.

In Figure 5, the ground distribution distances for different customers in the “$ptp$” and “$cyc$” modes are depicted. In the “$ptp$” mode, more DLPs will reduce the distances drastically, while the “$cyc$” will achieve almost equal distances.

In Figure 6, the computing times of different customers and DLPs are depicted. Solving the “$cyc$” model is time-consuming comparatively. From the points when 45 customers, it can be concluded that the solving performance is sensitive to the data. However, as a general trend, more customers involved will cost more computing time. However, computing performances can be improved by using advanced computing devices and technologies.

Comparing the results in Figure 5 and Figure 6, the ground distribution distances increase gradually with the number of customers involved, while the computing times are also increasing much more sharply and fluctuate instantly.

6.3.3. Performance Studies on the Iterative Heuristics

In Algorithm 1 (the iterative heuristics), $f^{\left[\mathrm{IH}\right]}$ solves $\left[\mathrm{M}\left(ptp,{\overline{z}}^{DLP}\right)\right]$, $\left[\mathrm{M}\left(cyc,I_d^{DLP}\right)\right]$, and $\left[\mathrm{M}\left(cyc\right)\right]$ once or several times. So, $f^{\left[\mathrm{IH}\right]}$ must increase the computation complexity, increasing computing time. $f^{\left[\mathrm{IH}\right]}$ tries to improve the result of $f^{cyc}$ and should affect the values of $z^{df},z^{gd}$, and computing time (CT). Therefore, define the following metrics ($\mathsf{\Delta }{z}^{df}, \mathsf{\Delta }{z}^{gd},\mathsf{\Delta }CT$) to evaluate their variance ratios of them.

$$\mathsf{\Delta }{z}^{df}=\frac{z^{df}\left(cyc\right)-z^{df}\left(IH\right)}{z^{df}\left(cyc\right)}\cdot 100\%$$

$$\mathsf{\Delta }{z}^{gd}=\frac{z^{gd}\left(cyc\right)-z^{gd}\left(IH\right)}{z^{gd}\left(cyc\right)}\cdot 100\%$$

$$\mathsf{\Delta }CT=\frac{CT\left(IH\right)-CT\left(cyc\right)}{z^{df}\left(cyc\right)}\cdot 100\%$$

In

Table 7. Performances of the iterative heuristics ($f^{\left[\mathrm{IH}\right]}$).

DLPs	$\mathit{I}$	${\mathit{z}}^{\mathit{d}\mathit{f}}$$\left(\mathbf{k}\mathbf{m}\right)$	$\mathbf{\Delta }{\mathit{z}}^{\mathit{d}\mathit{f}} (\%)$	${\mathit{z}}^{\mathit{g}\mathit{d}}$$\left(\mathbf{k}\mathbf{m}\right)$	$\mathbf{\Delta }{\mathit{z}}^{\mathit{g}\mathit{d}} (\%)$	CT (s)	$\mathbf{\Delta }\mathit{C}\mathit{T} (\%)$	Iterations
2	20	498.71	33.17	147.91	0.00	0.20	18.72	2
2	25	643.84	32.87	176.55	0.00	0.38	36.47	2
2	30	765.32	29.96	169.62	0.00	0.71	67.78	2
2	35	828.19	36.84	187.60	0.00	0.76	48.79	2
2	40	886.98	36.73	165.13	0.00	2.51	59.52	2
2	45	1026.08	37.96	213.78	0.00	4.38	68.62	2
2	50	1107.85	40.55	207.94	0.00	6.25	71.06	2
2	55	1231.74	38.52	223.64	0.00	6.77	68.71	2
2	60	1286.63	39.44	231.49	0.00	9.59	63.07	2
3	20	612.21	11.70	135.89	0.00	0.31	59.88	2
3	25	775.50	23.38	164.53	15.30	0.53	70.12	4
3	30	899.25	23.42	157.60	15.34	1.13	76.10	4
3	35	981.41	30.23	175.57	13.99	1.46	77.68	4
3	40	1042.47	29.25	152.39	5.98	2.75	79.35	4
3	45	1287.11	21.23	203.46	0.00	5.67	69.97	2
3	50	1402.62	22.27	201.05	0.00	4.65	67.69	2
3	55	1313.95	32.50	217.06	1.20	9.26	83.92	4
3	60	1302.57	36.86	226.44	3.83	27.86	92.54	6
4	20	661.75	11.15	128.31	1.26	0.34	69.95	3
4	25	841.54	14.83	156.95	1.03	0.44	66.70	3
4	30	967.37	15.90	150.03	0.31	0.67	56.53	3
4	35	1088.63	20.69	160.55	9.74	2.59	84.10	4
4	40	1128.47	21.44	141.04	9.23	3.12	73.64	4
4	45	1372.57	17.29	206.38	−2.74	6.65	77.90	3
4	50	1439.68	20.82	197.56	2.64	7.02	78.53	5
4	55	1500.38	24.12	210.78	5.73	8.23	86.05	5
4	60	1509.15	28.09	223.60	2.49	13.27	85.10	5
Avg.		1051.92	27.08	182.70	3.16	4.72	68.83	3.15

, the three metrics’ values of different customers and DLPs are presented, as well as the iteration times before the algorithm $\left[\mathrm{IH}\right]$ stops. Through two to five iterations (an average of 3.15 times), the drone flying distances can always be reduced by an average of 27.80%. The ground distribution distances can be reduced by 3.16% on average. Notably, some datasets incur no improvement, when ${\overline{z}}^{DLP}=2$ or both customers and DLPs are not many. In addition, the ground distribution distance may be increased when saving the drone flying distance, as the result of dataset $I45P4$. The computing time ($CT$) will be increased by an average of 68.83%.

6.4. Discussions and Managerial Implications

Drones are applied to various application scenarios, including transportation and logistics. Generally, drones are capable of last-mile or final-stage logistics because their capacities are limited [45]. However, in many logistics scenarios, the packages are limited in size and quantities, and drones can be used in various stages of logistics systems. Island distribution is a special application scenario of drones because the sea separates the land and the island. Passenger and cargo transportation primarily depends on ships, while waterborne transportation is sensitive to climate conditions and cannot be suitable for urgent demands. In COVID-19, drones can be used to distribute viral tests to potentially infected patients [13]. Under these considerations, drone-based island distribution should be an effective solution.

We make the following generations as managerial implications based on the experimental results.

(1) As studied in Section 3, when considering drone-based distribution and logistics solutions, we can meet various solutions. Drones can help solve certain critical problems, but not all. We can integrate drones with other means of logistics to develop a holistic solution. Furthermore, the present facility and variable costs should all be considered and balanced.

(2) As studied in Section 4 and Section 5, the direct distribution mode uses drones only to provide terminal-to-customer package delivery. It is simple but involves great investment in DLPs that are expensive and fixed facilities, not so suitable for application scenarios when the technologies are still fast developing. The bi-stage distribution mode includes basic and improved modes, and their models can be solved by three methods. They use drones, DLPs, and general ground distribution to combine the different advantages and contribute to adjustable costs.

(3) As studied in Section 6, the proposed models and methods can achieve promising performances for small and medium-scale test instances. They can be used on small islands. Different drone-based distribution modes present different evaluation values in the proposed metrics (as studied in Section 3.3). The cyclic bi-stage distribution mode can decrease ground logistics costs. It can be further improved by iterative heuristics.

The regulatory agencies, numerous States, and entities are involved in the creation of safe integration with manned aviation [5]. Recent conceptual and regulatory advancements in the field of Urban Air Mobility (UAM) in Europe were elaborated to outline the digital ecosystem in which aviation and non-aviation actors would exchange information to ensure operations’ efficiency, safety, and regulatory compliance [46].

7. Conclusions

We developed the study on a drone-based logistics system for island emergency distribution, considering the specialties of islands. First, the study was conducted after examining the relations studies, technologies, and applications in the area of drone-based logistics. We investigated many experimental systems on drone-based logistics, which were created by the leading companies in package delivery companies. Moreover, drone bases are testing transporting cargo from land to islands. During the COVID-19 epidemic days, some islands face shortages of resources and medicals, which encouraged the study. Second, three drone-based distribution modes are proposed to categorize the various possibilities of using drones in islands: the direct distribution mode, the point-to-point, and the cyclic bi-stage distribution modes. The bi-stage modes use drones for the first-stage land-to-island transportation, and use ground transportation to serve the second-stage distribution, while the DLPs mediate the two stages. We also developed four metrics to evaluate the modes and instructed the modeling, algorithm devices, and comparison studies. Third, solution methods were developed to solve the models, especially for the bi-stage mods. An iterative heuristic was designed to improve the cyclic bi-stage distribution mode. Finally, using the famous “Mount Putuo” island as a case, we developed dataset generation methods and conducted a series of experiments to demonstrate and investigate the proposed modes, models, and solution methods.

This study can be extended in the following aspects in future studies. First, in Section 4.2, the point-to-point bi-stage distribution mode is formulated as a multi-objective program $\left[\mathrm{M}\left(ptp\right)\right]$. We further formatted a single-objective linear program $\left[\mathrm{M}\left(ptp,{\overline{z}}^{DLP}\right)\right]$ based on it, where ${\overline{z}}^{DLP}$ sets the DLPs. It is beneficial to develop an algorithm for solving the multi-objective program and investigate the tradeoffs among the three objectives. Second, in Section 4.3, the cyclic bi-state distribution mode is formulated as separated TSP models for the determined DLPs and assignments of customers to the DLPs. In future studies, we can formulate multi-depot TSPs considering drone traveling distances to the depots. Third, the iterative heuristics developed in Section 5 is effective, while their optimality cannot be ensured. Furthermore, the terminal criteria of the heuristic algorithm consider the optimization of the two costs, drone traveling and ground distribution distances. Indeed, these two costs can be united, or we can test their tradeoffs. Moreover, new algorithms can be devised to find optimal solutions by using global optimization, other than multi-stage solution methods.