Multivisit Drone-Vehicle Routing Problem with Simultaneous Pickup and Delivery considering No-Fly Zones

,


Introduction
For a long time, the "last mile" has consistently been the bottleneck in the efciency of the distribution process. Improving the operational efciency of the "last mile" is crucial to logistics companies and even society. Currently, in the "last-mile" distribution, trucks often have low distribution efciency due to urban trafc congestion or rugged mountain roads, which pose signifcant challenges to logistics companies. However, the application of drones in the distribution feld has the advantages of not being afected by ground obstacles, fast delivery, and low cost; hence, it has attracted increasing attention from the industry and research scholars. Although the load capacity of drones is small and battery technology cannot support long-range delivery, the joint delivery of trucks (ground vehicles) and drones can form a novel and efcient delivery mode [1] which has led to research on the drone-vehicle routing problem (DVRP). Furthermore, introducing pickup services into this novel combination model could ofer a practical and promising last-mile solution, given the growing demand and practical importance of reverse logistics.
Te drone-vehicle routing problem with simultaneous pickup and delivery (DVRPSPD) is a complex route planning problem. Tis problem aims to meet the pickup and delivery requirements of a set of customers while minimizing the route cost. Te characteristics of this problem are as follows. First, joint delivery by trucks and drones is considered. Te truck feets start and end at the depot, and the drones are launched and retrieved from the depot or the truck they are transported to. Te trucks are responsible for long-distance pickup and delivery tasks, while the drones handle the same short-distance tasks. Tis service method can greatly reduce labor and time costs while enhancing service efciency. Additionally, it considers simultaneous pickup and delivery for meeting the customer needs. Unlike traditional parcel delivery, this method enables trucks and drones to provide pickup and delivery services during the same trip, which greatly reduces delivery time and improves delivery efciency.
In 2013, Amazon's drone delivery plan marked the earliest instance of drone involvement in delivery services. Tree years later, the Prime Air drone by Amazon successfully delivered its frst order. In subsequent tests, drone delivery services have shown to be efective in certain scenarios. For instance, in Shenzhen, China, drones have been utilized to deliver goods and food from stores or restaurants to customers, ofering residents near business districts a novel "3 km 15 min" service experience in response to the boom in instant retail [2]. Furthermore, some companies began exploring synchronous collaboration systems that combine trucks and drones. Mercedes-Benz's Vision Van, for instance, includes a fully automated cargo space and two rotor drones for autonomous delivery.
While the purpose of these tests was to address the technical hindrances associated with drone delivery services, certain researchers have delved into the issue from an operational standpoint. For example, Kim et al. [3] placed emphasis on the planning of drone-assisted healthcare services, which involves drones being utilized for delivering medicine to patients and collecting samples such as blood or urine tests. Ham [4] extended the parallel drone scheduling traveling salesman problem (PDSTSP), considering the drone's cargo and pickup functions. Tis means that the drone can either return to the depot after completing a delivery or directly move to another customer to pick up goods. Wikarek et al. [5] also explored the capacitated vehicle routing problem with drones, which involves drones delivering packages to customers while also retrieving packages from them using multiple mobile depots to launch and retrieve drones. However, their study only used drones for pickup and delivery independently rather than in synchronous conjunction with trucks. Current research on combinatorial problems primarily focuses on delivery services [6][7][8][9][10][11][12][13]. In light of this, Zhang and Li [14] introduced the pickup and delivery service into the collaborative distribution system of trucks and drones to maximize drone usage. However, their study does not take into account the capacity constraint of the truck and restricts the launch and retrieve node of the drone to be the same node, which simplifes the problem and increases the waiting time. Meng et al. [15] broke this limitation and extended the problem to the vehicle routing problem with simultaneous pickup and delivery (VRPSPD) for the study.
Despite the advantages of using drones for delivery and pickup operations, implementing such services is still in its infancy and limited to specifc scenarios. Tere are some obstacles to wider adoption, as numerous cities have instituted no-fy zones to ensure the safety of drone deliveries and protect the public interest. Tese no-fy zones restrict drones from entering sensitive areas, such as airports, nuclear power plants, and government buildings, to avoid potential security threats, infringement of the legal rights of others, and interference with the normal operation of public facilities. In the joint distribution path planning problem of drones and trucks, the existence of no-fy zones limits the movement range and path selection of drones and increases the complexity and difculty of the problem, which is suitable for the actual drone delivery scene. Terefore, the existence of no-fy zones in the joint delivery of trucks and drones is necessary.
To increase the utilization rate of drones, some technically feasible operations can be considered, such as allowing drones to carry multiple packages per trip and providing simultaneous pickup and delivery services, not just delivery, which can increase the utility of the drone. Currently, there is no literature that considers the existence of no-fy zones in the study of the DVRPSPD. Terefore, this study considers the impact of no-fy zones on path planning based on the DVRPSPD and allows drones to provide pickup and delivery services for multiple customers within their capacity and endurance. Tis problem is known as the multivisits drone-vehicle routing problem with simultaneous pickup and delivery considering no-fy zones (MDVRPSPDNF). We built a discrete optimization model with the goal of minimizing the route cost, aiming to provide theoretical support for improving the logistics company's response to customers' various needs and the normal delivery of drones. Furthermore, we propose a two-stage heuristic algorithm based on the simulated annealing (SA) algorithm, which includes a path encoding method that can represent the branch structure, an adjustment method for infeasible solutions, and an algorithm acceleration strategy. Tis algorithm is fast and efective in solving the MDVRPSPDNF problem. Finally, we observe through computational experiments that a combined truck-drone system allowing multiple visits per trip saves signifcant delivery costs. Te results of these numerical experiments indicate that logistics enterprises would not only achieve cost reduction and increase efciency from route optimization but also consider the rationality of fexible and efective distribution vehicles and resource allocation.

