Unmanned aerial vehicle pick-and-place path planning method with time window and stop capacity constraints

Abstract

The application discloses a UAV pick-and-deliver path planning method with time window and stopover capacity limit, and belongs to the technical field of UAV logistics path planning. The method constructs a mathematical model containing multiple constraints with the goal of minimizing the total flight distance and delay penalty; under the condition of relaxing the stopover capacity constraint of the station, an adaptive large neighborhood search algorithm is used to generate a basic path set of multiple constraints, and then the final extended path pool is obtained by expansion; then the mathematical model is decoupled into a main problem and a sub-problem based on a logic-based Benders decomposition framework, the stopover capacity constraint of the station is re-included in the solution, and the distribution scheme satisfying all the constraint conditions is solved through the iterative interaction of the main problem and the sub-problem. The application realizes two-stage collaborative optimization, improves the calculation efficiency while ensuring the solution quality, and can effectively support the path planning and scheduling decision-making in the urban UAV distribution scenario.

Description

Technical Field

[0001] This application belongs to the field of drone logistics path planning technology, specifically involving a drone pickup and delivery path planning method with time window and docking capacity limitations. Background Technology

[0002] In recent years, with the deepening of the national low-altitude economic strategy and the continuous optimization of related airspace management policies, logistics drones, as an important component of the smart logistics system, are undergoing profound changes in their application scenarios. Drone delivery models have gradually expanded from the relatively closed, single-route deliveries within industrial parks to more open and networked urban scenarios. Under this trend, drone delivery paths have also evolved from simple point-to-point travel between fixed locations to a networked path planning problem that requires comprehensive consideration of various hard infrastructure constraints and coordinated scheduling within a complex network of sites, placing higher demands on the problem-solving capabilities of optimization algorithms. In urban last-mile logistics scenarios, drone path planning faces multiple complex constraints, primarily including: urban drone stations are limited by physical space and safety control requirements, with only a fixed number of docking spots, creating a hard capacity constraint on drone landing, loading / unloading, and battery replacement operations; the delivery network contains various types of spatial nodes, including multi-functional station nodes that can simultaneously serve as pickup points, delivery points, origin / termination stations, and battery swapping stations, as well as independent customer nodes that exist independently of stations and only function as pickup or delivery points. The spatial distribution and functional combinations of these heterogeneous nodes significantly increase the complexity of path planning; furthermore, given the high time-sensitivity requirements of urban drone delivery, the necessary pickup and delivery time windows further increase the difficulty of drone scheduling and planning under the hard capacity constraint of stations. Simultaneously, when battery levels meet range requirements, drones can continuously visit multiple nodes without returning to the station. While this operational flexibility helps improve delivery efficiency, it also leads to a more complex routing structure and significantly increases computational difficulty during path identification.

[0003] For the UAV path planning problem, existing technical solutions mainly fall into two categories: one is to use commercial exact solvers such as CPLEX and Gurobi to solve small-scale examples and obtain theoretical solutions by establishing a mixed integer programming model; the other is to design genetic algorithms, ant colony algorithms, tabu search algorithms and various hybrid heuristic algorithms, and design heuristic rules for specific scene characteristics to quickly obtain feasible path solutions.

[0004] However, existing technical solutions still have the following shortcomings: First, although commercially available exact solvers can obtain theoretical solutions, their computational complexity increases exponentially with the problem size, resulting in insufficient solution speed and making it difficult to obtain feasible solutions for large-scale city-level delivery tasks within a reasonable timeframe. Second, existing heuristic algorithms are mostly designed for specific simplified scenarios and cannot simultaneously achieve coupled optimization of UAV scheduling and path planning while considering station docking capacity constraints. They also lack unified modeling capabilities for complex factors such as multi-functional nodes, continuous service, and battery swapping decision coupling. Third, heuristic rules are strongly coupled with problem characteristics, resulting in poor scenario adaptability, high algorithm migration costs, and difficulty in adapting to the differentiated operational needs of different enterprises and cities. Fourth, it is difficult to balance solution quality and computational efficiency. Traditional metaheuristic algorithms experience a sharp drop in efficiency under large-scale instances, and the coupled optimization of continuous service and battery swapping decisions is insufficient, failing to fully unleash the potential for extended range. Fifth, existing methods lack effective constraint modeling and verification mechanisms when dealing with station docking capacity constraints, making it difficult to simultaneously handle capacity feasibility verification and conflict correction issues during path optimization. Summary of the Invention

[0005] To address the aforementioned problems in the existing technology, this application provides a UAV delivery path planning method with time window and docking capacity limitations. The technical problem to be solved by this application is achieved through the following technical solution: A method for planning unmanned aerial vehicle (UAV) pickup and delivery routes with time windows and docking capacity limitations includes: S1: Construct a mathematical model for urban drone delivery route planning. The optimization objective of the mathematical model is to minimize the sum of the total drone flight distance and the penalty cost for delayed delivery. The constraints of the mathematical model include drone route structure constraints, drone endurance feasibility constraints, drone time constraints, and station docking capacity constraints. S2: An adaptive large neighborhood search algorithm is used to generate a basic path set that satisfies the UAV path structure constraints, UAV endurance feasibility constraints, and UAV time constraints. The basic path set is then expanded to obtain the final expanded path pool. S3: Based on the final extended path pool, the mathematical model is decoupled from the main problem and sub-problems using a logic-based Benders decomposition architecture, and the station docking capacity constraint is reintroduced into the solution. The delivery scheme that satisfies all constraints is solved through iterative interaction between the main problem and the sub-problems.

[0006] This application has the following beneficial effects: (1) This application constructs a mathematical model that includes constraints on the path structure of delivery drones, constraints on the feasibility of drone endurance, constraints on drone time and constraints on station docking capacity. The multi-functional integrated drone station is split and represented independently according to its function, and a heterogeneous node network model is constructed. This effectively solves the complex combination optimization problem of deep coupling between station scheduling requirements and path planning caused by docking capacity constraints. At the same time, this application sets hard time windows for pickup nodes and soft time windows for delivery nodes. While meeting the strict requirements for pickup time, it also takes into account the economic trade-off of delivery delay costs, which significantly improves the algorithm's ability to cope with real-world operational scenarios.

[0007] (2) This application adopts adaptive large neighborhood search to generate a basic path set that satisfies the endurance constraint and time window constraint and expands the path. By searching under the condition of relaxing the station docking capacity constraint, it overcomes the problem of limited search space caused by the tight coupling between capacity constraint and path structure in traditional methods. It can fully explore the diverse potential of path structure without being constrained by docking capacity. At the same time, through the combination design of three destruction operators and three repair operators and the adaptive weight adjustment mechanism, the algorithm's generalization ability to different structural examples is guaranteed. Combined with hash table deduplication exploration and simulated annealing acceptance criteria, the algorithm's ability to escape local conditions is further enhanced, ensuring the full generation of high-quality paths.

[0008] (3) This application further supplements the basic path pool generated by ALNS with paths that have the potential for improvement by using short path expansion strategy and negative test number path expansion strategy. The short path expansion strategy rapidly increases path diversity by combining short paths of station-task-station, while the negative test number path expansion strategy selects paths with negative test numbers based on duality theory. The two expansion strategies complement each other and progress step by step, providing a high-quality path resource reserve for path combination optimization.

