A bi-objective optimization model for the medical supplies' simultaneous pickup and delivery with drones.

In the COVID-19 pandemic, it is essential to transport medical supplies to specific locations accurately, safely, and promptly on time. The application of drones for medical supplies delivery can break ground traffic restrictions, shorten delivery time, and achieve the goal of contactless delivery to reduce the likelihood of contacting COVID-19 patients. However, the existing optimization model for drone delivery is cannot meet the requirements of medical supplies delivery in public health emergencies. Therefore, this paper proposes a bi-objective mixed integer programming model for the multi-trip drone location routing problem, which allows simultaneous pick-up and delivery, and shorten the time to deliver medical supplies in the right place. Then, a modified NSGA-II (Non-dominated Sorting Genetic Algorithm II) which includes double-layer coding, is designed to solve the model. This paper also conducts multiple sets of data experiments to verify the performance of modified NSGA-II. Comparing with separate pickup and delivery modes, this study demonstrates that the proposed optimization model with simultaneous pickup and delivery mode achieves a shorter time, is safer, and saves more resources. Finally, the sensitivity analysis is conducted by changing some parameters, and providing some reference suggestions for medical supplies delivery management via drones.

Introduction

The global outbreak of COVID-19 causes serious problems such as limited material distribution and increased transportation time. All these results in unnecessary deaths of infected persons (Zhu et al., 2020). As social contact is the leading cause of COVID-19, medical doctors suggest that maintaining social distance is the best way to limit its spread (Mitrokhin et al., 2020). Thus, transporting medical supplies quickly and accurately to specific locations in the isolation area is crucial in health care services. On the other hand, it is challenging to take test samples such as COVID-19 to the testing laboratory for a quick inspection to timely identify and isolate the infected patients, and provide medical services as soon as possible (Koshta et al., 2021). In short, it is imperative to improve the efficiency of the pickup and delivery of medical supplies in times of pandemics such as COVID-19.

In the medical field, drones (or unmanned aerial vehicles, UAVs) play a vital role in delivering medical supplies safely and quickly (Khan et al., 2021). Drones can help medical staff perform tasks more efficiently and save lives (Nyaaba and Ayamga, 2021). The benefits of using drones to transport medical supplies include (Koshta et al., 2021, Moshref-Javadi and Winkenbach, 2021): (1) Break the restrictions of ground transportation and effectively minimize the impact of road restrictions and community closures during the lockdown. (2) Realize contactless delivery, avoid face-to-face contact between delivery personnel and medical staff, and reduce the risk of cross-infection. (3) Reduce delivery time and deliver emergency supplies to designated locations in the shortest time.

A drone is a promising technology to combat life-threatening public health emergencies such as COVID-19, and it is increasingly used to deliver medical supplies (Glick et al., 2021). A representative application case of using drones to transport medical supplies was Zipline, which used its fixed-wing drones to transport blood, medicines and vaccines to remote areas (Magdalena and Lora, 2021). Ghana applied drones to transport test samples to the testing laboratory, as shown in Fig. 1 a (Aryn, 2020). UPS cooperated with Matternet to transport medical samples with drones in the WakeMed Park in North Carolina (Hazel, 2019). Since the outbreak of COVID-19, SF Technology has built a non-contact aerial transport bridge in the epidemic area in China, as shown in Fig. 1 b (SF-Technology, 2021). They invested drones in five cities, including Wuhan, to transport medical supplies, such as protective clothing, gloves, food, and medicines. The total amount of medical supplies transported in a single day was 1.8 tons, and the total flight mileage exceeded 1,500 km.

Fig. 1

Examples of using drones to transport medical supplies.

From the perspective of actual demands and application prospects, the application of drones in medical supplies delivery is necessary and feasible. However, many problems remain unresolved in dispatch operation, such as pickup and delivery route planning. In this paper, we propose a new route planning model and algorithm to optimize the use of multiple drones to simultaneously deliver medical supplies (such as protective clothing, gloves, vaccines, blood, etc.) from the central hospitals to the epidemic prevention nodes and pickup test samples from the epidemic prevention nodes to the central hospitals.

The core of the decision is to assign a sequence of pickup and delivery tasks to each drone while meeting the dual time requirements of medical supplies delivery and test samples recovery. Considering the characteristics of drone, such as limited endurance time and small payload capacity, the decisions involved must include the followings: which central hospitals should be equipped with drones, how many drones should be deployed in each central hospital, and how to plan the trips between each central hospital and the epidemic prevention node, whether drones are reused, and how to take into account operating costs and the efficiency of pickup and delivery in the delivery plan.

The rest of this paper is organized as follows. Section 2 reviews previous research related to this study. Section 3 describes our problem and establishes a mixed integer programming (MIP) model for drone delivery. Section 4 introduces the solution algorithm based on NSGA-II. Section 5 presents numerical experiments and analyses. Next is Section 6, which summarizes this paper.

Literature review

This section summarizes current work relevant to the research problem in this paper, including the application of drones in medical supplies delivery, vehicle routing problem with drone (VRPD), and vehicle routing with simultaneous pickup and delivery (VRPSPD).

The application of drones in the medical supplies delivery

In recent years, distributing medical supplies safely and efficiently against major public health emergencies has become a hot issue. It was found that timely relief supplies are the key to reducing loss of life and other impacts in response to public health emergencies. There were severe problems, such as limited distribution of supplies between hospitals and the untimely transportation of medical supplies, which have reduced the safety of medical supplies and the cure rate of patients. To address this problem, Liu et al. (2021a) proposed a method for dispatching medical supplies in major public health emergencies to improve the efficiency of medical supplies distribution. Wang et al. (2021a) constructed a Markov Decision Process model to study emergency medical supplies scheduling strategies for COVID-19 and other public health emergencies. The model showed how to dispatch limited emergency medical supplies to optimize the service rate of the entire system. Patil et al. (2021) also pointed out that medical supplies are critical in humanitarian relief efforts. They analyzed 20 potential obstacles to the sustainable development of the medical supply chain, and the results showed that material, operations, and logistical are the main influencing factors. Martins et al. (2021) introduced the concept of agile optimization for distributing post-disaster medical supplies because every second can decisively save lives in humanitarian logistics.