Related Literature
Tis section provides a synopsis of the pertinent literature that previously explored concepts in relation to the research presented in this study. For a more extensive summary and survey, we direct the reader to studies conducted by Otto et al. [16], Khouf et al. [17], Macrina et al. [18], and Madani and Ndiaye [19].
Murray and Chu [6] initially proposed the fying sidekick traveling salesman problem (FSTSP), which serves as a cornerstone for the joint truck-drone route optimization problem, garnering considerable attention from both industry and academia. In this problem, a truck carries a drone for delivery purposes, with each drone fight serving only one customer. Based on FSTSP, Ha et al. [7] solved the FSTSP problem by modifying the objective function from minimizing delivery time to minimizing delivery costs. Tey proposed two heuristic methods-TSP-local search (TSP-LS) and greedy randomized adaptive search procedure (GRASP)-to address this issue. Agatz et al. [8] proposed the traveling salesman problem with drone (TSPD), where they set the drone retrieve node limit and the drone launch node limit to the same node, and designed a hybrid algorithm that combines local search and dynamic programming to solve the problem. Later, Mario et al. [20] explored the possibility of along-route operations in the TSPD, allowing the drone to serve a broader area. Chang and Lee [21] highlighted the importance of fnding a better route based on the TSPD to achieve a larger drone delivery area. Tey established a nonlinear planning model and demonstrated its efectiveness. Yurek and Ozmutlu [22] proposed a decomposition-based iterative algorithm for moderately sized TSPD problems. Kim and Moon [23] introduced drone stations to the TSPD, creating the traveling salesman problem with drone stations. Wang et al. [24] examined the TSPD from a manager's perspective, with biobjectives of operating cost and completion time, and proposed an improved nondominated sorting genetic algorithm to solve the problem. Additionally, several studies focused on exploring precise methods for the FSTSP/TSPD [25][26][27][28] and heuristic algorithms [29,30].
Kitjacharoenchai et al. [9] extended the FSTSP problem, proposing a new variant of a truck carrying multiple drones. To this end, they established a mixed integer programming (MIP) model and heuristic algorithm. Murray and Raj [31] conducted an in-depth study of the multiple fying sidekicks traveling salesman problem (MFSTSP), considering the energy consumption model of drones with respect to the package weight, speed, and operation time. Tey proposed a three-stage iterative heuristic algorithm that can handle up to 100 customer instances. Furthermore, Peng et al. [32] proposed a truck-assisted multidrone package delivery method and provided a hybrid genetic algorithm to solve the problem. Freitas and Penna [33] proposed a novel MFSTSP for parcel delivery services using parking-assisted truckdrones as an alternative to FSTSP. Dell'Amico et al. [34], Gavani et al. [35], Lu et al. [36], and Tinic et al. [37] studied exact and heuristic algorithms to solve MFSTSP problems.
Gonzalez et al. [10] expanded on the FSTSP problem, also known as the truck-drone team logistics problem, by considering the number of customers a drone visits during a single trip. Tis approach allows drones to visit multiple customers, which reduces delivery costs. Luo et al. [38] delved further into the issue with their study of the multivisit traveling salesman problem with multidrones, considering drone energy consumption factors such as fight time, selfweight, and payload. Tey proposed a multistart tabu search algorithm to solve the problem for up to 100 customers. Windras Mara et al. [39] introduced a new mathematical formulation and a heuristic algorithm based on adaptive large-scale neighborhood search (ALNS) techniques to optimize the route for combined systems. Mahmoudi and Eshghi [40] explored the energy-constrained multivisit  traveling salesman problem for multiple drones at   noncustomer rendezvous locations, providing a mathematical model and heuristic algorithm for large-scale examples. Wang et al. [11] expanded the single-truck feet problem to encompass multiple trucks and designated it as the vehicle routing problem with drone (VRPD) to better ft the practical distribution scenarios of enterprises. Masmoudi et al. [12] analyzed the operational expenses of VRPD and devised an ALNS algorithm to handle issues with up to 250 customers. Ulmer and Tomas [13] centered their VRPD research on the same-day delivery. Wang and Sheu [41], Tamke and Buscher [42], and Zhen et al. [43] each employed exact algorithms to solve the VRPD model. Euchi and Sadok [44] proposed a hybrid genetic sweeping algorithm for addressing the VRPD problem involving up to 200 customers. Lei et al. [45] adopted a novel dynamic artifcial bee colony algorithm to minimize overall operating costs for the VRPD. Schermer et al. [46] broadened the VRPD to include multiple drones per truck to serve all customers. Masmoudi et al. [12] extended the VRPD to a feet of drones equipped with payload bays for serving more customers during a single trip. In Kuo et al. [47] study, a biobjective mathematical model of the VRPD was established to minimize delivery completion time and carbon emissions. Tey proposed a nondominated sorting-based genetic algorithm to solve this problem.
Kuo et al. [48] explored the VRPD with customer time windows by minimizing travel costs in their study that integrates multiple factors. Tey developed a variable neighborhood search algorithm to solve a ffty-customer problem. Dayarian et al. [49] proposed a delivery service based on the VRPD, where trucks meet the needs of customers who placed orders, while drones fulfll newly emerging dynamic needs. Kim et al. [3] focused on droneassisted medical services that deliver medicine to patients and retrieve test kits such as blood or urine samples. Ham [4] extended the PDSTSP to include the delivery and pickup of the drone. Wikarek et al. [5] studied the capacity vehicle routing problem with drones that not only deliver but also accept packages from customers and used several mobile depots for drone launch and retrieval. Lu et al. [50] performed an advantageous study on the cooperative distribution system of trucks and drones in the multiobjective humanitarian routing problem. In their study, customers are defned as demand points or replenishment points. In contrast, in Zhang and Li's study [14], the customer is both a demand point and a replenishment point, and it is restricted that the launch and retrieve node of the drone must be the same. But their study did not take into account the capacity constraints of trucks. However, Meng et al. [15] allowed the launch and retrieve node of the drone to be diferent and extended this problem to the VRPSPD for research. Jeong et al. [51] considered two critical practical issues, namely, the impact of package weight on drone energy consumption and no-fy zones, and extended the previous vehicle routing model to hybrid delivery systems.
However, there are no VRPD studies that simultaneously addressed pickup and delivery in no-fy zones. Numerous cities in the world have designated no-fy zones, preventing drones from visiting customers within them. Tis imposes Discrete Dynamics in Nature and Society greater demands and obstacles on the route planning of both trucks and drones. Our approach allows each drone to serve multiple customers on a single trip, making joint route optimization between trucks and drones more practical. Hence, the fexibility of the MDVRPSPDNF yields new practical pickup and delivery models. Recent research results in this feld are summarized in Table 1.
T represents the number of trucks, D represents the number of drones carried on each truck or the number of drones when serving customers independently with the truck, and m is the abbreviation of "Multiple." SV/MV represents the number of customers per trip (i.e., single visit or multiple visits). SPD represents whether the problem considers both demand for pickup and delivery at the same time. NF represents whether the problem considers no-fy zones. Te above literature suggests that the joint distribution problem of drones and trucks can be regarded as adding drones to the VRP. Terefore, there are also numerous types of research worth learning in this regard [52,53].