[0009] (4) This application adopts a logic-based Benders decomposition architecture to decouple the mathematical model from the main problem and sub-problems. The first sub-problem (which only includes the station docking capacity constraint) is used to quickly determine the capacity feasibility of the delivery scheme. The second sub-problem (which aims to minimize the delayed delivery time and adjusts the docking time of each UAV at each node while keeping the node access sequence unchanged) is used to accurately correct the capacity conflict. The feasible region is gradually shrunk by combining the iterative feedback mechanism of cut and feasible cut, which takes into account both the solution efficiency and the feasibility of the scheme, and realizes the collaborative optimization of path selection and station scheduling.

[0010] (5) This application breaks through the bottleneck of exponential growth in problem size of traditional commercial solvers by using the branch and bound algorithm to solve the main problem integer programming model and the iterative interaction mechanism between the main problem and subproblems. It can solve the large-scale dense distribution task at the city level in a short time, and provides technical support for the large-scale commercial operation of low-altitude logistics network.

[0011] The present application will be further described in detail below with reference to the accompanying drawings and embodiments. Attached Figure Description

[0012] Figure 1 This is a flowchart illustrating a method for unmanned aerial vehicle (UAV) pickup and delivery path planning with time windows and docking capacity limitations provided in this application. Figure 2 This is a schematic diagram of the overall algorithm framework provided in this application; Figure 3 This is a schematic diagram of the ALNS algorithm corruption-repair process provided in this application. Detailed Implementation

[0013] The present application will be described in further detail below with reference to specific embodiments, but the implementation of the present application is not limited thereto.

[0014] Example 1 like Figure 1 As shown, this embodiment provides a method for planning unmanned aerial vehicle (UAV) pickup and delivery routes with time windows and docking capacity limitations, including steps S1 to S3, as follows: S1: Construct a mathematical model for urban drone delivery route planning.

[0015] In urban last-mile logistics scenarios, drone path planning faces multiple complex constraints. Urban drone stations, limited by physical space and safety control requirements, are configured with only a fixed number of docking spots, creating a hard capacity constraint on drone landing, loading and unloading, and battery replacement operations. The delivery network contains various types of spatial nodes, including multi-functional station nodes that can simultaneously function as pickup points, delivery points, origin / termination stations, and battery swapping stations, as well as independent customer nodes that exist solely as pickup or delivery points. The spatial distribution and functional combinations of these heterogeneous nodes significantly increase the complexity of path planning. When battery levels meet range requirements, drones can continuously visit multiple nodes without returning to the station. While this operational flexibility helps improve delivery efficiency, it also leads to a more complex routing structure and significantly increases computational difficulty during path identification.