Therefore, it is urgent to adopt new technologies and ideas to improve the efficiency of medical supplies distribution to control the safety problems and losses caused by the untimely delivery of medical supplies. There is increasing evidence that drones play a crucial role in providing medical health, and drones can save more lives by reducing emergency response time (Euchi, 2021, Nyaaba and Ayamga, 2021, Otto et al., 2018). In order to combat the impact of COVID-19, Koshta et al. (2021) believed that drone is a viable option to improve the efficiency and effectiveness of the humanitarian aid supply chain, and one of the main application areas is the medical health supply chain. Ayamga et al. (2021) also stated that drones could provide logistics solutions for personal protective equipment, test kits, vaccines, drugs, and laboratory samples in the COVID-19 era. However, expected negative consequences of urban innovation involving advanced technologies include failure of algorithmic decision making (Yigitcanlar et al., 2021) may lead to drone malfunctions and delays, which can be fatal for medical supplies. Other problems include falling objects from height, which can cause the deaths of humans and animals on the ground. Therefore, Glick et al. (2021) developed a modelling framework to measure the delivery reliability of drones with random demand and meteorological conditions. They analyzed the trade-offs among drone reliability, fleet size, population size, and meteorological conditions. Yakushiji et al., 2020a, Yakushiji et al., 2020b carried out a series of drone transportation tests to prove using drones to transport emergency supplies. Their experimental results indicated that drones' transportation of medical supplies positively impacts the medical health system. The advantages of drones in improving efficiency and reducing costs in various applications were introduced by Ling and Draghic (2019), including blood delivery, laboratory testing, medical equipment delivery, and drug delivery. Haidari et al. (2016) applied a simulation model to evaluate the impact of using drones for vaccine distribution. They concluded that drone systems could save up to $0.21 in logistics costs per dose of vaccine and increase vaccine availability by 2% compared with the traditional multi-tiered land transport system.

Technologies like drones develop rapidly and bring many opportunities (Li, 2018, Luo et al., 2022). Nevertheless, immature robotics technology causes unstable movement performance (Yang et al., 2022), and the use of drones as a safer means of delivery than traditional ones is a research focus for future applications (Otto et al., 2018). In addition, the COVID-19 pandemic still impacts human life and safety on a global scale, and the need for social distance has created space for new applications of drones (Koshta et al., 2021). Safety, timeliness and essential economic requirements of medical supplies delivery force the application of drones to be of utmost importance. Distributing medical supplies quickly and accurately through new technologies in significant public health emergencies is an important field that will continue to be studied. Based on this, this paper discusses the application of drones in constructing urban aerial transport bridges for medical supplies delivery.

Vehicle routing problem with drone

In the existing research, the logistics model of VRPD can be divided into pure-play drone-based models, unsynchronised multi-modal models, synchronised multi-modal models, and resupply multi-modal models (Moshref-Javadi and Winkenbach, 2021). Only the first model type applies to the drone-only delivery problem (DDP). The drones are used to deliver packages directly from the warehouse to the customer. The remaining three models are all aimed at the truck-drone problems, which means drone delivery needs to cooperate with trucks (Jeong et al., 2019, Kuo et al., 2022, Murray and Chu, 2015). This section mainly reviews the research on DDP related to this paper.

Studies on DDP typically assume that there are multiple drones in the warehouse, and each drone can serve one or more nodes at a time. Dorling et al. (2017) proposed two variants of drone delivery's vehicle routing problem model. The first is to minimize total operating costs with delivery time constraints, and the second is to optimize delivery time with budget constraints. They introduced a linear approximation function to calculate the power consumption that varies linearly with the payload and battery weight. Furthermore, a simulated annealing heuristic algorithm was designed to solve the model. Cheng et al. (2020) further expanded the research of Dorling et al. (2017). They expressed the power consumption of the drone as a nonlinear function of payload and travel time in the multi-trip drone routing problem model considering the time window. The logical cut and subgradient cut were introduced to process the nonlinear power function, and the branch-and-cut algorithm was used to solve the drone routing problem. An automated drone delivery system was studied by Choi and Schonfeld (2017), in which all customers' demands were the same. They utilized the relationship between battery capacity, payload and flight range to optimize the drone fleet size. Song et al. (2018) constructed a mathematical model aiming at the maximum weighted sum of the total number of covered tasks and total travel distance to support drone delivery logistics. Considering that drones are limited by loading capacity and flight time, the model allowed drones to share multiple sites. Ham (2020) defined a drone scheduling problem for a warehouse material transfer system to improve efficiency and reduce costs. Yakici (2016) considered an integer linear programming model to optimize drone delivery based on location routing problems. The optimization goal of the model was to maximize the sum of the importance values ​​corresponding to the covered nodes, and an ant colony optimization metaheuristic algorithm was designed to solve the problem. However, the nodes with low importance values ​​may not be served in their research, which is a situation that cannot occur in the delivery of medical supplies.

Operations research techniques and optimization modelling have been heavily used in scenarios where drones are used in healthcare and humanitarian logistics (Macrina et al., 2020, Rejeb et al., 2021). Rabta et al. (2018) studied the application of drones in the last mile delivery of humanitarian logistics in a single-depot multi-node network. The goal of their optimization model was to minimize the total travel distance of drones subject to the constraints of payload and energy. In the application of maritime search and rescue, Cho et al. (2021) analyzed the coverage path planning problem for a multi-UAV area. They introduced a two-stage approach, where the first stage was used for area decomposition and the second stage achieved the optimal coverage path to minimize time. Kim et al. (2017) described drone-aided healthcare services implemented in rural areas and introduced two models for finding the optimal number of drone centers and route planning. Focused on the blood supply problem in an emergency, Wen et al. (2016) proposed a UAV-based capacitated vehicle routing problem model with distance cost and flight times as optimization objectives. Chowdhury et al., 2021, Chowdhury et al., 2017 considered many drone trajectory-specific factors, such as battery recharging, drone hovering and many others in detail. Chowdhury et al. (2021) proposed a mixed-integer linear programming model with the objective of the post-disaster inspection cost due to these factors for the heterogeneous fixed fleet drone routing problem. Chowdhury et al. (2017) used a continuous approximation model to determine the optimal location of a distribution center, corresponding emergency inventory and service area. Gentili et al., 2022, Kim et al., 2019, Shavarani et al., 2021 studied the capacity location problem with a drone. (Shavarani et al., 2021) considered two optimization objectives: minimizing the total cost of facility construction and drone procurement, and the number of uncovered customers. Kim et al. (2019) probed the uncertainty of flight distance caused by battery consumption fluctuations in the modelling framework. The optimization problem in (Gentili et al., 2022) aimed to minimize the total disutility value based on the perishability of emergency medical supplies. They assumed that each platform only had one drone and could serve one node at a time. Considering the battery life of drones during medical supplies transport, Dhote and Limbourg, 2020, Ghelichi et al., 2021, Macias et al., 2020 jointly optimized the location selection of charging stations or tactical centers and the route of drones.