Problem Description.
We examine the scenarios of forward and reverse terminal logistics distribution in the presence of no-fy zones as the subject of research. Te distribution route is optimized through the use of a cooperative distribution mode of trucks equipped with drones, as shown in Figure 1. Tis optimization aims to achieve wider distribution and more efcient package delivery.
Te problem is defned as follows: there is a depot and multiple customers, some of whom are situated within no-fy zones. Te demand and recycling volume of each customer, as well as the distance between the depot and the customers and between any two customers, is known. Te depot possesses an unlimited number of homogeneous trucks and drones. A single drone can transport goods to multiple customers within its maximum loading capacity and farthest fying distance. Customers can be serviced only once by one truck or one drone. Te depot dispatches several trucks and drones to transport goods to customers. Te trucks are sent out from the depot, each carrying a drone, and the cargo load of the truck must not exceed its maximum loading capacity. Te drone can be launched and retrieved from the depot or the truck carrying it, and it returns to the depot or its truck after delivering and receiving goods to several customers. Te objective of this study is to devise a routing scheme for truck and drone services that minimize the routing cost.
MDVRPSPDNF can be described as an edge selection problem in a directed and connected network graph. Let G � (V, A) be a directed graph, where V � 0, 1, 2, · · · , n, n + 1 { } represents the node set, node 0, n + 1 represents the depot, the rest of the nodes represent customers, and A represents the arc set, A � (i, j) | i ∈ V, j ∈ V, i ≠ j . On directed graph G, a reasonable delivery route must start from node 0 and end at node n + 1. C � V\ 0, n + 1 { } represents the customer set, K � 1, 2, · · · , k { } represents the set of delivery trucks, and P � 1, 2, · · · , p represents the set of drone trips. F � 1, 2, · · · , f refers to the gathering of customers in nofy zones. According to the customer's demand and the weight of the recycled volume, all customers are divided into two categories: customer set only served by trucks, and customer set C U : C U � C/C T served by both trucks and drones. In set C U , when customer i's demand q i and recycling volume r i do not exceed the maximum loading capacity Q u of the drone and the travel distance after adding customer i does not exceed its farthest fying distance D max , drone delivery is preferred. If the above conditions are met but the route cost of using drones for distribution is higher than that of using trucks, priority is given to trucks for delivery.
Te parameters and decision variables involved in the MDVRPSPDNF model are shown in Table 2.
Before building a mathematical model, some basic assumptions must be established: (i) If either the truck or the drone arrives earlier, they must wait for the other to arrive before restarting the pickup/delivery process (ii) Because the drone fies almost in a straight line, Euclidean distance is used for calculation; the truck is limited by the road network such that it is assumed that the distance of the truck is d t ij � δ · d u ij , where δ is a constant, which represents the road tortuosity (iii) Drones are not allowed to enter no-fy zones (iv) When the drone is delivering, the truck proceeds to the next customer node to complete the pickup and delivery needs (v) Te battery is replaced immediately after the drone meets the truck to ensure endurance in the next drone fight (vi) Te truck and drone maintain a constant speed during the delivery process