[0016] Based on the aforementioned problem scenario, the mathematical model constructed in this step aims to minimize the sum of the total drone flight distance and the penalty cost for delayed delivery. The constraints of the mathematical model include: drone path structure constraints (drones originate and terminate at stations, multiple drone paths must cover all delivery tasks, pick-up and delivery nodes for the same task must be visited consecutively by the same drone, and the number of drones at each station remains balanced after delivery); drone endurance feasibility constraints (the flight time between any two battery swaps does not exceed the maximum supported flight time of the battery); drone time constraints (the docking time at each node meets a preset time window range); and station docking capacity constraints (the number of drones docked at the same station at the same time does not exceed the station's docking capacity). The drone time constraints include soft and hard time windows.

[0017] S2: The ALNS (Adaptive Large Neighborhood Search) algorithm is used to generate a basic path set that satisfies the UAV path structure constraints, UAV endurance feasibility constraints, and UAV time constraints. The basic path set is then expanded to obtain the final expanded path pool. The core idea of ​​the ALNS algorithm is to iteratively generate diverse feasible paths through destruction and repair operations under relaxed station docking capacity constraints. Specifically, the algorithm first constructs an initial solution using a greedy strategy, then continuously improves the path quality through the synergistic effect of various destruction and repair operators. Simulated annealing acceptance criteria are used to control solution updates, and an adaptive mechanism dynamically adjusts the weights of each operator. The construction of the path pool is integrated throughout the entire ALNS process: all paths contained in solutions that pass the simulated annealing acceptance criterion are immediately added to the pool; after iteration terminates, unaccepted solutions are sorted in ascending order of objective function value, and paths from the top preset proportion of solutions are added to the path pool to achieve full exploitation and storage of path resources.

[0018] Building upon the basic path pool generated by ALNS, this method further expands the path pool using short-path expansion and negative test number path expansion strategies. The method defines a path containing only one or two pairs of delivery task pickup and delivery points as a short path. The short-path expansion strategy generates short paths by traversing all combinations of station nodes and delivery task pickup and delivery node pairs, connecting the origin station, delivery task pickup node, delivery task delivery node, and destination station sequentially. New paths that satisfy the endurance constraint and are not yet included in the path pool are added to the path pool. The negative test number path expansion strategy, based on duality theory, searches for paths with a negative test number (reduced cost) and adds them to the path pool. The test number reflects the potential improvement of the candidate path to the current objective function. For minimization problems, a negative test number indicates that introducing the path may further reduce the objective function value, thus the path has the potential to improve the current delivery plan.

[0019] S3: Based on the final extended path pool, the LBBD (Logic-based Benders Decomposition) architecture is used to decouple the mathematical model from the main problem and sub-problems, and the station docking capacity constraint is reintroduced into the solution. The delivery scheme that satisfies all constraints is solved through iterative interaction between the main problem and sub-problems.

[0020] The above steps S2 to S3 constitute the overall framework of the algorithm in this application, as follows: Figure 2 As shown in the figure, this is a schematic diagram of the overall algorithm of this application, illustrating the complete solution process from path pool generation to delivery plan output. After the algorithm starts, the ALNS algorithm is first executed to generate a basic path pool, and then the path pool is expanded through short path expansion and negative test number path expansion strategies to obtain the final expanded path pool. Entering the main problem-solving phase, the branch-and-bound algorithm is used to solve the integer programming model of the main problem. If the solution of the current node does not meet the integer requirement, a branch operation is performed on that node, generating child nodes, and the search for child nodes continues. If an integer solution is obtained, the subproblems (first subproblem SP1 and second subproblem SP2) are initiated for feasibility testing and adjustment: Subproblem SP1 is used to determine whether the current delivery plan meets the station stopping capacity constraint. If it does, a feasible solution is obtained; if not, stopping time is further adjusted through SP2. If SP2 can find an adjustment plan that meets all constraints, and the objective function value remains unchanged after adjustment, a feasible solution is directly obtained. If the adjustment generates additional delayed delivery costs, a cut is generated and returned to the main problem. If no feasible solution is found, a feasible cut is generated and the branch is pruned. If there are still remaining nodes to be searched in the main problem, branch solving continues; otherwise, the algorithm ends and outputs the delivery plan with the smallest objective function value in the currently explored feasible region. The entire process gradually shrinks the feasible region through iterative interaction between the main problem and subproblems and iterative feedback from the cut plane, achieving efficient and accurate solution.

[0021] Example 2 This embodiment further defines the process in S2 of generating the basic path set that satisfies the UAV endurance feasibility constraint and the UAV time constraint using the Adaptive Large Neighborhood Search (ALNS) algorithm, specifically including: S21a: Initial feasible solutions are constructed using a greedy algorithm to generate an initial path set.

[0022] Let the number of available drones be... ,based on Launch drones to generate An initial path is generated, consisting of only one node: the starting point node, which is the initial deployment site node for each drone. To quickly generate the initial solution, only a relatively simple and singular path form is considered when generating the initial path: the drone starts from its initial deployment site, visits the pick-up and delivery nodes to make deliveries, and after completing each pair of pick-up and delivery node visits, the drone immediately goes to the site closest to the delivery node to swap batteries, and then continues to perform other order delivery tasks, or docks at the site to end its journey.

[0023] To generate this type of path, a service sequence is defined. For delivery tasks The path segment consists of three nodes arranged sequentially: the pickup node, the delivery node, and the station node closest to the delivery node. Since the current states (current node positions and arrival times) of different drones vary during path generation, to facilitate the calculation and recording of path costs, [the following is omitted as it is not directly related to the path segment description]. Different Differentiate the same delivery task The corresponding different service sequences.

[0024] The iterative process is as follows: First, calculate all feasible service sequences. The earliest completion time. In each iteration, among all service sequences that satisfy the battery capacity constraint... In this process, the service sequence with the shortest completion time is selected. and add it to the corresponding path At the end. Then, recalculate the path. The maximum completion time for the remaining unfinished service sequence is determined. This process is repeated until all tasks are completed. The resulting set of paths constitutes the initial path set, which contains multiple paths corresponding to multiple drones.

[0025] Step S22a: Take the initial path set as the current path set, remove some delivery tasks from all paths contained in the current path set using the destruction operator, and reinsert the removed delivery tasks into the paths of the current path set using the repair operator to generate a new path set.

[0026] At the destruction rate In the case of, there are The task is removed and reinserted in each iteration. In each iteration, a destructive operator is first randomly selected using a roulette wheel method, and then selected according to the rules of the destructive operator. The task is removed from the path, and then a repair operator is randomly selected using a roulette wheel method. The previously removed task is then reinserted into the path according to the operator's rules. The destruction-repair process is as follows: Figure 3 As shown.

[0027] Figure 3 This diagram illustrates the destruction-repair process of the Adaptive Large Neighborhood Search (ALNS) algorithm in this application, demonstrating the complete process of generating a new solution from the current solution through destruction and repair operations in one iteration. , This indicates the pickup nodes for different delivery tasks. , This indicates the delivery node for the corresponding delivery task. The dashed line with an arrow indicates the pairing relationship between nodes or the direction of the insertion operation. Figure 3 (a) shows the current path, displaying a complete drone flight path including station nodes. Pickup nodes for two delivery tasks , and their corresponding delivery nodes , (a) Nodes are arranged in the order of access; (b) shows the path after destruction, displaying the destruction rate. Below, the disruption operator randomly selects a portion of the delivery tasks (such as delivery tasks). , (a) The state after removal from the path: The removed node is removed from the current path and awaits subsequent repair; (b) The repaired path: This shows the result after the repair operator re-inserts the removed delivery task into the path according to the insertion rules. The removed node is rearranged to a suitable position on the path, forming a complete new path. In each iteration, ALNS first randomly selects a destruction operator to perform a destruction operation on the path in the current solution using the roulette wheel method, and then randomly selects a repair operator to perform a repair operation on the destroyed path using the roulette wheel method, generating a new set of paths as the basis for the next iteration, thus realizing continuous exploration and optimization of the solution space.

[0028] Step S23a: Update the current path set according to the objective function value of the new path set, and update the weights of the destruction operator and the repair operator according to the search effect of the new path set.

[0029] If the objective function of the new solution is better than the current solution, then the current solution is updated to the new solution. Otherwise, if the new solution is better than the current solution or worse than the current solution but passes the simulated annealing acceptance criterion, then the current solution is updated to the new solution. In the next iteration, the ALNS algorithm performs destruction and repair based on the new current solution.

[0030] The "adaptive" in adaptive large neighborhood search refers to the fact that the weights of each repair and destruction operator are adaptively adjusted during the iteration process based on the search results of the operators. The algorithm includes three destruction operators and three repair operators. and These represent the initial selection probabilities for each destruction and repair operator. Each destruction and repair step yields a feasible solution. A new solution was obtained. Update the scores of the corresponding destruction and repair operators according to the rules in Table 1: Table 1. Parameters related to ALNS operator weight update

[0031] Table 1 shows that if For the current solution, the operator score increases. ;like The objective function is inferior to the current solution but superior to it. Operator score increases ;like The objective function is not as good as However, the operator score increased after passing the acceptance criterion. .in, .

[0032] ALNS per iteration This is followed by a weight adjustment phase, during which the scores of each operator are reset. Let the operator... No. The weight before the update was: The updated weights are The operator weights are updated using the following formula: ; in, Operator In the current period The total score in the game, Indicates the operator in the current period The number of times it is selected This is a custom weight adjustment coefficient used to adjust the proportion of the two parts in the weight calculation formula. And the operator... The probability of being selected is updated to .

[0033] S24a: Repeat S22a to S23a until the preset iteration termination condition is reached. Add the current path set, each path in the current path set, and some paths that have not passed the test to the basic path set to obtain the basic path set that satisfies the UAV path structure constraint, UAV endurance feasibility constraint, and UAV time constraint.

[0034] The destruction-repair and weight update steps are repeated until the total number of iterations of the ALNS algorithm or the number of iterations for no improved solution reaches the upper limit, at which point the algorithm stops iterating. Path pool construction is carried out throughout the entire ALNS process: all paths contained in solutions that pass the simulated annealing acceptance criterion are immediately added to the pool; after iteration terminates, unaccepted solutions are sorted in ascending order of objective function value, and paths from the top-preset proportion of solutions are added to the path pool. The top-preset proportion is a fixed numerical ratio pre-defined in the algorithm, ranging from 0 to 1, and in most practical applications, it can be set to 5%-10%.

[0035] Example 3 This embodiment further defines the specific types of the destruction operator and the repair operator, as follows: Destruction operators include random removal, Shaw Removal (similarity-based removal), and Worst Removal. Repair operators include greedy insertion, regret insertion, and regret-m insertion.

[0036] The random removal operator is used to randomly remove a specified number of pickup and delivery nodes corresponding to delivery tasks from all paths contained in the current path set.

[0037] The Shaw Removal (similarity-based removal) operator is used to remove delivery tasks with high similarity based on their similarity coefficient. Let... This is the set of all currently available paths. Define one of the paths. For two delivery tasks and Similarity coefficient between The smaller the similarity coefficient, the higher the similarity between the two delivery tasks. The calculation method is as follows: ; The above formula measures the similarity of delivery tasks by the difference between their coordinate locations and time windows, using distance. Dimension weight parameters and time window Dimension weight parameters Weighting is applied. In the formula... and These represent delivery tasks. The corresponding pickup and delivery nodes, Indicates delivery task and The distance between pickup nodes, Indicates delivery task and The distance between delivery nodes. and These represent delivery tasks. The earliest pickup time window at the pickup node and the latest delivery time window at the delivery node. and It is a pre-defined value, preset by the user according to the problem scenario.

[0038] The Worst Removal operator removes the delivery task with the highest insertion cost in the current path. It calculates the insertion cost of each delivery task in the current path. , Representing a path The total cost in the current delivery solution, Representing a path Remove delivery task The objective function value is calculated as the sum of the total distance traveled along the path and the late delivery penalties for delivery tasks within that path. This disruptive operator prioritizes removing delivery tasks with higher costs from the path.

[0039] The repair rules for all three repair operators need to refer to the insertion cost of the delivery task, that is, the difference between the total cost of the path after inserting the delivery task and the path cost before insertion. Let... For delivery tasks In the path Insertion cost in Indicates in path The Nodes and the Nodes Insert delivery tasks in between Pickup nodes and the same delivery task Delivery nodes Inserted at the pickup node The next position after the insertion point forms a new node sequence. The total cost of the new path is the increase in the total cost of the original path after the insertion operation.

[0040] The greedy insertion operator is used in each iteration to calculate the feasible insertion positions of each delivery task to be inserted in all available paths and calculate the corresponding insertion cost, and select the delivery task with the minimum insertion cost and its corresponding optimal insertion position for insertion.

[0041] Each delivery task to be inserted In the path The optimal insertion cost is Delivery task The globally optimal insertion cost is defined as That is, the minimum optimal insertion cost of the task across all available paths: In each iteration of the greedy insertion operator, select from all delivery tasks to be inserted. Minimum delivery task Insert it into the corresponding optimal insertion position (i.e. Get the minimum value The path where the time is located (and the corresponding insertion positions); then, recalculate the optimal insertion cost for other delivery tasks that have not yet been inserted on the affected path. Repeat the above steps until all delivery tasks have been re-inserted.

[0042] The regret insertion operator is used in each iteration to calculate the insertion cost difference between the optimal and second-best insertion positions for each delivery task to be inserted. The delivery task with the largest insertion cost difference is then inserted into its optimal insertion position. Optimal insertion cost among all paths Sort in ascending order, and record... For the sorted number The path number corresponding to each optimal insertion cost, where Indicates the optimal insertion position. This indicates the second-best insertion position, and so on. . Indicates delivery task In the path The optimal insertion cost is the insertion path for the task (paired pickup and delivery nodes). After a certain position, the increase in the total cost of the new path relative to the original path; Indicates delivery task After sorting the optimal insertion costs on all feasible paths in ascending order, the th The path index (or location identifier) ​​corresponding to the optimal insertion cost, where This indicates the optimal insertion position (lowest cost). This indicates the suboptimal insertion position, and so on; let the difference be... In each iteration, the delivery task with the largest difference is inserted into its insertion position, and the insertion cost and difference of other delivery tasks to be inserted in the path are recalculated. The above steps are repeated until the re-insertion of all delivery tasks is completed.

[0043] The regret-m insertion operator is used in each iteration to calculate the difference between the sum of the insertion costs of the previous m optimal insertion positions and the insertion cost of the optimal insertion position for each delivery task to be inserted. The delivery task with the largest difference is then inserted into its optimal insertion position. Regret-m is a generalization of the regret operator. The difference is defined as follows: ; In the formula, Indicates the delivery task number; Indicates before consideration One optimal insertion position; Indicates delivery task In its first Optimal insertion position Insertion cost; Indicates delivery task forward The sum of insertion costs at each optimal insertion position; Indicates delivery task At the optimal insertion position The insertion cost (i.e., minimum insertion cost) on the [top / bottom].

[0044] therefore, Indicates delivery task forward The difference between the sum of costs for each optimal insertion position and the optimal insertion position itself is the value of the task. A larger value indicates a greater cost loss if the task is not inserted at the optimal position, and therefore, it should be prioritized for insertion at the optimal position. The regret-m operator selects the optimal position in each iteration. Insert the largest uninserted delivery task into its optimal insertion position. .

[0045] Similar to the regret heuristic, regret-m inserts the delivery task with the largest m-difference at its insertion position in each iteration, and recalculates the insertion cost and cost difference of other delivery tasks to be inserted in that path. This process is repeated until all delivery tasks have been re-inserted. If there are insufficient insertable paths for some delivery tasks... Article (i.e.) If the number of insertable paths is the smallest and the total difference of remaining insertable paths is the largest, then the delivery task with the fewest insertable paths and the largest total difference of remaining insertable paths will be inserted into its position first.

[0046] Example 4 This embodiment further defines the weight update process, simulated annealing acceptance criteria, and iteration termination conditions, specifically including: In S23a, when updating the weights of the destruction and repair operators, the weights of the destruction and repair operators used in the current iteration are only updated based on the search performance of the new path set if the new path set is a new path set that has not been explored in previous iterations. Specifically, in order to explore high-quality feasible solutions within a larger solution space, solutions stored in a hashtable (a data structure used to store explored solutions) are updated if and only if... When a new solution is found that has not been explored in previous iterations, its weights are updated.

[0047] In S24a, when adding a path to the base path set, the simulated annealing acceptance criterion is used to determine whether to accept the new path set to replace the current path set. The simulated annealing acceptance criterion accepts the new solution with the probability shown in the following formula. : ; In the formula, To simulate the temperature coefficient of annealing, its initial value is set to... The objective function value is required to be different from the initial solution in the first iteration. The solution is The probability is accepted. The temperature coefficient is based on a fixed cooling rate. Declined generation by generation ( That is, after each iteration, the temperature coefficient is updated to... . When understood Better than the original feasible solution Time (i.e.) hour), This means that a new solution superior to the original feasible solution must be determined by the acceptance criterion; otherwise, the solution is accepted. The probability is between 0 and 1, and The lower the objective function value, the higher the likelihood that it will pass the acceptance criterion.

[0048] If the new path set passes the simulated annealing acceptance criterion, it is used as the current path set. If the new path set fails the simulated annealing acceptance criterion, the paths in the new path set are sorted in ascending order of objective function value, and the top preset proportion of these paths are added to the base path set along with the paths in the current path set. The iteration terminates when the total number of iterations reaches a preset upper limit or the number of iterations without improved solutions reaches a preset upper limit.

[0049] Example 5 This embodiment further defines the process of expanding the basic path set in S2, specifically including: S21b: Perform a short path expansion strategy on the basic path set to obtain a path pool after short path expansion.

[0050] Path expansion strategy input: known site set and delivery task collection It iterates through all station nodes and delivery tasks, forming a sequence starting from the originating station. Delivery tasks Pickup nodes Same delivery task Delivery nodes Destination Station Shortest path formed by connecting them sequentially ,in and All are collections of sites The site nodes in the [location]. If For path pool A new path not included in the plan, and the path is feasible (meets the drone's endurance feasibility constraint, i.e., the flight time from takeoff to the first battery replacement is within the maximum flight time supported by the battery), will join in In the shortest path, s 1 indicates the starting station. s 2 represents the destination station; both belong to the station set. S .

[0051] For short routes containing two delivery tasks, to reduce unnecessary calculations caused by route infeasibility due to failure to switch batteries midway, the delivery tasks are pre-screened. The specific screening process is as follows: Iterate through all delivery task pairs that contain two different delivery tasks. ,in and These represent two different delivery task numbers. Calculate the first delivery task. Delivery nodes With the second delivery task Pickup nodes The distance between them. If this distance does not exceed the preset distance limit, then... Define them as similar task pairs and add them to the similar task pair set. The preset distance limit is a fixed value parameter set in advance to determine whether two delivery tasks constitute a similar task pair. Its value is preset by the user based on the station distribution density and the drone's endurance, and this application does not impose specific limitations on it. In practical applications, this value is usually set within the range of 800 meters to 1500 meters, which can cover most urban last-mile delivery scenarios.

[0052] Subsequently, all sites and the set of similar task pairs are traversed. The task pair is used to find a feasible short path that includes both delivery tasks, and new paths that meet the drone's endurance feasibility constraints are added to the path pool. The short path takes the form of... That is, from the starting station Depart and visit the delivery tasks in sequence. Pickup nodes , Delivery nodes Delivery tasks Pickup nodes , Delivery nodes Finally arrived at the destination station ,in and All belong to the site collection .

[0053] S22b: Apply the negative test number path expansion strategy to the path pool after short path expansion to obtain the final expanded path pool.

[0054] According to duality theory, the test number of each path reflects the potential improvement in the objective function value relative to the current dual solution. That is, for a minimization problem, paths with negative test numbers have the potential to improve the delivery solution. The linear relaxation model of the path combination set covering model is as follows, with the existing path pool forming the basis of this model: ; in, System represents path The 0-1 variable indicating whether a product is selected for the delivery plan is relaxed to a continuous variable not less than 0 in the relaxation model. The route cost consists of two parts: the drone's flight distance and the penalty for delayed delivery. and Representing paths Is it based on the site? For the starting point and path Does it serve delivery tasks? For all existing paths in the path pool, their corresponding and Given parameters, This represents the maximum number of available drones.

[0055] This integer programming model aims to minimize path costs and combines paths in the path pool into delivery schemes under the constraints that the number of drones called up does not exceed the number of available drones and the paths cover all delivery tasks. This represents the set of all paths in the path pool.

[0056] Based on the derivation of duality theory, the path can be obtained. The test numbers are shown below: ; in, In the original problem The dual variables corresponding to the capacity constraints of each starting site can be viewed as the site's dual variables. The shadow price as the starting point; In the original problem The dual variables corresponding to the service constraints of each delivery task can be regarded as completing the delivery task. The shadow price. and Representing paths Is it based on the site? For the starting point and path Does it serve delivery tasks? For all existing paths in the path pool, their corresponding and The parameters are known.

[0057] The steps for expanding the negative test result path are as follows: The first step is to construct a zero-test-number path pool. First, based on the basic path pool, a corresponding set covering problem is constructed. Then, the dual problem model of the set covering problem, dualRP, is constructed, and the dualRP is solved to obtain the values ​​of the dual variables. Next, the calculation... Calculate the test count of all paths in the path pool, and select paths with a test count of 0 to add to the path pool. .

[0058] The second step involves performing a local search on paths with zero test scores and calculating the test score, then adding paths with negative test scores to the path pool. For... In each path, the five designed local search operators are called sequentially: (1) Change the starting point. Change the path The original starting point Change to another site If the changed path meets the drone's endurance feasibility constraint (i.e., the drone's range from the new starting point after the change), If the flight time from takeoff to the first battery replacement does not exceed the maximum flight time supported by the battery, then the path cost is calculated as follows: Become After changing the starting station, the path verification count changed from 0 to: ; in, Indicates the original path The total cost consists of two parts: the drone's flight distance and penalties for delayed delivery. Indicates the path The starting station is from Replace with Then, the total cost of the new path; Indicates the original starting point station The corresponding dual variable value, Indicates the new starting point station The corresponding dual variable value.

[0059] The change in this test number reflects the change in path cost after changing the starting point. ) Excluding the changes in dual variables caused by changes in the starting station ( If this value is less than 0, the new path has a negative test number.

[0060] (2) Replace the battery swapping station or the destination station. Traverse all battery swapping stations and destination stations along the path, and try to replace them with other stations. If the path is feasible after the replacement, the path check count becomes .

[0061] (3) Delete a delivery task The corresponding pickup and delivery nodes. Traverse all pickup and delivery nodes in the path, and sequentially attempt to delete pairs of pickup and delivery nodes corresponding to the same delivery task. According to the triangle inequality, the new path after deleting nodes has a shorter flight distance between two battery swaps than the original path's flight distance, thus satisfying the power feasibility requirement. The verification number is... .

[0062] (4) Insert a delivery task that is not originally included in the path. When inserting a new location for pickup and delivery nodes, it is necessary to ensure that the existing paired pickup and delivery points are not separated. The new path check count after insertion becomes... .

[0063] (5) Replace one delivery task in the path with another delivery task that was not originally included in the path. Iterate through each delivery task included in the path. Replace its pickup and delivery node with another delivery task. If the pickup and delivery nodes are changed and the route remains feasible, the number of checks becomes [number missing]. . Indicates delivery task The shadow price.

[0064] In the local search, for each new feasible path found, its test number is calculated using the dual variable and path parameters. New paths with negative test numbers are added to the negative test number path pool. After completing a local search for all paths with zero test scores, then... Add to current path pool Clear and Repeat the negative test result path generation steps above until the maximum number of repetitions is reached. until.

[0065] Example 6 This embodiment further defines S3, specifically including: S31: The LBBD (Logic-based Benders Decomposition) architecture is adopted to decouple the mathematical model from the main problem and sub-problems. The UAV path structure constraints are retained in the main problem, while the site docking capacity constraints, UAV endurance feasibility constraints and UAV time constraints are handled in the sub-problems.

[0066] Given the significant increase in problem-solving complexity due to the coupling of various constraints such as station docking capacity, endurance, and time window constraints, this step employs the Dantzig-Wolfe Decomposition approach to decouple the original problem. Specifically, the station docking capacity constraint, which is difficult to model and solve directly, as well as the time window and endurance constraints affected by station docking scheduling, are decomposed into sub-problems, constructing a master-sub-problem collaborative solution architecture. The master problem uses a branch and bound algorithm to solve an integer programming model, with the objective function being to minimize the sum of the total drone flight distance and the late delivery penalty cost. Path combinations are selected from the path pool to form a complete delivery plan. The drone path structure constraints cover the following three aspects: task coverage constraint, ensuring that all pending delivery tasks are covered; drone quantity constraint, ensuring that the total number of drones called does not exceed the deployment scale; and station traffic balance constraint, ensuring that the number of drones docked at each station remains dynamically balanced after delivery is completed.

[0067] The sub-problem serves as a verification module, modeled using constraint programming (CP). It is invoked only when a feasible solution to the main problem is obtained, to verify the feasibility of the delivery plan's station docking capacity. Due to the specific nature of the station docking capacity constraint, it is separated into a verification sub-problem. The core mechanism of the constraint programming model is as follows: interval variables containing start and end times represent the drone's docking behavior at each node. Multiple docking interval variables are arranged sequentially to form a path sequence, and then the station docking capacity limit is expressed based on the cumulative resource constraints of these interval variables.

[0068] S32: The station stopping capacity constraint is reintroduced into the solution. Based on the final expanded path pool, the delivery scheme that satisfies all constraints is solved through iterative interaction between the main problem and sub-problems.

[0069] Given the computational overhead of accurately calculating station docking times, a two-level checking subproblems are constructed to achieve hierarchical verification. After obtaining the delivery plan from the main problem, it is first input into the first subproblem (SP1, a constraint satisfaction model without an objective function, containing only the cumulative resource constraint of station docking capacity) for verification. If the capacity constraint is satisfied, the solution is retained and returned to the main problem, driving the subsequent branch and bound process based on the objective function value; if not, the solution is input into the second subproblem (SP2) for in-depth adjustment. The second subproblem aims to achieve capacity constraint satisfaction by adjusting the drone station docking times. It uses minimizing the delayed delivery time as the objective function, adjusting the docking time of each drone at each node while keeping the node access sequence unchanged, so that the adjusted docking time period satisfies the station docking capacity constraint, the drone endurance feasibility constraint, and the drone time constraint.

[0070] Regarding the cutting plane generation mechanism: If adjusting the second subproblem leads to an increase in delivery delay time, an optimization cut is added to the main problem to provide lower bound information for the delivery delay penalty term of the solution; if the second subproblem determines that the solution cannot satisfy the capacity constraint through time adjustment, a feasibility cut is added to the main problem to permanently eliminate subsequent occurrences of the solution. The above master-subproblem interaction gradually shrinks the feasible region through iterative feedback, achieving efficient and accurate solution.

[0071] Example 7 This embodiment S32 is further defined, specifically including: S321: Using each path in the final expanded path pool as a decision variable, construct and solve the integer programming model of the main problem to obtain a delivery solution consisting of several paths. The delivery solution satisfies the UAV path structure constraints in the main problem.

[0072] Specifically, the main problem uses minimizing the sum of the total drone flight distance and the penalty cost for delayed delivery as the objective function. Each path in the path pool is used as a decision variable, which is a 0-1 variable indicating whether the path is selected into the delivery plan. The drone path structure constraints cover the following three aspects: task coverage constraint, ensuring that each delivery task is covered by at least one selected path; available drone quantity constraint, ensuring that the total number of selected paths does not exceed the number of available drones; and station traffic balance constraint, ensuring that the number of selected paths as starting points and the number of selected paths as ending points at each station are dynamically balanced.

[0073] S322: Input the delivery plan into the subproblem, and determine whether the delivery plan meets the station parking capacity constraint, the drone endurance feasibility constraint, and the drone time constraint by solving the subproblem. If it is determined that they are met, then the delivery plan meets all the constraints in S1.

[0074] The subproblems are solved using a Constraint Programming (CP) model. The CP model uses interval variables to represent the drone's docking behavior at each node. These interval variables include the docking start time and the docking end time, and multiple interval variables are arranged sequentially to form the docking sequence for each path. The constraints of the subproblems include: cumulative resource constraints such as the number of drones docking at the same site at the same time not exceeding the site's docking capacity; drone time constraints such as the docking time at each node meeting a preset time window; and drone endurance feasibility constraints such as the flight time between any two battery swaps not exceeding the battery's maximum supported flight time.

[0075] S323: If the delivery plan satisfies all the constraints in S1, then the delivery plan is returned to the main problem as a feasible solution; if the delivery plan does not satisfy all the constraints in S1, then the delivery plan is adjusted through subproblems, and an optimization cut or feasibility cut is added to the main problem based on the adjustment results. The subproblems include the first subproblem (SP1) and the second subproblem (SP2).

[0076] Specifically, if the first subproblem determines that the delivery plan satisfies the station docking capacity constraint, then since the delivery plan has already satisfied the drone path structure constraint in the main problem, it is determined that the delivery plan satisfies all constraints and is returned to the main problem as a feasible solution. If the first subproblem determines that the delivery plan does not satisfy the station docking capacity constraint, then the docking time is adjusted through the second subproblem. If the second subproblem can find an adjustment plan that satisfies all constraints, then a cut is added to the main problem based on the adjusted delayed delivery time; if the second subproblem cannot find an adjustment plan that satisfies all constraints, then a feasible cut is added to the main problem.

[0077] S324: Return to S321 until the main problem is solved, and output the delivery plan with the minimum objective function value.

[0078] The branch and bound algorithm systematically divides the feasible region and solves for linear relaxation lower bounds, gradually converging to an integer solution: It solves the linear relaxation problem at each node of the search tree; if the solution does not meet the integer requirement, it selects a fractional variable to branch, generating child nodes to continue the search; simultaneously, it updates the global upper bound using the node objective value and prunes branches with lower bounds worse than the current solution, until all nodes are traversed or pruned. After iterative convergence, the delivery solution that minimizes the objective function value of the main problem within the currently explored feasible region and satisfies all constraints as verified by subproblems is output as the delivery solution.

[0079] Example 8 This embodiment further defines S322, specifically including: S3221: Input the delivery plan into the first subproblem (SP1), and determine whether the delivery plan satisfies the station parking capacity constraint by solving the first subproblem. The first subproblem is a constraint satisfaction model that only includes the station parking capacity constraint.

[0080] Given the computational overhead of accurately calculating station stop times, this embodiment constructs a two-level checking subproblems to achieve hierarchical verification. The first subproblem (SP1) is a constraint satisfaction model without an objective function, containing only a cumulative resource constraint of station stop capacity. The input is multiple path sequences corresponding to the delivery plan, and SP1 only determines whether the plan satisfies the station capacity constraint.

[0081] S3222: If the first subproblem determines that the delivery plan satisfies the station docking capacity constraint, then considering that the delivery plan has already satisfied the UAV path structure constraint in the main problem, it is determined that the delivery plan satisfies all the constraints in S1.

[0082] S3223: If the first subproblem determines that the delivery plan does not meet the station parking capacity constraint, then it is determined that the delivery plan does not meet all the constraints in S1.

[0083] If the first subproblem determines that the delivery plan meets the station parking capacity constraint, then the solution is retained and returned to the main problem, and the subsequent branch and bound process is driven by the objective function value; if the first subproblem determines that it does not meet the constraint, then the solution is input into the second subproblem for in-depth adjustment.

[0084] Example 9 This embodiment further defines S323, specifically including: S3231: If the first subproblem (SP1) determines that the delivery plan does not meet the site docking capacity constraint, then the delivery plan is input into the second subproblem (SP2). By solving the second subproblem, it is found whether there is an adjustment plan that meets the site docking capacity constraint, the drone endurance feasibility constraint, and the drone time constraint. The second subproblem aims to minimize the delayed delivery time. While keeping the node access sequence unchanged, it adjusts the docking time of each drone at each node so that the adjusted docking time period meets the site docking capacity constraint, the drone endurance feasibility constraint, and the drone time constraint.

[0085] Adjusting station stopping times involves the interconnected effects of three types of constraints: First, hovering waiting constraints, where changes in stopping time may increase the hovering waiting time before landing, a process that continuously consumes power and jeopardizes the feasibility of the route's endurance; second, time window propagation constraints, where delayed stopping will trigger a chain reaction of delays in subsequent stops (including pick-up and drop-off stops), and delivery delays will trigger or exacerbate late delivery penalties; and third, core capacity constraints, where the adjusted stopping time must meet the station's capacity limit. Therefore, the second subproblem (SP2), in addition to the basic capacity constraints, fully incorporates time variable constraints and power endurance constraints to ensure feasibility, and uses minimizing late delivery time as the objective function (SP2 only adjusts the time parameters and does not modify the node access sequence, so it does not affect the flight distance term in the original objective function).

[0086] The Constraint Programming (CP) model for the second subproblem uses interval variables to represent the docking behavior of drones at each node. These interval variables include the start and end times of docking, and multiple interval variables are arranged sequentially to form the docking sequence for each path. The constraints of the second subproblem include: cumulative resource constraints such as the number of drones docking at the same site at the same time not exceeding the site's docking capacity; drone time constraints such as the docking time at each node meeting a preset time window; and drone endurance feasibility constraints such as the flight time between any two battery swaps not exceeding the maximum flight time supported by the battery.

[0087] S3232: If the second subproblem can find an adjustment scheme that satisfies all the constraints in S1, then obtain the delayed delivery time corresponding to the adjustment scheme, and add an optimization cut to the main problem based on the delayed delivery time.

[0088] If adjusting the second subproblem results in an increase in delivery delay time, a cut is added to the main problem to provide lower bound information for the delay penalty term of the solution. A cut is a lower bound constraint added to the main problem based on the delay time in the adjusted solution regarding the delay penalty cost of that delivery solution.

[0089] S3233: If the second subproblem cannot find an adjustment scheme that satisfies all the constraints in S1, then add a feasibility cut to the main problem.

[0090] If the second subproblem determines that the solution cannot satisfy the capacity constraint through time adjustment, a feasible cut is added to the main problem to permanently eliminate subsequent occurrences of the solution. A feasible cut is a constraint used to permanently eliminate the current delivery solution and its subsequent identical path combinations.

[0091] Example 10 This embodiment further restricts the time constraints for drones and the site traffic balance constraints, as follows: In the drone time constraints, the pickup node's time window is a hard time window, requiring the drone to arrive within the time window, which must be strictly satisfied; the delivery node's time window is a soft time window, allowing for delayed delivery but incurring a late delivery penalty cost that is included in the optimization objective. Delays will trigger the penalty cost term in the objective function.

[0092] The site traffic balance constraint is: for each site, the number of selected paths starting from that site is equal to the number of selected paths ending at that site, that is, the number of drones docked at each site after delivery is completed is kept in dynamic balance.

[0093] To verify the effectiveness of this application, the effectiveness of the algorithm in this application is verified through actual numerical examples.

[0094] The test cases were generated based on the actual geographic coordinates of Meituan drone stations in Shenzhen and areas with frequent food delivery orders in Shenzhen. The cases were named in the form of nmi, where n, m, and i represent the number of tasks, the number of stations, and the case number, respectively. The test results are shown in the table below. In Table 2, FT, TT, Obj, and runTime represent the drone path flight time (in minutes), delayed delivery time (in minutes), objective function value, and algorithm runtime (in seconds), respectively.

[0095] Table 2. Algorithm solution results for different site capacities.

[0096] As shown in Table 2, the reduction in objective function value due to decreased site capacity conflicts is very significant with the increase in site capacity. However, the benefits of increasing site capacity also exhibit diminishing returns; further increasing site capacity when it is sufficient offers little improvement to the objective function value. For task-intensive examples with relatively large task scales within a small area, this application can also complete the solution in a short time, demonstrating that it can effectively overcome the bottleneck of traditional solvers in terms of problem scale and achieve efficient solutions for large-scale, dense urban delivery tasks.

[0097] The above description is merely a preferred embodiment of this application and is not intended to limit this application. Various modifications and variations can be made to this application by those skilled in the art. Any modifications, equivalent substitutions, improvements, etc., made within the spirit and principles of this application should be included within the protection scope of this application.

[0098] It is worth noting that the terms "first" and "second" in this application are used for descriptive purposes only and should not be construed as indicating or implying relative importance or implicitly specifying the number of technical features indicated. Therefore, a feature defined as "first" or "second" may explicitly or implicitly include one or more of that feature. In the description of this application, "multiple" means two or more, unless otherwise explicitly specified.

[0099] The above description, in conjunction with specific preferred embodiments, provides a further detailed explanation of this application and should not be construed as limiting the specific implementation of this application to these descriptions. For those skilled in the art, various simple deductions or substitutions can be made without departing from the concept of this application, and all such modifications or substitutions should be considered within the scope of protection of this application.

Claims

1. A method for planning unmanned aerial vehicle (UAV) pickup and delivery routes with time windows and docking capacity limitations, characterized in that, include: S1: Construct a mathematical model for urban drone delivery route planning. The optimization objective of the mathematical model is to minimize the sum of the total drone flight distance and the penalty cost for delayed delivery. The constraints of the mathematical model include UAV path structure constraints, UAV endurance feasibility constraints, UAV time constraints, and site docking capacity constraints; the UAV time constraints include soft time windows and hard time windows. S2: An adaptive large neighborhood search algorithm is used to generate a basic path set that satisfies the UAV path structure constraints, UAV endurance feasibility constraints, and UAV time constraints. The basic path set is then expanded to obtain the final expanded path pool. S3: Based on the final extended path pool, the mathematical model is decoupled from the main problem and sub-problems using a logic-based Benders decomposition architecture, and the station docking capacity constraint is reintroduced into the solution. The delivery scheme that satisfies all constraints is solved through iterative interaction between the main problem and the sub-problems.

2. The UAV delivery path planning method with time window and docking capacity limitations according to claim 1, characterized in that, S2 uses an adaptive large neighborhood search algorithm to generate a set of basic paths that satisfy the UAV path structure constraints, UAV endurance feasibility constraints, and UAV time constraints, including: S21a: An initial feasible solution is constructed using a greedy algorithm to generate an initial path set, wherein the initial path set contains multiple paths corresponding to multiple drones; S22a: Using the initial path set as the current path set, remove some delivery tasks from all paths contained in the current path set using a destruction operator, and reinsert the removed delivery tasks into the paths of the current path set using a repair operator to generate a new path set; S23a: Update the current path set according to the objective function value of the new path set, and update the weights of the destruction operator and the repair operator according to the search effect of the new path set; S24a: Repeat S22a to S23a until the preset iteration termination condition is reached. Add the current path set, each path in the current path set, and some paths that failed the test to the basic path set to obtain the basic path set that satisfies the UAV path structure constraint, UAV endurance feasibility constraint, and UAV time constraint.

3. The UAV delivery path planning method with time window and docking capacity limitations according to claim 2, characterized in that, The destruction operators include the random removal operator, the Shaw Removal operator, and the Worst Removal operator; the repair operators include the greedy insertion operator, the regret insertion operator, and the regret-m insertion operator; The random removal operator is used to randomly remove a specified number of pickup nodes and delivery nodes corresponding to delivery tasks from all paths contained in the current path set; The Shaw Removal operator is used to remove delivery tasks with high similarity based on the similarity coefficient between them. The Worst Removal operator is used to remove the delivery task with the highest insertion cost in the current path; The greedy insertion operator is used to calculate the feasible insertion positions of each delivery task to be inserted in all available paths and calculate the corresponding insertion cost in each iteration, and select the delivery task with the minimum insertion cost and its corresponding optimal insertion position for insertion. The regret insertion operator is used to calculate the insertion cost difference between the optimal insertion position and the second-best insertion position for each delivery task to be inserted in each iteration, and insert the delivery task to be inserted with the largest insertion cost difference into its optimal insertion position. The regret-m insertion operator is used to calculate, in each iteration, the difference between the sum of the insertion costs of the previous m optimal insertion positions and the insertion cost of the optimal insertion position for each delivery task to be inserted, and insert the delivery task to be inserted with the largest difference into its optimal insertion position.

4. The UAV delivery path planning method with time window and docking capacity limitations according to claim 2, characterized in that, In S23a, when updating the weights of the destruction operator and the repair operator, the weights of the destruction operator and the repair operator used in this iteration are only updated according to the search effect of the new path set if the new path set is a new path set that has not been explored in previous iterations. In S24a, when adding a path to the base path set, a simulated annealing acceptance criterion is used to determine whether to accept the new path set to replace the current path set. If the new path set passes the simulated annealing acceptance criterion, it is used as the current path set. If the new path set fails the simulated annealing acceptance criterion, the paths in the new path set are sorted in ascending order of objective function value, and the top preset proportion of these paths are added to the base path set along with the paths in the current path set. The simulated annealing acceptance criterion is used to accept the new path set with a predetermined probability if, in each iteration, the objective function value of the new path set is inferior to that of the current path set. The probability decreases as the number of iterations increases; The iteration termination condition is when the total number of iterations reaches the preset upper limit or the number of iterations without a continuously improved solution reaches the preset upper limit.

5. The UAV delivery path planning method with time window and docking capacity limitations according to claim 1, characterized in that, In S2, the basic path set is expanded to obtain an expanded path pool, which includes: S21b: Execute a short path expansion strategy on the basic path set to obtain a path pool after short path expansion; wherein, the short path expansion strategy traverses all station nodes and delivery tasks, and takes the path formed by sequentially connecting the starting station, the pickup node of a single delivery task, the delivery node of the same delivery task, and the destination station as a short path; if the short path is a new path not included in the path pool and satisfies the UAV endurance feasibility constraint, then the short path is added to the path pool; for a short path containing two delivery tasks, first filter out the close task pairs where the distance between the pickup node and the delivery node does not exceed a preset distance limit, and then generate a short path that satisfies the UAV endurance feasibility constraint based on the close task pairs and add it to the path pool. S22b: Execute a negative test number path expansion strategy on the path pool after the short path expansion to obtain the final expanded path pool; wherein, the negative test number path expansion strategy constructs a linear relaxation model of the set covering problem based on the current path pool, solves the dual problem of the linear relaxation model to obtain the dual variable values ​​corresponding to each constraint, calculates the test number of all paths in the path pool based on the dual variable values, selects paths with a test number of zero to form a zero test number path pool, performs a local search operation on each path in the zero test number path pool to generate a new path, calculates the test number of the new path, and adds the new path with a negative test number to the path pool.

6. The UAV delivery path planning method with time window and docking capacity limitations according to claim 1, characterized in that, S3 include: S31, The mathematical model is decoupled from the main problem and the sub-problems using a logic-based Benders decomposition architecture, and the UAV path structure constraints are retained in the main problem, while the station docking capacity constraints, the UAV endurance feasibility constraints and the UAV time constraints are processed in the sub-problems. S32, the station stopping capacity constraint is reintroduced into the solution, and based on the final extended path pool, a delivery scheme that satisfies all constraints is solved through iterative interaction between the main problem and the sub-problems.

7. The UAV delivery path planning method with time window and docking capacity limitations according to claim 6, characterized in that, S32 includes: S321: Using each path in the final extended path pool as a decision variable, construct an integer programming model for the main problem and solve it to obtain a delivery scheme consisting of several paths; the delivery scheme satisfies the UAV path structure constraints in the main problem; S322: Input the delivery plan into the sub-problem, and determine whether the delivery plan satisfies the station docking capacity constraint, the drone endurance feasibility constraint, and the drone time constraint by solving the sub-problem. If it is determined that they are satisfied, then the delivery plan satisfies all the constraints in S1. S323: If the delivery plan satisfies all the constraints in S1, then the delivery plan is returned as a feasible solution to the main problem; if the delivery plan does not satisfy all the constraints in S1, then the delivery plan is adjusted through the sub-problems, and a cut or feasible cut is added to the main problem according to the adjustment result; wherein, the sub-problems include a first sub-problem and a second sub-problem; S324: Return to S321 until the main problem is solved, and output the delivery solution with the minimum objective function value.

8. The UAV delivery path planning method with time window and docking capacity limitations according to claim 7, characterized in that, S322 includes: S3221: Input the delivery plan into the first sub-problem, and determine whether the delivery plan satisfies the station parking capacity constraint by solving the first sub-problem; the first sub-problem is a constraint satisfaction model that only includes the station parking capacity constraint; S3222: If the first subproblem determines that the delivery scheme satisfies the station docking capacity constraint, then, considering that the delivery scheme has already satisfied the UAV path structure constraint in the main problem, it is determined that the delivery scheme satisfies all the constraints in S1. S3223: If the first subproblem determines that the delivery plan does not meet the station parking capacity constraint, then it is determined that the delivery plan does not meet all the constraints in S1.

9. The UAV delivery path planning method with time window and docking capacity limitations according to claim 7, characterized in that, S323 includes: S3231: If the first subproblem determines that the delivery plan does not meet the station docking capacity constraint, then the delivery plan is input into the second subproblem. By solving the second subproblem, it is found whether there is an adjustment plan that meets the station docking capacity constraint, the drone endurance feasibility constraint, and the drone time constraint. The second subproblem aims to minimize the delayed delivery time. While keeping the node access sequence unchanged, the docking time of each drone at each node is adjusted so that the adjusted docking time period meets the station docking capacity constraint, the drone endurance feasibility constraint, and the drone time constraint. S3232: If the second subproblem can find an adjustment scheme that satisfies all the constraints in S1, then obtain the delayed delivery time corresponding to the adjustment scheme, and add a cut to the main problem according to the delayed delivery time; S3233: If the second subproblem cannot find an adjustment scheme that satisfies all the constraints in S1, then add a feasible cut to the main problem.

10. The UAV delivery path planning method with time window and docking capacity limitations according to claim 1, characterized in that, In the drone time constraints, the pickup node's time window is a hard time window, requiring the drone to arrive within the time window; the delivery node's time window is a soft time window, allowing for delayed delivery but incurring a late delivery penalty cost, which is included in the optimization objective. The UAV path structure constraints include: for each station, the number of selected paths starting from that station is equal to the number of selected paths ending at that station.

Images