In the literature mentioned above, we detect that battery energy capacity, payload, and flight range are the main constraints considered in drone delivery research, which is also one of the reasons why the traditional vehicle routing problem model is not suitable for drone logistics. Scholars have tried to expand the research of DDPs by considering features such as multi-trip of drones, integrated optimization of location and routing problems, and multiple optimization objectives. However, models that include multiple features simultaneously are still rare. Optimization modelling is an accepted and general approach regarding drone utilization in medical supply or healthcare service delivery. Most of the existing studies analyzed the problem from the perspective of drone specificity, and the coordination between the timeliness of delivery and drone scheduling is ignored.

Vehicle routing with simultaneous pickup and delivery

VRPSPD is essential for the pickup and delivery problem (PDP) (Gutiérrez-Sánchez and Rocha-Medina, 2022). Research on VRPSPD originated in (Min, 1989) and has developed into a fruitful and active research area. The existing variants and extensions of VRPSPD include VRPSPD with time windows, heterogeneous VRPSPD, multi-depot VRPSPD, green VRPSPD, stochastic VRPSPD and others (Koç et al., 2020).

VRPSPD is widely used to analyze problems in logistics systems and various industries because of its practical importance and benefits. Du et al., 2021, Hornstra et al., 2020 proposed a mathematical planning model to solve the last mile delivery problem by combining the use of parcel lockers with VRPSPD. Hornstra et al. (2020) considered three types of customers according to the parcel pickup and delivery methods. Du et al. (2021) analyzed the impact of real-time capacity changes on the courier service flow. Parcel lockers, alternative delivery points and time windows were simultaneously considered in the studied VRPSPD of (Sitek et al., 2021). Liu et al. (2021b) introduced the VRPSPD that considered real-time traffic conditions, aiming to determine the optimal vehicle route planning with minimum total travel period. They transformed the dynamic problem into a sequence of re-optimization VRPSDP problems by splitting time periods. Park et al. (2021) proposed a waiting strategy based on a rerouting indicator for solving VRPSPD under dynamic demand constraints. Chaieb and Ben Sassi (2021) discussed the application of VRPPD with time windows in home healthcare, involving logistical actions including the delivery of medical services to patients and the transport of patients' unused medications and biological samples to hospitals. Fazi et al. (2020) explored a VRPSPD for optimizing inland container shipping to find the best container allocation, thereby improving barge utilization. Calculating the total fuel consumption cost or the carbon emission cost caused by fuel consumption is a common approach in green VRPSPD (Olgun et al., 2021, Qin et al., 2019). Some interesting extensions on VRPSPD are also two-echelon VRPSPD (Liu and Jiang, 2022, Luo et al., 2021), VRPSPD with inventory and location (Sherif et al., 2021, Zarrat Dakhely Parast et al., 2021), and VRPSPD with drone (Ham, 2018, Martins et al., 2021). Ham (2018) used a constraint programming method to extend the parallel drone scheduling traveling salesman problem to consider drop and pickup synchronization, which is applied to solving multi-truck, multi-drone, and multi-depot scheduling problem. Martins et al. (2021) studied a two-echelon VRPSPD using drones. However, their pickup and delivery activities only occurred at intermediate nodes, and the demands of all customer nodes have not been followed.

VRPSPD is an NP-hard optimization problem as it is developed from vehicle routing problem (VRP) (Luo et al., 2021). Since metaheuristics have the advantage of searching for near-optimal solutions with acceptable quality in a reasonable time (Cordeau et al., 2002), many researchers have applied metaheuristics to solve different VRPs, especially problems with multiple objectives. Yu et al., 2022, Sherif et al., 2021 applied the simulated annealing algorithm to VRPSPD. Chaieb and Ben Sassi, 2021, Shi et al., 2020 improved the tabu search metaheuristic. Liu et al. (2021b) proposed a hybrid algorithm combining the ant colony system and the virtual transformation method. Hornstra et al., 2020, Hof and Schneider, 2019 explored the optimal solution of VRPSPD by improving the adaptive large neighborhood search metaheuristic.

Among numerous known metaheuristics, the genetic algorithm (GA) has been successfully applied to solve single-objective and multi-objective VRPs (Koç et al., 2020). In the work of (Park et al., 2021, Sitek et al., 2021), GA was applied to solve the VRPSPD with minimized operating costs. Luo et al., 2021, Wang et al., 2021b integrated a clustering algorithm and an improved NSGA-II to analyze a two-echelon VRPSPD with multiple objectives. Shavarani et al. (2021) selected NSGA-II and NSGA-III for the bi-objective location problem of drone delivery. Maskooki et al. (2022) designed a customized genetic algorithm with NSGA-II as the base multi-objective evolutionary algorithm for the dynamic bi-objective routing problem. Available experimental results show that NSGA-II can solve complex models for real-world problems. Therefore, combining the features of the proposed bi-objective VRPSPD with drone, we design a modified NSGA-II based on the general framework of NSGA-II. Two classical multi-objective evolutionary algorithms, the multi-objective evolutionary algorithm based on dominance and decomposition (MOEA/DD) (Ishibuchi et al., 2017, Wen et al., 2016) and the improved decomposition-based evolutionary algorithm (I-DBEA) (Ali et al., 2021, Anwar and Younas, 2020), were used for performance comparison.