Mathematical Model.
Based on the mathematical model constructed by Luo et al. [38] and combined with the specifc characteristics of MDVRPSPDNF, the following model is obtained: Discrete Dynamics in Nature and Society

Route Constraints.
Route constraints are given as follows: A collection of customers that can be served by trucks and drones k∈K p∈P Recycling volume of customer i ∈ C Q t /Q u Maximum loading capacity of trucks/drones Customer i ∈ C is in the no-fy zones L 0k Load of truck k ∈ K when it leaves the depot Loading capacity of truck/drone after service to customer i ∈ C u i /d i Truck/drone access sequence for node i ∈ V 6 Discrete Dynamics in Nature and Society Te objective function (1) endeavors to minimize the cost of the route, which is the summation of the truck and drone route costs. Constraint (2) afrms that each customer can only receive a single service. Constraint (3) assures that customers served by trucks are exclusively served by them due to load capacity restrictions. Constraint (4) ensures that customers serviced by trucks are directly served by them. Constraint (5) ensures that customers requiring drone services are only reachable through drone trips. Every drone trip comprises an open route with a launch and retrieve node, indicating the commencement and conclusion of the trip, respectively. Te launch node of a drone has solely an outdegree, and the retrieve node has solely an indegree. Constraint (6) ensures that all nodes visited by the drone trip p have a fow balance. Constraints (7) and (8) limit the total outdegree and indegree of each node. Constraint (9) maintains the fow balance of the truck and ensures that the truck can service each customer at most once. Constraint (10) necessitates that each edge is traveled at most once during a drone trip. Constraint (11) indicates that each truck can be used at most once, starting and ending from the depot. Constraint (12) mandates that for each edge visited by the drone, at least one customer must be serviced on p to prevent the drone from fying unnecessarily. Constraints (13) and (14) designate the launch and retrieve nodes for the drone's trip p. Constraints (15)- (16) require that the launch and retrieval nodes be visited by trucks. Constraint (17) guarantees that the launch and retrieve nodes of a drone trip p difer. Constraint (18) ensures that an equal number of launch and retrieve nodes are present per drone trip and that there is at most one launch node per trip. Constraint (19) indicates that each truck can only be used after leaving the depot. Constraint (20) indicates that drones are not allowed to independently complete tasks. Constraints (21) and (22) represent the subtour elimination of trucks and drones, respectively.

Capacity Constraints.
Capacity constraints are given as follows: Discrete Dynamics in Nature and Society Constraint (23) is utilized to compute the primary load of the designated truck k situated at the depot. Constraint (24) represents the maximum loading limit that the truck is capable of withstanding, while constraint (25) represents the maximum loading limit of the drone. Constraint (26) manifests as the formula for calculating the loading capacity of truck k after the initial customer has been serviced along the route. In the same vein, constraint (27) corresponds to the formula for determining the loading capacity of truck k after any other customer (excluding the frst) has been serviced along the way. Similarly, constraints (28) and (29) are used to calculate the loading capacity of the drone during the trip. Constraints (30) and (31) indicate that the loading capacity of the truck and the drone after serving any customer on the route must be equal to or greater than its maximum loading capacity.

No-Fly Zones and Farthest Flying Distance Constraints.
In the domain of logistics distribution, the incorporation of a drone can enhance distribution efectiveness and concurrently reduce costs. To ensure the safety of drone fights, a crucial approach entails the prudent delineation of obstacle avoidance routes and fight routes during the fight control of drones. Hence, it is of utmost importance to devise a rationalized fight plan and route planning while also considering the no-fy zones and the transportation requirements of drones.
Te current set of no-fy zones for drones includes certain designated areas, such as airports, military bases, prisons, schools, and government properties, that permanently prohibit any sort of fying activity, alongside temporary no-fy zones that are caused by unforeseen emergencies such as forest fres. In light of this, these no-fy zones must be considered when designing the fight route, and the drone must not enter these zones.
Tis study delineates the no-fy zones as an area where drones are strictly prohibited from fying. Within this domain, drones are neither permitted to launch nor allowed to enter from other areas. Figure 2 highlights three distinct prohibited routes for drones traversing no-fy zones: (a) customers in no-fy zones cannot be used as the launch node of the drone; (b) customers within no-fy zones cannot be served by drones; and (c) customers in no-fy zones cannot be used as the retrieval node of the drone.
In summary, this study introduces a series of constraints to ensure the rationality and accuracy of the path planning model design and restrict the drone from passing through no-fy zones.
i∈V j∈V 8 Discrete Dynamics in Nature and Society Constraint (32) means that when customer i is in a no-fy zone, its corresponding binary variable must be 1. Constraints (33) and (34) indicate that if the drone must retrieve or launch, the retrieve/launch node must be set at a location out of the no-fy zone. Constraints (35) and (36) mean that if the path transfers from node i to node j and either node i or node j is in the no-fy zone, the decision variable y kp ij must be 0. Constraint (37) denotes the farthest fying distance constraint of the drone. Constraints (38) -(40) represent the value range of variables.

