3.1 Constraints for UAV logistics path planning

When performing specific tasks, UAVs require three-dimensional trajectory planning to ensure that they fly safely and efficiently along a predetermined route. This planning process involves a number of constraints, such as limits on turning angles, climbing and descending angles, the minimum length of a segment, the minimum and maximum altitudes to be flown, and the maximum length of a range. These constraints vary depending on the mission requirements and the actual flight environment, and aim to ensure that the UAV accomplishes its flight mission smoothly and safely ^[16]. Let each path in finding the corresponding path node as $P_{ij} =(x_{ij} $, $y_{ij} $, $z_{ij} )$, the next path node is denoted as $P_{i,j+1} =(x_{i,j+1} $, $y_{i,j+1} $, $z_{i,j+1} )$, then the cost function about the path length is shown in Eq. (1).

In Eq. (1), the Euclidean distance between two points is $\left\| \overrightarrow{P_{ij} P_{i,j+1} }\right\| $. UAVs will encounter various types of obstacles on the way of flight, which will have an impact on the operating condition of UAVs. Therefore, the way of modeling obstacles can be used to guarantee the normal flight of the UAV. The obstacle collision is shown in Fig. 1.

Fig. 1. Schematic diagram of obstacle collision.

In Fig. 1, The green area represents the collision space. $D$ denotes the ring width of the collision area. $K$ denotes the set of all obstacles, and each obstacle is regarded as a cylinder with center projection coordinates denoted as $C_{k} $ and radius denoted as $R_{k} $. For a given path $\left\| \overrightarrow{P_{ij} P_{i,j+1} }\right\| $, there exists an associated cost proportional to its distance from $d_{k} $ to$C_{k} $. The simulation of the obstacle modeling is shown in Fig. 2.

Fig. 2. Simulation diagram of three-dimensional obstacles.

In Fig. 2, the obstacle is viewed as a bar, and the starting point is set to be $P_{s} (x_{1} $, $y_{1} $, $z_{1} )$ and the end point to be $P_{e} (x_{m} $, $y_{m} $, $z_{m} )$. The cost function of the threat generated by the obstacle is shown in Eq. (2).

In Eq. (2), $S$ denotes the ring width of the hazardous area. If the feasibility and safety constraints of obstacles in the environment are considered, there exists a total cost function as shown in Eq. (3).

In Eq. (3), $F(X_{i} )$ denotes the total cost function; $b_{1} $ and $b_{2} $ denote the weight coefficients; $X_{i} $ denotes the decision variable, which is determined by the path node $P_{ij} =(x_{ij} $, $y_{ij} $, $z_{ij} )$ [14]. The UAV itself also has some constraints, such as the constraints of the range, as shown in Eq. (4).

In Eq. (4), $L_{\max } $ denotes the farthest distance for safe return; $l_{i} $ denotes the navigation distance of each segment; $n$ denotes the total navigation distance. The UAV flight speed expression in the speed constraint is shown in Eq. (5).

In Eq. (5), $V_{\max } $ refers to the maximum speed of UAV flight and $V_{\min } $ refers to the minimum speed. The expression for the height constraint of the UAV is shown in Eq. (6).

In Eq. (6), $H_{\max } $ refers to the maximum flight altitude and $H_{\min } $ refers to the minimum flight altitude. The constraints of the angle are shown in Eq. (7).

In Eq. (7), $\varphi $ denotes the corner constraint of the UAV. The minimum flight distance also needs to be controlled within a certain range, as displayed in Eq. (8).

In Eq. (8), $L$ denotes the actual flight distance. $L_{\min } $ denotes the minimum flight distance from the starting point to the target point. In summary, multiple constraints are considered in the UAV flight path planning study, including the minimum flight distance, angle limitation, height limitation, speed limitation, and range limitation ^[12].

3.2 Optimization of PSO algorithm based on population dynamic update approach

PSO algorithm is an optimization strategy based on group intelligence, by simulating the behavior of biological groups such as flocks of birds or schools of fish, the optimization problem is abstracted into the search behavior of particles in the group. The optimal solution is searched through continuous iteration ^[17]. The mathematical language is used to describe the $i$ th particle position, as displayed in Eq. (9).