Available VRPSPDs are inadequate for planning DDP for medical supplies delivery: either they do not allow for drone reuse, resulting in solutions that use too many drones, or they do not consider the impact of battery and payload weight on energy consumption, resulting in costly or infeasible routes; either the time window is ignored, affecting the efficiency of material delivery, or the use of drones depends on trucks. In addition, minimizing cost as the optimization objective is still the popular modelling approach. Multi-objective problems have been considered in some studies. However, there are no studies integrate multiple real-world complexities, including multi-depot, multi-drone, multi-trip, dual time windows, simultaneous pickup and delivery, and joint optimization of location and routes.

From what has been reviewed and discussed above, we are convinced of the critical role of new technologies and ideas in strengthening the safety and timeliness of medical material delivery. There have been research attempts to discuss and analyze drones to transport medical supplies in recent years. However, the application of drones in medical supplies delivery is still in its infancy stage. The existing optimization model is not sufficient to match the demand for delivery of medical supplies in public health emergencies. Optimizing the scheduling plan of drones is the key to improving the pickup and delivery efficiency and reducing operating costs. This paper will comprehensively consider the characteristics of the drones' multi-trip, simultaneous pickup and delivery requirements of nodes, and time requirements for medical supplies delivery to construct a bi-objective location and routing joint optimization model. The research provides a more general solution for medical supplies delivery with drone, where drones can be reused, rather than simply assuming that there is no limit on the number of drones in each depot. Each drone operation can cover the simultaneous pickup and delivery demands of one or more nodes, which takes the payload and energy consumption of the drone as the constraint premise.

The main scientific contribution of this paper can be summarized as follows.

First, this paper provides a new solution for drones’ route planning, which indicates how drones can be used as an aerial transport bridge for the medical supplies delivery in public health emergencies, thereby improving the efficiency of medical supplies delivery, ensuring the timeliness of sample testing, and reducing potential social contact and the risk of infection.

Second, the DDP studied in this paper can be defined explicitly as the multi-trip drone location routing problem with simultaneous pickup and delivery (MT-DLRP-SPD), an extension of VRPD and VRPSPD. It is of academic significance owing to the specific characteristics of VRPD and VRPSPD summarized from the literature.

Third, a modified NSGA-II (M−NSGA−II) is introduced. Numerical experiments are performed with multiple data sets to solve and verify the proposed model. Moreover, the M−NSGA−II is compared with the multi-objective solution algorithms MOEA/DD and I-DBEA. Experimental results prove that the new algorithm has better performance.

Optimization model

Model description

The investigated problem in this paper is how to dispatch drones in multiple central hospitals (CHs) and plan the routes of multiple drones to cover the demands of epidemic prevention nodes (EPNs) for medical supplies delivery and samples pickup in a safe, efficient and economical manner. Here, we assume that all medical supplies are sent from the CHs, and the pickup samples are sent back to the CHs for inspection and quarantine. Drone delivery has the advantage of speed so that they can be reused in a scheduling task. The CH also works as a charging station. Each drone may return to the CH many times to load medical supplies, unload samples and replenish energy. Fig. 2 illustrates the main process of medical supplies delivery by drones.

Fig. 2

Schematic diagram of drone delivery.

Therefore, the proposed MT-DLRP-SPD is defined on a directed graph . Where represents the set of network nodes. is the set of candidate CHs, and is the set of EPNs. Each node in set has two attributes: the quantity to be delivered (kg) and the quantity to be picked up (kg). and respectively correspond to a latest delivery time (s) and (s), where must be sent to the EPN before , and must be returned to a CH before . The arc set is defined as , each arc is associated with a travel distance (m), a travel time (s) and a travel cost . represents the set of available drones in the CH. Each drone can perform multiple trips and visit multiple EPNs in one trip. It is assumed that the drones deployed in CHs are homogeneous. When the drone departs from the CH, it is fully charged. We assume the drones are travelling between nodes at a constant speed, and the driving state is not affected by the external environment such as weather.

The problem consists in choosing which CHs to participate in the service or can be defined as which CHs to deploy drones, deciding each CH services for which EPNs, and designing the drone dispatch plan and routes such that the objective functions are optimized, and the following constraints are satisfied: (1) Each route starts from a CH and ends at the CH where it started. (2) Each EPN is only visited once to cover its pickup and delivery demands. (3) The drone payload constraint, battery energy constraint, and the latest delivery time requirements of medical supplies and samples must be respected.

Model formulation

In order to improve the economy and timeliness of medical supplies delivery by drones on the premise of ensuring safety, we construct a bi-objective MIP model for the proposed MT-DLRP-SPD (BOMT-DLRP-SPD). The main objective of this MIP model is to determine the optimal location routing plan with minimum total operating cost and total drone travel time for the MT-DLRP-SPD. The main sets, parameters, and decision variables used to build this model are given in Table 1 .

Table 1

Notations in the model.

Notations	Definition
Set
H	Set of candidate CHs,H=1,2,...,m.
P	Set of EPNs,P=1,2,...,n.
L	Set of network nodes,L=H∪P.
K_m	Set of available drones in the CHm,K_m=1,2,...,k,m∈H.
Parameter
K	Number of available drones in the CHs.
Q	Largest payload of drones.
q_i	Quantity to be delivered at EPNi,q_i⩽Q.
p_i	Quantity to be picked up at EPNi,p_i⩽Q.
t_q_i	Latest delivery time of q_i.
t_p_i	Latest delivery time of p_i.
θ_i	Service time of EPNi.
d_ij	Travel distance on arc i,j.
v	Cruise speed of the drones.
t_ij	Travel time on arc i,j,t_ij=d_ij/v.
E_max	Maximum battery energy capacity of drones.
c_ij	Unit travel cost.
c_mf	Fixed cost of opening a CH.
c_df	Fixed cost of using a drone.
c_e	Unit energy cost.
f_i	Cumulative energy consumption when the drone arrives at the EPNi. In particular, f_im represents the cumulative energy consumption when the drone returns to the CH m from the EPNi, and f_mi represents the cumulative energy consumption before the drone departs from the CHm.
w_b	Battery weight of the drones.
qijmk	Load of the drone traveling on the arc i,j.
α	Energy consumption parameter.
β	Energy consumption parameter.
t_i	The time when the drone arrives at the EPNi. In particular, T_im represents the time when the drone returns to the CH m from the EPNi.
Daimk	Remaining delivery quantity of the drone when it arrives at the EPNi.
Dlimk	Remaining delivery quantity of the drone when it leaves from the EPNi.
Paimk	Picked up quantity of the drone when it arrives at the EPNi.
Plimk	Picked up quantity of the drone when it leaves from the EPNi.
B_m	The time to replenish energy, load and unload goods in the CHm.
M	An infinite number.
Decision variable
o_m	o_m=1 if CH m participates in the service, 0 otherwise.
u_mk	u_mk=1 if drone k of CH m is used, 0 otherwise.
xijmk	xijmk=1 if arc i,j is traversed by drone k of CHm, 0 otherwise.
rijmk	rijmk=1 if drone k of CH m returns to the CH from EPN i to load medical supplies, unload samples and replenish energy, and then starts a new route with EPN j as the first customer, 0 otherwise.