Solution Approach
As a variant of the VRP, the MDVRPSPDNF is likewise a problem that belongs to the NP-hard class. Te traditional optimization method is difcult to solve; hence, we choose the SA algorithm with a simple structure that is capable of jumping out of the local optimal solution to solve it. Nevertheless, traditional SA algorithms are not capable of optimizing the truck and drone routes simultaneously. Terefore, a two-stage heuristic algorithm, based on the SA algorithm, is proposed in this study. Te algorithm follows the "Route First, Cluster Second" concept. First, the problem is treated as a VRPSPD, and the SA algorithm is employed to devise the entire route. Ten, customers are categorized based on their demand and recycling volume. Customers in set C U are then excluded from the route, and drones are used to serve them. For every drone service, the minimum cost route is calculated, which satisfes the maximum loading capacity and the farthest fying distance constraints of the drone. Te search process permits the drone to cater to multiple customers while adhering to the maximum loading capacity and farthest fying distance constraints. Moreover, an adjustment mechanism for infeasible solutions is devised to comply with the no-fy zone constraints. To enhance the iteration speed of the algorithm, an acceleration technique is proposed. Tis approach greatly reduces the complexity of the calculation problem and improves the efciency of the algorithm search. Te specifc algorithm is elaborated as follows.

Te First Stage: Simulated
Annealing Algorithm to Solve the VRPSPD. Te central idea of the SA algorithm is to obtain the optimal solution by performing several neighborhood operations on the current solution. Compared with other optimization algorithms, SA allows the acceptance of a solution that is worse than the current one with a specifc probability during the search process. Tis attribute enables the algorithm to avoid being trapped in the local optimum. When employing the SA algorithm to solve the VRPSPD, three vital stages must be incorporated: coding and objective function, neighborhood structure, and acceptance criteria and annealing.

Coding and Objective Function.
In this study, the depot and customers are coded in the solution simultaneously, i.e., the depot is inserted into the solution in the form of a number greater than the number n of customers. Te length of the solution is n + m − 1, where n is the number of customers and m is the number of trucks.
Considering the constraints of loading capacity, we adopt the method of imposing penalties on routes, violating the constraints to make each segmented route meet the constraints of loading capacity. Terefore, the objective function is as follows: where s is the delivery plan, f(s) is the value of the objective function, t(s) is the total distance cost of the vehicle, c(s) is the sum of the violations of the loading constraints when the truck leaves each node on each route, and α is the weight of the violation of the loading constraints.

Acceptance Criteria and Annealing.
Te core concept of SA is to accept a solution worse than the current solution with a certain probability. Te formula for the probability of accepting a new solution is as follows: where f(S curr ) represents the objective function value of the current solution and f(S new ) represents the objective function value of the new solution.
As the number of iterations increases, the probability of acceptance must also decrease. Tis is because the search time must be reduced, i.e., the temperature decreases with the increase in the number of iterations, and the formula is as follows: where β is the cooling factor and 0 < β < 1.

Te Second Stage: Routing Algorithms for Truck-Drone.
Te multivisit drone-vehicle routing problem for pickup and delivery conundrum encompasses three distinct decisions. Te frst step is to determine which customers are served by trucks and drones, respectively. Te subsequent decision involves determining the order in which the truck will visit the customers. Te third decision involves designating the sequence in which each drone will visit the customers as well as deciding upon launch and retrieve nodes. Tese nodes can either be the depot or any customer location for the truck service. Imperatively, these decisions have a consequential impact on one another and therefore must be made in tandem. Te ensuing discourse delves into each of these three decisions. Figure 3(a), along with its corresponding encoding solution is depicted in Figure 3(b). Te encode comprises two components: the frst part, denoted as Part 1, represents the order in which the customers are visited by the truck. Tis section commences and concludes with 0, which denotes the depot. Te digits between two zeros indicate the customer visitation sequence for the truck service. Te second part, referred to as Part 2, represents the sequence of customers visited by the drone. Te frst and last digits in this section signify the launch and retrieve nodes for each drone trip, respectively. Te k-th row pertains to the route taken by the k-th truck, while the p-th trip of the drone carried by the k-th truck is represented by the drone route.

Truck-Drone Routing Representation. Te combined route for a truck and drone is depicted in
To facilitate ease of handling, solutions are presented as multidimensional arrays. Given n customers, k trucks, and k drones, a single drone may undertake a maximum of p trips, and the solution is coded in kq dimensions. q � max k p + 1, where max k p represents the maximum number of trips among k drones. Tis arrangement allows for the maximum number of trips made by the truck and drone. Te decoded route shown in Figure 3(b) comprises the following: the route taken by the frst truck is (0⟶1⟶3⟶4⟶0), while the frst trip of the drone carried by the frst truck is (1⟶2⟶3), and the second trip of the drone carried by the frst truck is (4⟶5⟶6⟶0); the route taken by the second truck is (0⟶7⟶8⟶10⟶0), and the frst trip of the drone carried by the second truck is (8⟶9⟶10). Evidently, the coding of multidimensional arrays lucidly demonstrates the collaborative distribution route of multiple trucks and drones. Tis coding methodology can also be expanded to other combinatorial optimization problems featuring branching structures.