In Eq. (9), $N$ denotes the population size; $d$ denotes the spatial dimension. The velocity is shown in Eq. (10).

In Eq. (10), $V_{i} $ refers to the velocity of the first $i$ particle. The updated velocity and position equation after iteration is shown in Eq. (11).

In Eq. (11), $w$ refers to the inertia weight factor; $c_{1} $ and $c_{2} $ refer to the learning factor; $r_{1} $ and $r_{2} $ satisfy the uniformly distributed random numbers of $[0$, $1]$; $X_{ipbest}^{t} $ denotes the individual optimal position. $X_{gbest}^{t} $ denotes the global optimal position. The linearly decreasing inertia weights are used to balance the global and local search ability, as displayed in Eq. (12).

In Eq. (12), $k_{\max } $ refers to the maximum iteration times. The research sets the maximum inertia weight factor $w_{\max } $ to 0.9 and the minimum inertia weight factor $w_{\min } $ to 0.4. This setting is a more reasonable choice and can provide good performance on most issues. The inertia weight within this range enables particles to have strong exploration ability in the early stage of search and good development ability in the later stage. This method adjusts the search speed according to the inertia weight factor during the search. The $w$ is displayed in Eq. (13).

In Eq. (13), $w$ has been updated. The size of $c_{1} $, $c_{2} $ shows two patterns in the pre and post period, and their change equations are shown in Eq. (14).

In Eq. (14), $c_{1i} $ denotes the initial value of $c_{1} $ and $c_{1f} $ denotes its termination value; $c_{2i} $ denotes the initial value of $c_{2f} $ and $c_{2f} $ denotes its termination value. The motion diagram of a single particle is shown in Fig. 3.

Fig. 3. Particle motion diagram.

In Fig. 3, in the particle swarm algorithm, each particle adjusts its speed relying on individual and group experience, and keeps moving iteratively to the local and global optimal position to change the state. $p^{k} $ denotes the particle position at the $k$th generation, and $V_{gbest} $ denotes the global optimal particle velocity. The study introduces a population dynamic update method to address the issue that the population is prone to fall into the local optimum. When the particle population in the global has not been updated many times, the global optimal particle velocity is changed at this time, as shown in Eq. (15).

In Eq. (15), $\alpha $ and $\beta $ are constants; $r$ refers to a random number in $[0,1]$; and $p_{j}^{m} $ denotes the $j$ th particle position when it gets the global optimum ^[18]. The operation flow of PSO algorithm is organized based on the above improvement steps, as displayed in Fig. 4.

In Fig. 4, first, the algorithm starts with initializing the velocity and position of the particle swarm. Next, the inertia weights $w$are calculated and based on this the fitness of each particle is determined. Once this step is completed, the new individual and population optimal solutions are determined. Subsequently, the velocity and position are updated. During this process, if it is detected that the optimal solution of the population has not been updated in $p$ consecutive iterations, the algorithm will use Eq. (15) to adjust the global optimal velocity, thus avoiding the particles from falling into the local optimum. Eventually, when the iteration reaches the preset upper limit, it will output the final global optimal solution.

Fig. 4. Improved PSO algorithm flowchart.

3.3 Path planning strategy based on A* with improved PSO

To plan the flight path, the choice of the fitness function affects the algorithm to get the optimal value. In the standard algorithm, the fitness function is only related to the length of the flight path. However, according to the real environment of UAV flight, when there are a lot of turns during the flight, it will affect the flight. Therefore, the paper adds the corner of the flight $\zeta $ to the fitness function, as shown in Eq. (16).

In Eq. (16), Direction $\phi _{(p-1,p)} =x_{p-1} -x_{p} $; $\zeta _{(p-1,p)} =y_{p-1} -y_{p} $; $\gamma _{(p-1,p)} =z_{p-1} -z_{p} $. The path and corner functions together form the fitness function shown in Eq. (17).

In Eq. (17), $l$ denotes a linear function and $l$ is proportional to the number of corners. A* Algorithm is an algorithm for path planning in the graph plane with multiple nodes and is commonly used to solve shortest path problems. It aims to reach the goal point at minimum cost by exploring all possible paths from the initial point. When dealing with large-scale search spaces, the A* algorithm may not be efficient enough because it may need to explore many nodes in the graph. Especially when the goal node is far from the start node, the algorithm takes longer to find the best path. During conventional A* algorithms, path planning often only considers extending the search from the current node to the direct neighboring nodes of the target node, as shown in Fig. 5.