Objective function 1. The operating cost function consists of two parts: Fixed costs , including the fixed costs of opening CHs and the fixed costs of using drones. Variable costs C, including the driving costs and energy costs of drones.

(1)

(2)

(3)where the fixed costs of opening a CH are the costs associated with the exclusive sites and staff involved in deploying drones in the CH. The fixed costs of drones are related to the number of drones used and refer to the drone depreciation (including administrative procedures and capital cost), insurance premium, maintenance costs (Dhote and Limbourg, 2020). Since we discuss multiple trips of the drone, the number of drones used is calculated by subtracting the number of drones reused from the number of drones departed, which is expressed as in Eq. (1).

Based on the analysis of Dorling et al. (2017), the energy consumption of the n-rotor drone is related to the battery weight and the load on the arc , which can be approximately calculated as:

(4)

Objective function 2. The travel time can be calculated as the sum of the working hours of all drones:

(5)

The working time of a drone is calculated from the time of departure from the CH until the drone returns to CH and is no longer in use. If the drone is used only once, its working time is . In the multi-trip problem, the working time of a drone is calculated continuously based on its first return time to CH , denoted by . The working time of a drone in Eq. (5) is composed of two parts: the working time of the drone for a single use and multiple trips. Therefore, if the drone is reused, the time of its first return to the CH should be subtracted to avoid double counting, expressed as in Eq. (5).

Constraints. We organize the constraints into five groups:

1) Route feasibility constraints

(6)

(7)

(8)

(9)

(10)

(11)

Constraint (6) forbids the unopened CHs send out drones. Constraint (7) forces that CHs cannot serve each other and there is no connection path between CHs. Constraint (8) ensures that each node can only be accessed once. Constraint (9) indicates that if a drone serves a node, it should leave from the same node to ensure the continuity of the route. Constraint (10) specifies that the number of trips for a drone leaving a CH equals to the number of trips back to the CH. Constraint (11) imposes that a node is accessed by the drone which departs from the CH.

(2) Load constraints

(12)

(13)

(14)

(15)

(16)

(17)

Constraints (12) (13) state the load situation when the drone departs from the CH. Constraint (14) confirms the relationship between the remaining delivery quantity on the arcs before and after the drone passes each node. Constraint (15) denotes the relationship between the picked up quantity on the arcs before and after the drone passes each node. Constraint (16) represents the sum of the remaining and picked up quantity on each arc — the total loading of the drone cannot exceed the maximum payload of the drone. Constraint (17) indicates that if the drone passes through the arc , then the net quantity of node h is less than the remaining capacity of the drone.

(3) Reusability constraints

(18)

(19)

(20)

Constraint (18) indicates that if the drone returns to the CH from the node, it can be used again to drive to another location. Constraint (19) ensures that if the reused drone leaves the CH, it has previously arrived from another location. Constraint (20) limits the number of drones used.

(4) Time constraints

(21)

(22)

(23)

(24)

Constraint (21) denotes the time relationship between nodes. If the drone traverses arc , the time to reach j equals to the sum of the time to reach i, the service time at node i, and the travel time of arc . Constraint (22) establishes the time relationship when the drone continues on another trip after returning to the CH. If the drone is reused and travels to node j, the time to reach node j equals to the sum of the time when the drone first returns to CH, the time to replenish energy and load and unload cargo at CH, and the travel time of the arc . Constraints (23) (24) limit the latest delivery time of goods delivered and picked up.

(5) Energy constraints

(25)

(26)

(27)

(28)

(29)

(30)

Constraint (25) indicates that the cumulative energy consumption at the beginning of each trip is 0. Constraint (26) ensures the energy relationship between nodes. If the drone traverses an arc , the cumulative energy consumption to reach j equals to the sum of the cumulative energy consumption to reach i and the energy consumption while travelling on the arc .Constraint (27) requires that the constraint for maximum battery energy capacity must be observed. Constraints (28) - (30) specify the variable domain.

In summary, the BOMT-DLRP-SPD model can be expressed as:

(31)

(32)

Algorithm design

NSGA-II is a metaheuristic algorithm based on non-dominated sorting and is regarded as an essential milestone in multi-objective evolutionary optimization (Deb et al., 2002).

A fast non-dominated sorting procedure is proposed. On the one hand, it reduces computational complexity. On the other hand, it merges the parent population with the child population, so that the next generation population is selected from the double space, thus retaining all the best individuals.

A crowded-comparison approach is designed. The crowded distance is based on the average distance between points. No user-defined parameters are required to maintain the diversity among population members.

The crowded-comparison operator, which combines the results of non-dominated sorting and crowding distance, is introduced. It is used as the comparison standard between individuals in the population so that the individuals in the quasi-Pareto domain can be uniformly extended to the entire Pareto domain, ensuring the population's diversity.

NSGA-II has recognized advantages in solving bi-objective optimization problems (Bandyopadhyay and Bhattacharya, 2013). Hence, we introduce a modified NSGA-II (M−NSGA−II) for the proposed BOMT-DLRP-SPD model. First, different from the polynomial mutation strategy used by NSGA-II, M−NSGA−II adopts a Gaussian mutation operator (Sun and Gao, 2019) to improve the local search performance for the focal search region. Second, the NSGA-II algorithm with a fixed mutation rate does not always find the final solution to the optimization problem. The adaptive strategy of M−NSGA−II addresses this problem by adjusting the mutation rate (Yi et al., 2018). Third, M−NSGA−II introduces a diversity maintenance strategy based on dynamic congestion distance (Jeyadevi et al., 2011) to obtain a Pareto front with high uniformity. The calculation process of M−NSGA−II is shown in Fig. 3 .