Joint Truck-Drone Routing Planning.
We use the pseudocode shown in Algorithm 1 (Figure 4) to plan the routing of trucks and drones. First, we classify all customers according to their demand and recycling volume, where customer set C T is only served by trucks and customer set C U is served by both trucks and drones. Ten, the customers in set C U are removed, served by drones, and added to drone trip p1. Next, we put the former node of the removed node into trip p1 as its launch node. Ten, we determine the type of the subsequent node of the removed node. If it belongs to set C T , it is used as the retrieval node of drone trip p1; if it belongs to set C U , it is added to trip p1, and we continue to judge the next node type of the subsequent node until customer or depot 0 in set C T is found. Notably, when constructing a drone path, it must meet the drone capacity and endurance. If the constraint is not met, the path is discarded. So far, a complete truck-drone routing solution has been generated.

Infeasible Solution Adjustment
Strategy. Te aforementioned procedures enable retrieval of the routes taken by both the truck and the drone; however, the no-fy zone constraint is not considered. Hereafter, a strategy for adjusting the infeasible solution with regard to the no-fy zones is presented herein. Te specifc strategy is as follows: if the customer who uses the drone service is in the no-fy zone, the customer will be added to the place where the path cost is the lowest in the carrier truck path of the drone; if the launch/ retrieve node is in the no-fy zones, the customers served by the drone in the drone trip are added to the truck route or other drone trips. All cases are compared, and the route with the lowest route cost increment is chosen. Figures 5(a) and 5(b) illustrate a visual representation of this procedure.

Algorithm Acceleration Strategy.
To improve the iteration speed of the algorithm, the following acceleration strategy is proposed to search for the optimal solution.
When searching for drone service customers, it is benefcial to maximize the efectiveness of drones by considering the drone loading critical value, driving capacity critical value, and whether it can continue to serve the following customer node. In Figure 6, the drone fies from node A to the nearest customer B when the capacity and range allow it. When customer B has been served, the capacity and endurance reach a critical node. How to determine whether to continue serving customer C becomes a key question. Tis study proposes conditions for selecting vehicle models, to avoid the algorithm searching for many nonfeasible solutions. Table 3 lists the selection criteria. If the conditions are met, the drone proceeds to serve customer C and then returns to node D; otherwise, it returns to node D directly.

Algorithm Flowchart.
Combining the above strategies, the two-stage heuristic algorithm process is shown in Figure 7.

Computational Experiments
In this section, we evaluate the efectiveness of the proposed two-stage heuristic algorithm and investigate the advantages of the DVRPSPDNF by conducting various experiments on numerical examples of diferent problem sizes. We frst explore the efectiveness of the proposed algorithm without considering the no-fy zone. Ten, we compare the performance of the proposed algorithm and mathematical model considering the no-fy zones. Finally, by testing problem instances of diferent scales, the advantages of the DVRPSPDNF over the truck-only service and single visit by drone are analyzed, and a sensitivity analysis is performed on two key parameters. Te algorithm is implemented via Matlab 2018b, and all evaluations are conducted on a single computer with the Windows 10 operating system and an Intel(R) Core (TM) i5-10200H CPU @ 2.40 GHz processor.
Because DVRPSPDNF is a new problem, no existing benchmarks exist. Terefore, we choose the P1 data set   In numerical experiments, the numerical part of the truck-and drone-related parameters refers to published reports and current practices [54,55]. Assuming that the maximum payload of the drone is 3 kg and the maximum fight distance is 20 km, the cost per km is approximately $0.078. Assuming that the maximum loading capacity of the truck is 100 kg, the cost per kilometer is approximately $0.78. Te parameter settings related to the algorithm are shown in Table 4.

Comparison Experiment for Determining the Efectiveness of the Algorithm without considering No-Fly Zones.
To assess the efcacy of the two-stage heuristic algorithm, we initially evaluated its performance on the P1 dataset without considering no-fy zones. We then compared the results with those obtained using the depth-frst search (DFS) [56] and maximum payload-improved simulated annealing (MP-ISA) [15] in terms of objective function values. Te summarized fndings can be found in Table 5. Specifcally, the Obj column represents the total route cost, the Time column denotes the running time, and the Gap is defned as the gap between the two-stage heuristic algorithm and DFS, calculated using formula Gap � (Obj Two−stage −Obj DFS )/Obj DFS * 100%. Moreover, the Gap 0 is defned as the gap between the two-stage heuristic algorithm and MP-ISA, computed using formula Gap 0 � (Obj Two−stage −Obj MP−ISA )/Obj MP−ISA * 100%.
Te experimental fndings demonstrate the superiority of the two-stage heuristic algorithm over DFS and MP-ISA in addressing the DVRPSPD problem. Specifcally, while the DFS algorithm ensures solution quality, its search efciency becomes limited as the problem scale expands, thereby hindering the discovery of optimal solutions. In contrast, both MP-ISA and the two-stage heuristic algorithm possess global search capabilities, enabling them to acquire reasonably good solutions within a relatively short time. However, as the problem size grows, the two-stage heuristic algorithm surpasses MP-ISA in terms of solution time and quality. On average, the two-stage heuristic algorithm achieves approximately a 4% reduction in the objective function value compared to MP-ISA. Hence, the two-stage heuristic algorithm has good adaptability and search ability when solving the vehicle routing problem with drones considering pickup and delivery, especially for situations dealing with large-scale problems.