Fig. 5. The limitations of the A* algorithm.

Fig. 5 illustrates the limitations of standard raster maps in UAV path planning. Due to its fixed grid structure, the turn angles are limited and usually vary between the number of directions that can be traveled. This results in paths generated by the traditional A* containing redundant nodes and impractical turns, making the paths non-optimal. In practice, this limitation may result in planning paths that are longer than the shortest paths that are actually reachable. When faced with the task of planning paths for multiple independent mobile bodies at the same time, the A* algorithm alone may no longer be applicable and needs to be combined with other methods. The UAV logistics and distribution path planning based on the A* algorithm with the improved PSO is displayed in Fig. 6.

Fig. 6. Process of UAV logistics distribution path planning model based on A* and improved PSO.

In Fig. 6, the study proposes an improved method that combines the A* algorithm with the PSO. The main idea is to first perform a local path search using the A* algorithm to determine the partial shortest paths. Then, the positions of the particles are determined based on each set of path points. Finally, an iterative search is performed globally using the improved particle swarm algorithm until the globally optimal path is found. The advantage of this approach is that it can fully utilize the exact path search property of the A* algorithm, while at the same time find the global optimal solution in a large-scale environment through the fast convergence of the PSO. Therefore, combining these two algorithms can accomplish the task of path planning for global environments more effectively.

4.1 Performance test of improved PSO algorithm

To lower the experimental error, the study uses the same device to start the testing and validation experiments, the operating system of the experiments is CentOS 7, Python version is 3.8, Pytorch framework is used, CUDA version is 11.4, GPU is NVIDIA RTX 3090 Ti, CPU is Intel Core i9-12900K, and memory is 128GB. 128GB, and 40 GB of video memory. The experiments use the OpenStreetMap (OSM) dataset to obtain city map data. These data usually contain information such as road networks, buildings, geographic coordinates, etc., which are applied to simulate UAV path planning in urban environments. To verify the superiority of the PSO for research optimization, PSO, Chaotic PSO (CPSO) and Distributed PSO (DPSO) are introduced for comparative analysis. As the iteration increases, the fitness values of the four methods are displayed in Fig. 7.

Fig. 7. Changes in fitness values of four methods.

In Fig. 7, the PSO algorithm without any method of optimization has the slowest convergence speed and higher fitness value. CPSO converges at 70 iterations with a value of 172.85, DPSO converges at 75 iterations with a value of 161.62, and the improved PSO algorithm of the study converges at 50 iterations with a value of 134.38. In summary, compared to the remaining two models, the proposed method converges faster and its performance is significantly better than DPSO and CPSO. Fig. 8 shows the fit of the research PSO algorithm.

Fig. 8. Fitting of the network.

Fig. 8(a) exhibits how well the predicted values of the PSO algorithm fit the actual expected values on the training set; Fig. 8(b) illustrates the fitting condition on the validation set; Fig. 8(c) represents the matching condition of the predicted with the expected values on the test set; and, finally, the correlation between the predicted outputs and the target values on the overall training dataset is given in Fig. 8(d). Among them, 'Y=T' indicates that the predicted value of the PSO algorithm is equal to the actual expected value. Usually, the algorithm is considered to have a good performance when the determination coefficient R${}^{2}$ exceeds 0.85. As seen in Fig. 8(d), the determination coefficient R${}^{2}$ reaches 0.93, which implies that the PSO algorithm has the ability to predict with high accuracy.

4.2 Comparative analysis of UAV path planning models based on improved PSO and A*

To verify the constructed model, three advanced models in the field of intelligent path planning and optimization were selected for comparative analysis. Introduce three models: Competitive Learning based Grey Wolf Optimizer (CLGWO), Dynamic Rapid Exploration Random (DRRT), and TreeExploration Enhanced Rapid Exploration Random Tree (ERT). Among them, CLGWO simulates the social hierarchy and hunting strategy of gray wolves for search optimization, combined with competitive learning mechanisms, to enhance the algorithm's global search ability and local development ability. DRRT utilizes a randomly extended tree structure for path planning, which can dynamically respond to environmental changes. By adding exploration enhancement mechanisms, ERT can quickly find the optimal path. The four models are simulated for 1000 times, and their success rate of obstacle avoidance is shown in Fig. 9.