Fig. 3

The flowchart of the M−NSGA−II.

Step1. Chromosome coding and decoding. The real-coded method is adopted, and the coding includes two layers. The first layer (F1) is the trip and route plan, which aims to allocate CHs and EPNs to different trips and determine the service order of EPNs on each trip. The second layer (F2) is the drone trip allocation, which aims to match each trip with a drone.

F1: trip and route planning. If we assume that the number of EPNs is and the maximum number of trips is , then is the coding dimension of the problem, and the value range of each dimension is . The random-key decoding method is used, that is, the smallest value in the code corresponds to 1, and the second smallest value corresponds to 2. By analogy, the code can be mapped to an arrangement on , representing EPN and the trip segmentation symbol. The delivery routes of trips can be obtained by using the segmentation symbol. Since the model requires the drones to return to its starting CH, the distance from each CH to the first EPN in each trip is calculated after obtaining trips to select the starting CH according to the principle of proximity.

For example, Fig. 4 displays the coding and decoding process for ,, where the coding dimension is . With known code, arrangement (4,1,13,2,6,8,9,11,3,12,10,5,7) can be obtained by decoding, where 1 ∼ 10 are EPN, and 11 ∼ 13 are trip segmentation symbol. There are four trips and their routes. By calculating the distances from 4, 2, 3, and 10 to each candidate CH, the CH to which the four trips belong can be determined.

Fig. 4

Examples of coding and decoding.

F2: drone trip allocation. The coding dimension is , and the value range of each dimension is . If the code value is greater than 0.5, the trips are merged, otherwise the trips are not merged. For instance, the coding dimension for is 3x1. Given the code (0.15, 0.08, 0.62), trips 1 and 2 are not merged, trips 2 and 3 are not merged, and trips 3 and 4 are merged. Trips 3 and 4 are performed by the same drone in the same CH.

Step 2. Generate initial parent population . We combine the ordered generation strategy and the random generation strategy to generate the initial population. First, the greedy strategy is used to generate 30% of the initial individuals in an orderly manner to ensure the effectiveness and quality of the initial population. Second, the random generation strategy is used to generate the remaining 70% of the initial individuals to maintain the diversity of the initial population.

Step 3. Generate offspring population through a genetic operation. The genetic operation is performed by the binary tournament selection operator (Deb et al., 2002), the simulated binary crossover operator (Zhao et al., 2019) and the Gaussian mutation operator (Sun and Gao, 2019). The selection criteria are executed based on the crowding distance. If individuals belong to different non-dominated ranks, the individual with a better rank (lower value) is preferred. If the individuals belong to the same non-dominated rank, the individual with a smaller crowding distance is selected first. Moreover, an adaptive mutation strategy (Yi et al., 2018) is introduced into the mutation operator, and the update rule of mutation probability p is shown in Eq. (33).

(33)where and represents the current generation and the maximum generation respectively. is a fixed real number, is the dimension of the problem, O = 2 in this paper.

Step 4. Non-dominated sorting. The combined population is sorted based on non-dominated conditions. The size of is 2 N. The individuals in are divided into several different non-dominated ranks .

Step 5. Individuals are selected to generate a new population according to non-dominated sorting and crowding distance results. In order to better maintain the horizontal diversity of Pareto-front, Jeyadevi et al., 2011, Luo et al., 2008 introduced dynamic crowding distance (DCD) based on the research of Deb et al. (2002). When selecting individuals to enter the new population, the individual with the lowest DCD value is deleted each time, and the DCD of the remaining individuals are recalculated until .The DCD of the individuals in the population can be calculated as:

(34)

(35)

(36)where is the crowding distance of individual . is the objective of the individual after sorting the population according to the crowding distance, and is the objective of the individual . is the variance of the crowding distance for adjacent individuals of individual , which provides information about the difference variations of crowding distance in different objectives.

Step 6. Termination conditions. If , then stop the process. Otherwise, increase generation () and go to step 3.

Experiment and analysis

This paper uses MATLAB R2018a to simulate the data experiments in Windows 10 i7-1.99 GHz 8 GB 64 -bit operating system.

Algorithm contrast

According to the instance generation framework proposed in research (Cheng et al., 2020, Deb et al., 2002, Dorling et al., 2017), we generate a set of instances with the number of customers (EPNs) in the interval [10,80], and each group of EPNs corresponds to 2–9 candidate CHs respectively. Considering the payload of the drone and the characteristics of medical supplies, each EPN is given a uniform random delivery demand of 0.5–2 kg and a random pickup demand of 0.1–1.5 kg. The delivery demands of 40% of EPNs are drawn uniformly from [0.5, 1.3], and the delivery demands of the remaining EPNs are drawn uniformly from [0.5, 2]. The pickup demands of 60% of EPNs are drawn uniformly from [0.1, 0.7], and the pickup demands of the remaining EPNs are drawn uniformly from [0.1, 1.5]. Regarding the scope of drone delivery, we take the H4 four-rotor drone with a maximum range of 15 km as an example (SF-Technology, 2021), and its cruise speed is 12 m/s. Based on the test of Dorling et al. (2017), we set ,,,,.

To verify the performance of the M−NSGA−II for solving the proposed BOMT-DLRP-SPD, the results of M−NSGA−II are compared with the results of the MOEA/DD and I-DBEA which are multi-objective solving algorithms from Tian et al. (2017). Algorithm performance metrics include the minimum value of the objective functions (the minimum operating cost , the minimum travel time ), the number of Pareto solutions , the running time , and the comprehensive performance metric—Hypervolume (HV) (Zitzler and Thiele, 1999) that reflects the convergence and distribution of the algorithm. The results of the three algorithms to solve different instances are shown in Table 2 .

Table 2

Summary results of three algorithms to solve different size instances.

Instance	M−NSGA−II					MOEA/DD
	C.Opt($)	T.Opt(s)	N.Prt	R.T(s)	HV	C.Opt($)	T.Opt(s)	N.Prt	R.T(s)	HV