Comparison Experiment for Determining the Efectiveness of the Algorithm considering No-Fly Zones.
To verify the efectiveness of the two-stage heuristic algorithm for the nofy zone problem, we conducted a performance comparison experiment with a mathematical model. Considering the DVRPSPDNF problem as an NP-hard problem, the model can only solve small-scale problems within a reasonable time frame. Terefore, we conducted detailed statistical analysis experiments for cases with sizes of fve customers. In this example, we regarded all possible node combinations as nofy zones, and the proportion of no-fy zones did not exceed 60%. We also conducted comparative experiments in medium and large-scale calculation examples; however, due to space limitations, we changed the selection method of the no-fy zone to randomly select several nodes to be placed in the no-fy zone. In the experiments, we used the Gurobi solver and the two-stage heuristic algorithm for testing, where the time limit of the Gurobi solver was set to 2 h. Te results are summarized in Table 6, where no-fy denotes customers in no-fy zones. For the solver solution, Table 6 records the upper and lower bounds (Obj_UB, Obj_LB) of the optimal solution value obtained by Gurobi, as well as its running time (Time). For the solution of the two-stage heuristic algorithm, the table records the optimal solution value (Obj) and computation time (Time) in 20 runs and defnes the gap between them as Gap 1 � (Obj two−stage − Obj solver )/Obj solver * 100%.
Te experimental results demonstrate that the two-stage heuristic algorithm achieves optimality for all the cases involving fve customers, albeit with longer computation times when compared to Gurobi. As the number of customers increases to 10, Gurobi still manages to obtain optimal solutions within the prescribed timeframe; however, the computational time escalates considerably. While the two-stage heuristic algorithm boasts relatively shorter computation times, it does sufer from a loss of computational accuracy ranging from 6% to 7%. In the cases involving 16-20 customers, both Gurobi and the two-stage Table 3: Selection criteria.  heuristic algorithm produce approximate solutions. However, compared with the two-stage heuristic algorithm, the Gurobi solver obtains a smaller objective function value, albeit at the cost of increased solving duration. For cases involving more than 25 customers, Gurobi fails to acquire either optimal or approximate solutions, whereas the twostage heuristic algorithm can promptly provide an approximate solution. Tis discrepancy arises because even a slight increment in the instance size can cause Gurobi to consume more time and memory resources. It can be seen that the solution of the solver is still unacceptably and computationally expensive; therefore, it is necessary to provide efcient heuristics for large problems. Furthermore, we found that no-fy zones have a certain impact on the total route cost by setting diferent numbers and locations of nofy zones. Because no-fy zones cannot be traversed, if there are too many no-fy zones or the location layout is unreasonable, detours may be required to avoid these areas, thereby increasing the route cost. Terefore, when optimizing the route scheme in actual delivery, the location and number of no-fy zones must be considered to reduce the route cost in a targeted manner.

Comparative Analysis of Diferent Service Modes.
To thoroughly investigate the advantages of the MDVRPSPDNF, it is imperative to conduct numerical experiments encompassing diverse service modes. Specifcally, when solely considering the utilization of trucks to service customers, the scenario is referred to as the T mode. On the contrary, if the trucks operate in conjunction with drones and a single drone trip caters to only one customer, the scenario is referred to as the SV mode. Finally, when trucks and drones operate in collaboration and a single drone trip visits multiple customers, the scenario is called the MV mode.
In this section, we use three service modes to solve 10 calculation examples under dataset A and 10 calculation examples under dataset P1. In each experiment, several nodes are randomly selected to be placed in the no-fy zone, and it is assumed that the paths between other nodes do not pass through this area. Te experimental results are shown in Table 7. In Table 7, Gap 2 represents the gap between MV and T, and the calculation formula is Gap 2 � (Obj MV − Obj T )/Obj T * 100%; Gap 3 represents the gap between SV and T, and the calculation formula is Gap 3 � (Obj SV − Obj T )/Obj T * 100%; Gap 4 represents the gap between MV and SV, and the calculation formula is Gap 4 � (Obj MV − Obj SV )/Obj SV * 100%. To observe the solution in each mode more clearly, we show the path graph obtained by the two-stage heuristic algorithm in Figures 8(a)-8(c), taking the A-n32-k5 as an example.
Upon analysis of the experimental data, results indicate that in most instances, the MV service mode's route cost is approximately 13.27% lower than that of the T service mode, while the SV service mode's route cost is about 0.52% lower than that of the T service mode. Tis highlights that in comparison to the pure truck service model, employing a joint service of trucks and drones greatly reduces route costs. Tis critical conclusion must be considered when evaluating last-mile services with delivery requirements.
Furthermore, the MV service mode's route cost is about 12.83% lower than that of the SV service mode. Tis can be attributed to the fact that drones can service more customers in the MV mode, thereby decreasing the routing costs and the distance covered by trucks. Tis fnding further corroborates the efcacy of the truck-mounted drone collaborative service in curtailing route costs within a brief period and the drones' ability to serve multiple customers in one trip.
Additionally, utilizing a combination of trucks and drones may not always decrease the delivery distance route. Tis is due to the infuence of the customer base and depot location on the outcome. Specifcally, if the depot is situated near the periphery of the customer region, the transport distance may not be reduced; however, some degree of the routing cost can still be saved. Conversely, if the depot is close to the center of all customers, the transportation distance and routing cost can be moderately reduced.
In summary, regardless of whether the transport distance is diminished, adopting a fexible distribution approach can efectively lower route costs, which is benefcial for enhancing customer satisfaction and aligning with current energy-saving and emission-reducing ecological principles. Simultaneously, this provides a potential solution for reducing labor costs in the future logistics sector.