Fig. 9. Comparison of obstacle avoidance effects among four models.

The data in Fig. 9 present the obstacle avoidance success rates of the four different models under different obstacle numbers and speed conditions. Specifically, Figs.9(a)-9(c) correspond to scenarios with 4, 8 and 12 obstacles. In these graphs, the horizontal axis presents the movement speed of the obstacles, while the vertical axis shows the corresponding success rate of obstacle avoidance. In the simulation process, if the model can successfully plan a path that can avoid all dynamic obstacles and safely reach the target position, it is considered that the model has successfully avoided obstacles. In Fig. 9(a), it is observed that all models achieve 100% obstacle avoidance success rate when the obstacle velocity is below 10 m/s. However, as the obstacle speed adds, the success rate of each model shows a decreasing trend. In this process, the research model and the DRRT model show a slower rate of decline, especially when facing high-speed moving obstacles, the obstacle avoidance success rate decreases from 100% to 90.02% in the DRRT model, while the success rate of the research model decreases to 96.23%. Further analysis of Fig. 9(b) and 9(c) reveals a similar trend to Fig. 9(a), i.e., obstacle avoidance success rate decreases in most of the models as the number and speed of obstacles increase. These observations suggest that the obstacle avoidance ability of the other three models, except the research model, is significantly affected when faced with a greater number of obstacles and higher velocities. In contrast, the research model maintained a high obstacle avoidance success rate of 94.85% even under the extreme conditions of obstacles with speeds up to 25 m/s and a number of 12. In the urban area UAV logistics path simulation, MatlabR2021a is used to construct a mountain model for testing, and the UAV flight space is set to $100$ m $\times$ $100$ m $\times$ $2000$ m, and the four types of models are applied to the three-dimensional path square. The results are presented in Fig. 10.

Fig. 10. Three dimensional diagram of four path planning methods.

In Fig. 10, all four methods mentioned above can find feasible flight paths in urban environments, and it can be seen that in the ERTT method in the middle part of the map, there is a great corner in order to avoid obstacles, compared to the research model and the DRRT model, the curves are smoother. The research model is the red line, which is the shortest and smoothest. It indicates that the research build model is more suitable for UAV path planning. To verify the applicability in urban environments, the urban environment will be simulated and the four models will be applied to simple and complex urban environments, as shown in Fig. 11.

Fig. 11. Results of simple and complex urban environmental planning.

In Fig. 11(a), in a simple urban environment, the terrain is relatively flat, the distribution of obstacles is concentrated, and there is relatively little interference with each other, which is a requirement for UAV flight to respond quickly. It can be seen that the paths of the four models are similar, none of them have collision with obstacles, but the research model planning path is the most gentle, the distance is the shortest, and it is the most close to the obstacles of the path planning. The four algorithms were subjected to 100 repetitions of experiments in a simple urban environment and a complex city. Table 1 displays the results.

In the data analysis in Table 1, the simulations for simple urban environments show that both the proposed model and the DRRT algorithm of the study effectively shorten the path lengths compared to ERRT and CL-GWO, thus improving the planning efficiency. Specifically, the DRRT algorithm presents path lengths of 1001 m, 1024 m, and 1042 m under the optimal, average, and worst conditions, respectively, with an average number of iterations of 64 and an average runtime of 94.14 ms. In contrast, the research model exhibits superior performance, realizing an optimal path length of 992 m, an average path length of 923 m, and a worst path length of 1011 m. In addition, the average iteration of the model is reduced to 45, and the runtime is shortened to 8.61 ms. In the more complex urban environment simulation, the research model continues to maintain its leading position, achieving an optimal path length of 1002 m, an average path length of 1031 m, and an average path length of 1016 m, with an average number of iterations of 58 and a runtime of 12.37 ms. These data fully demonstrate the efficiency and adaptability of the research model in urban area UAV logistics and distribution path planning.