EPC-10	8496.99	1280.82	14	13.45	0.574	9052.52	1333.99	10	28.65	0.552
EPC-20	17522.59	2777.71	14	19.49	0.502	17612.67	2892.35	12	34.64	0.491
EPC-30	22539.34	3868.98	25	29.11	0.472	23485.18	4003.17	14	49.66	0.384
EPC-40	36673.63	5069.58	15	37.12	0.533	37934.38	5550.60	4	53.09	0.413
EPC-50	38257.70	6164.75	18	36.27	0.516	41152.55	6707.64	6	47.39	0.426
EPC-60	54973.81	7762.56	17	47.82	0.456	56294.13	8541.28	6	66.56	0.385
EPC-70	69137.33	9395.86	13	51.30	0.422	72492.24	9902.16	3	71.59	0.358
EPC-80	75310.29	10026.07	14	58.42	0.468	79349.85	11178.87	9	81.00	0.390
Average	40363.96	5793.29	16	36.62	0.493	42171.69	6263.76	8	54.07	0.425
t-test	–	–	–	–	–	−3.57	−3.55	5.97	−13.17	5.34
p-value	–	–	–	–	–	9.09E-03	9.29E-03	5.59E-04	3.39E-06	1.08E-03

As can be seen from Table 2, for the BOMT-DLRP-SPD proposed in this paper, the solution performance of M−NSGA−II is better than MOEA/DD and I-DBEA under the same computing environment. It is mainly manifested in five aspects: First, the convergence and distribution of the population obtained by M−NSGA−II is better than MOEA/DD and I-DBEA. In the results of the eight instances, the HV values of M−NSGA−II are greater than the results of the other two algorithms, especially when the size of the instance increases and the advantage of M−NSGA−II is more prominent. Second, M−NSGA−II has the shortest running time among the three algorithms, and this conclusion will not change with the increase of the instance size. Third, in acquiring Pareto solutions, M−NSGA−II can always obtain a Pareto solution set containing more solutions. Fourth, our model expects to obtain a solution with lower operating costs and shorter travel time. According to the minimum values of the objective functions obtained by the three algorithms, the results of M−NSGA−II have lower minimum operating costs and minimum travel time. Fifth, according to the statistical analysis of the two sets of results (M−NSGA−II vs. MOEA/DD, and M−NSGA−II vs. I-DBEA), the calculated t-test and p-value confirmed a significantly irrelevant relationship between the groups. Therefore, from the perspective of the basic metrics and comprehensive performance metrics of the algorithms, M−NSGA−II is conducive to obtaining better results in a reasonable time when solving the BOMT-DLRP-SPD model.

Case study and model analysis

Based on the instance EPN-50 described in Section 5.1, this section further analyses the impact of model features and parameter changes.

Results of the case study

There are 36 Pareto solutions in the Pareto solution set with the most considerable HV value of the instance EPN-50. The optimal solutions in the Pareto solution set are all feasible. It is essential to select a final delivery plan in the set by an arbitrary method (Bortolini et al., 2016, Lu et al., 2012). As a result, we introduce Eq. (37) to converge to such a final solution:

(37)where n is the index of the n-th solution laying on the Pareto front. and are the single objective optimal solutions for operating cost and travel time, respectively. The solution with minimum can be regarded as a delivery plan that is ultimately used. It is worth noting that we provide a selection method that reasonably considers the two objective functions. In practical application, different objective weights can be set according to the emergency of medical supplies delivery. For example, the solution with the least travel time is adopted because the safety of life and supplies is the most important in an emergency. Fig. 5 a depicts the distribution of 36 Pareto solutions on the Cartesian coordinate system and the detailed plan of the final solution selected according to formula (37).

Fig. 5

Schematic diagram of the results.

According to the detailed plan shown in Fig. 5 b, CH1 is not involved in the delivery, which means that there is no need to deploy drones and related supporting equipment in CH1. 14 drones are required to complete this batch of delivery tasks, of which the drone-10 in CH6 performs two trips. Besides that, we have compared the results under the scenarios of simultaneous pickup and delivery (SPD) and separate pickup (P) and delivery (D), which aims to illustrate the importance of considering simultaneous pickup and delivery to improve delivery efficiency and ensure the timeliness and safety of medical supplies. The results are displayed in Table 3 .

Table 3

Results of simultaneous delivery and pickup vs separate delivery and pickup.

Result	Condition
	SPD	D	P	D + P	Δvs.SPD
Cost ($)	37299.73	37311.62	38864.45	76176.07	104.23%
Time (s)	6436.76	6410.25	6520.83	12931.08	100.89%
Number of drones	14	14	15	29	107.14%
Travel distance (m)	51137.10	50221.43	53924.75	104146.18	103.66%
Energy consumption (kWh)	6.20	7.50	6.67	14.17	128.71%

The only conclusion we can reach from Table 3 is that the simultaneous pickup and delivery mode has better timeliness, safety and resource saving performance than the separate pickup and delivery modes. In terms of optimization objectives, separate pickup and delivery require 104.23% more operating costs than simultaneous pickup and delivery, and it increases the travel time of the drones by 100.89%. In addition, separate pickup and delivery increase the travel distance of the drones by 103.66% and consume 128.71% more energy than simultaneous pickup and delivery. If the delivery and pickup tasks are performed separately, 107.14% more drones need to be deployed than simultaneous pickup and delivery, which causes idle resources. More importantly, it is undesirable in major public health emergencies that require reduced contact. Even if the drones are reused for pickup and delivery, the number of deployed drones is 15. However, this may not meet the time requirements for delivering medical supplies and samples.

Sensitivity analysis

This section attempts to carry out sensitivity analysis from the two aspects of drone cruise speed and battery energy capacity to observe the impact of model parameters change on the results of the entire delivery plan (including total operating cost, travel time, the number of drones used, and total travel distance). The solutions of the bi-objective optimization are all selected from the Pareto solution sets according to Eq. (37).

Fig. 6 reveals the impact of changes in drone cruise speed, which shows that the total operating cost, total travel time, the number of drones used, and the travel distance of the delivery plan have all declined as the increase of drone speed. It can be seen from Fig. 6 a that when the drone speed increases from 6 m/s to 21 m/s, the total cost and travel distance have a decreasing trend, but the changing trend is not apparent. As shown in Fig. 6 b, when the drone speed is increased from 6 m/s to 21 m/s, the total time of the drone is reduced by 98.13%, and the number of drones used is decreased by 33.33%.