Sensitivity Analysis.
To further study the infuence of key parameters on the experimental results, we conducted a series of tests on A-n32-k5, A-n34-k5, and A-n36-k5 examples, where we changed the drone's farthest fying distance and maximum loading capacity. In these three calculation examples, 10 levels were set for each parameter, and a total of 60 scenarios were generated. In comparison, we only change one parameter in each scenario and record the test results. First, we set the range of the farthest fight distance of the drone to start from 10 km, increase it by 10 km each time up until 10 gradients at the end of 100 km, and obtain the corresponding objective function value. Te line graph is shown in Figure 9. It is not the case for all parameters that bigger values are better; for example, when the maximum loading capacity is large, it may lead to a low actual load rate. Hence, we set the maximum loading capacity of the drone to 3 kg (the drone serves one customer at a time), 12 kg (drones serve multiple customers at a time),  Te MP-ISA * algorithm is written according to the algorithm architecture in the literature [15] and is coded by using Matlab 2018b. Its objective function is to calculate the total transportation cost of trucks and drones. and 30-100 kg (starting from 30 kg and increasing by 10 kg each time up until 100 kg, for a total of eight gradients, in which the drone can provide services for any customer except for those in the no-fy zones). We obtain the corresponding objective function value, which is depicted by a line graph in Figure 10. Figure 9 indicates that by holding other factors constant, raising the farthest fying distance of the drone will decrease the objective function value to some extent. Tis is due to the increased number of drone service nodes, which further reduces the overall route cost. Nevertheless, there is a lower limit on the route length and cost reduction, which is infuenced by the restrictions of the drone's maximum loading capacity and no-fy zones.
Based on Figure 10, we can draw a conclusion that given other constant variables, increasing the drone's load capacity leads to a reduction in the objective function value to a certain extent, as the increased capacity allows for more  service nodes and further reduces the total route cost. However, beyond a certain critical value for the drone capacity and farthest fying distance, the objective function value will not be further reduced. Tis is because the optimization of one parameter is restricted by other parameters. For example, a maximum loading capacity of 100 kg would enable the drone to serve any customer except those in no-fy zones. Although this would increase the   number of customers served, the drone's fight capability and no-fy zones would constrain it. Tus, companies developing logistics drones must consider the fight distance and load capacity, rather than solely pursuing unilateral optimization of either the fight distance or load capacity.

Conclusions
Tis study examines an integrated system that combines trucks and drones for simultaneous pickup and delivery operations while accounting for no-fy zones. Tis represents an important expansion of the DVRPSPD as the system simultaneously considers constraints such as multivisit, farthest fying distance, maximum loading capacity, no-fy zones, and truck-drone collaboration. We formulate a MILP model with the objective of minimizing the total route cost for both trucks and drones. To address large-scale problems, we develop a two-stage heuristic algorithm comprising two main components: a SA algorithm to tackle the VRPSPD problem and a heuristic approach to handle the DVRPSPD problem with consideration for no-fy zones. In the latter part, we introduce a multidimensional array to represent the collaborative delivery path involving multiple trucks and drones. Additionally, we propose a scheme for adjusting infeasible solutions and an algorithm acceleration strategy for handling no-fy zones. Furthermore, we conduct a series of computational experiments to validate the efcacy of our proposed model and algorithm. Based on the experimental outcomes, as compared to pure truck-based services and the single-visit service mode, the multivisit service mode achieves an average route cost reduction of approximately 13.27% and 12.83%, respectively. While higher payload capacity and increased fight distances confer substantial cost advantages, it is important for enterprises to not solely pursue one-sided optimization. Instead, a balanced consideration of fight distance and payload should be employed to achieve optimal results. In conclusion, this study recommends that logistics enterprises should not only focus on cost reduction and efciency enhancement through route optimization but also consider the feasibility of fexible and efective distribution vehicles and resource allocation.
To enhance the quality of pickup and delivery services, future research should investigate the incorporation of customer service time constraints. Moreover, deploying multiple drones to provide more efcient service is a promising avenue for expansion. Finally, it is imperative to address the MDVRPSPDNF challenges in low-carbon logistics.

Data Availability
Te data in this article can be obtained from VRP Web (dorronsoro.es).

Conflicts of Interest
Te authors declare that there are no conficts of interest.