Fig. 6

Impact of changes in drone cruise speed.

Fig. 7 displays the impact of changes in the battery energy capacity of drone. There is no doubt that the size of the drone energy capacity has a significant impact on the results of the delivery plan. Fig. 7 a indicates that with the increase of energy capacity, the delivery plan's total cost and travel distance show an apparent decreasing trend. Comparing the results of energy capacity at 0.3kWh and 1.5kWh, it is found that the total cost and the travel distance have decreased by 43.19% and 45.75%, respectively. Fig. 7 b demonstrates that as the energy capacity increases, the number of drones used decreases and the average travel time per drone trip increases. When the energy capacity increases from 0.3kWh to 1.5kWh, the number of drones used is reduced from 19 to 12, there is a reduction ratio of 58.33%, and the average travel time per drone trip increases by 38.04%.

Fig. 7

Impact of changes in drone energy capacity.

In short, the data results manifest that the influence of drone energy capacity changes on the delivery plan results is more apparent than the influence of drone cruise speed. The size of the battery energy capacity of the drone is a pivotal factor in determining the number of drones used. In practical applications, if drones cannot consider both high battery energy capacity and high cruise speed, managers should use drones with different attributes according to different delivery environments. For example, when the outbreak of a pandemic such as the COVID-19, in urban communities where social distance is strictly restricted, drones with large battery energy capacity and strength endurance are more suitable for delivering medical supplies in large demand such as protective clothing, gloves, and medicines. Because this is conducive to reducing the number of drones used and operating costs. However, when faced with emergency treatment needs of epidemic prevention nodes such as isolation warehouses and community hospitals, medical supplies such as blood, vaccines, and samples may have strict requirements for delivery time. At this time, drones with high cruise speed will be a more sensible choice. Because this can contribute to significantly reducing the travel distance and delivery time, which will come at a particular cost.

Conclusions and future research

The application of drones reduces the use of human resources in complex environments, thereby decreasing the infection rate. In major public health emergencies that threaten the safety of cities like the spread of the COVID-19 pandemic, social distance restrictions affect the timeliness of the medical supplies delivery, which directly imperils the people's lives. To improve the efficiency of medical supplies delivery and test sample recovery, a bi-objective MIP model, denoted as BOMT-DLRP-SPD, is proposed for the defined MT-DLRP-SPD to construct a drone-based urban aerial transport bridge for medical supplies. M−NSGA−II is used to solve the proposed MIP model. The results of multiple data instances indicate that M−NSGA−II has better solution performance than MOEA/DD and I-DBEA. Through the sensitivity analysis of the cruise speed and battery energy capacity of the drone, we find that: (1) The increase in the cruise speed of the drone leads to a reduction in the operating costs, total travel time, number of drones used, and travel distance of the final solution, among which the changing trend of total travel time is the most obvious. (2) The increase in energy capacity of drone results in the reduction of the operating costs, the number of drones used, and the travel distance of the final solution, among which the change in the number of drones used is the most obvious. The average travel time of each drone path gradually increases as the energy capacity increases. The above findings can support managers who use drones to pick up and deliver medical supplies.

The BOMT-DLRP-SPD with simultaneous pickup and delivery mode has better time, safety, and resource-saving performance compared with the separate pickup and delivery modes. To better match the transportation environment and reduce the resource waste, the BOMT-DLRP-SPD also considers the multi-trip optimization of drones and the joint optimization of location and routing problems under the constraints of drones' payload and energy capacity. Meanwhile, with the dual time constraints of medical supplies delivery and test samples recovery, the model takes operating cost and the travel time of drones as the optimization objectives to output delivery plans that balance economic and time benefits. In the decision-making process, if managers have strict requirements on transportation time, the solution with the least travel time can be selected in the Pareto solution set. In an epidemic prevention normalization sitting, managers must consider cost as an optimization objective. If cost impact is ignored, enough drones will be dispatched to complete the mission regardless of cost. We hold that jointly optimizing total operating cost and travel time is a meaningful consideration under the premise of meeting emergency demands.

From an academic research perspective, this type of operations research technique and MIP model can be extended for the research of VRPSPD. An extension is VRPSPD with demand splitting (Archetti et al., 2011, Maini and Sujit, 2015). The pickup and delivery need of some or all nodes in the logistics network may be so large that drones cannot satisfy the node demand in a single visit. Then it is worthwhile to study how to solve such problems by splitting and distributing the demand of nodes and allowing multiple visits to the nodes. In addition, the problem of fair relief distribution also cannot be ignored. In humanitarian logistics, when resources are scarce and in shorter supply than the needs, the overall degree of demands satisfaction should be considered, and a the relief supplies should be distributed to demand points as equally as possible (Anaya-Arenas et al., 2018, Gutjahr and Fischer, 2018). Taking hybrid approaches, such as combining constraint programming and MIP, exploring VRPSPD with distribution fairness is an interesting attempt (Ham, 2018). Furthermore, any disruption to the delivery network can harm material delivery (Dehghan et al., 2021). It is crucial to construct an MIP under disturbance management strategy, which will effectively mitigate the impact of disturbance factors on network stability.

The designed MIP can be extended to solve optimization problems in various fields in practical applications. A similar application to the studied drone delivery problem is the involvement of unmanned vehicles in last-mile logistics. In the future, with the advancement of autonomous driving technology, large-capacity unmanned vehicles can be widely used to deliver express parcels and pick up recyclable packaging boxes (Reed et al., 2022). The last extension could consider a future freight underground logistics system (ULS) (Fan et al., 2020). The system requires an effective logistics network design for future demand growth. ULS is characterized by multi-depot, high-cost, integrated pickup and delivery, and multi-modal transportation. Nevertheless, it cannot replace ground logistics in the short run. Therefore, how to utilize ULS for bulk delivery and then combine it with ground transportation for pickup is another interesting direction for future research.

CRediT authorship contribution statement

Yuhe Shi: Writing – original draft, Software, Data curation. Yun Lin: Methodology, Conceptualization. Bo Li: Investigation, Visualization. Rita Yi Man Li: Writing – review & editing.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.