Bi-level intelligent dynamic path planning for an UAV in low-altitude complex urban environment

Abstract

An offline and online bi-level structure-based dynamic path planning algorithm is proposed for an unmanned aerial vehicle (UAV) in low-altitude complex urban environment. First, an improved Hunger Games Search (HGS) algorithm is developed to generate an offline optimized path under the UAV’s performance constraints and the known static obstacles’ constraints. The individuals of the proposed algorithm will be divided into multiple groups to increase the population diversity. And then, a dynamic grouping strategy and a quantum-behaved behavior are proposed to solve the premature convergence’s problem and the imbalance problem between exploration and exploitation ability in HGS. To improve the dynamic obstacle avoidance efficiency of the algorithm, the dynamic obstacles are classified into three categories: newly added no-fly zone, known and unknown dynamic obstacles. Then, utilizing the information of the offline optimized path and the airborne sensors, three kinds of online planning strategies—an improved rapid-exploring random tree (RRT), a changing speed strategy, and a novel three-dimensional rolling windows—are introduced to dynamically update the path or speed of the UAV. Simulation results indicated that the improved HGS can enhance the performances of the traditional HGS and outperform other compared algorithms on the benchmark functions. Meanwhile, the online planning strategies can effectively achieve dynamic obstacle avoidance within the constraints of offline path. More specially, the planning time and angles of the local path to avoid the no-fly-zone’s influence are improved by 11.3% and 56.8% through utilizing the improved RRT.

Introduction

In the past few decades, much attentions have been paid to the applications of the unmanned aerial vehicle (UAV) in military and civil domains, such as the target tracking (Yao et al., 2015), logistics delivery (Song et al., 2018), and environment monitoring (Shen et al., 2004). During the missions, autonomous flight technologies have played a very important role (Li et al., 2022; Yang et al., 2015). To some extent, the realization of autonomy and intelligence mainly depends on the flight control and path planning techniques (Li et al., 2020). And the successful execution of the UAV’s missions is up to the reasonability and effectiveness of the planned path. In low-altitude complex urban environment, there exist not only various static obstacles, such as buildings and trees, but also a variety of dynamic obstacles, such as birds, cooperative/noncooperative agents or UAVs. The security and effectiveness is threatened by the incomplete environment information and limited detection range of airborne sensors, which also bring great challenges to path planning problem with obstacle avoidance for the UAV (Zhou et al., 2021). Therefore, it is crucial to plan a collision-free flight path with lower energy consumption under the related constraints in complex and changing environments. Two requirements should be satisfied for the path planning of the UAV in low-altitude complex environment: one is to find the optimal path under the UAV’s performance constraints and the constraints of the known static obstacles, and the other is to avoid the dynamic obstacles’ threat under the optimal path constraint, the UAV’s performance constraints and the sensors’ information.

According to the time-domain type, path planning with obstacle avoidance of the UAV can generally be classified into two categories, that is, offline planning and online planning (Zhao et al., 2018). If the information of global environment is completely known to the UAV, the problem is known as the offline path planning. The offline planning focus on the optimality and stability of the generated path. It is a complex global optimization problem with multiple specific constraints essentially. A series of methods have been proposed to deal with this complex global optimization problem, such as sampling-based techniques (Sucan and Kavraki, 2012), artificial neural networks (Duan and Huang, 2014), and heuristic methods (Mac et al., 2016). Whereas, it is difficult to consider the related constraints of the UAV into these above algorithms. Among all kinds of path planning technologies, swarm intelligence optimization algorithms (SIOAs) can conquer these difficulties. It is very effective to deal with the complex global optimization problems including offline planning problem (Heidari and Pahlavani, 2017). “No free lunch theorem” has pointed out that no algorithm can solve all types of optimization problems, and each algorithm has its superiority only for some specific problems (Wolpert and Macready, 1997). Therefore, a variety of new or improved algorithms are proposed to obtain more promising results with respect to different problems (Chrouta et al., 2018; Liu et al., 2020; Yang et al., 2021). Yang et al. (2021) indicated that there exists a research gap in all previous SIOAs, and it compels users to focus on improving operations of algorithms based on specific evolutionary process. A general-purpose algorithm named Hunger Games Search (HGS) was designed to focus more on performance rather than metaphor change. Moreover, it has already been proved to be well-performed when solving some engineering problems (Nguyen and Bui, 2021; Onay and Aydemır, 2022). However, offline planning cannot guarantee the security in the uncertain environment. The calculation time and overhead of them always increase exponentially with the scale and complexity of the environments. Thus, it is impossible to dynamically update the path of the UAV by only utilizing offline planning.

Online planning is a dynamic multi-objective optimization problem (Zhao et al., 2018), which is defined under the condition that the global environment information is partially or completely unknown. Some algorithms are widely used to address the online planning problem, such as rolling windows method (Xi and Zhang, 2003), rapid-exploring random tree (RRT) (Huang and Sun, 2020), artificial potential field (APF) (Lin et al., 2021) and so on. Xi and Zhang (2003) proposed a rolling windows method based on local information for dynamic path planning problem of mobile robots. To solve the real-time tracking path planning problem for a solar-powered UAV with the constraint of energy-optimal, a planning method based on particle swarm optimization (PSO), rolling optimization and model predictive control (MPC) was proposed in Huang et al. (2016). Huang and Sun (2020) proposed a greedy strategy–based bi-direction RRT algorithm for the UAV. They aimed to replan path in real-time according to the complex dynamic environment information. Online planning algorithms can make quick reactions to the changes of environmental information. Nevertheless, they usually overlook the optimization of the whole environment, which will lead to the failure of finding the ideal or feasible path.

Therefore, multiple hybrid path planning algorithms which combine offline and online planning are proposed to overcome the limitations or constraints of them. They aim to realize dynamic path planning of the UAV in low-altitude dynamic environment (Aldao et al., 2022; Chen et al., 2020; Elmokadem and Savkin, 2021; Wu and Hu, 2020). In Chen et al. (2020), a hybrid algorithm based on A* and cubic spline method was proposed to solve the path planning problem for the UAV in a dynamic environment. But the speeds of the dynamic obstacles are set to a constant value, and this algorithm has been studied in the two-dimensional environment. Wu and Hu (2020) proposed a A* and RRT based hybrid algorithm to realize dynamic planning in low-altitude dynamic urban area. And it pointed out that it was necessary to classify the dynamic obstacles and design corresponding online planning strategies for the UAV. In Elmokadem and Savkin (2021), a hybrid algorithm based on an improved RRT and a sliding mode control’s reactive control law based was proposed. It aims to realize the UAV’s autonomous navigation in low-attitude dynamic environment with partially unknown information. Aldao et al. (2022) presented a hybrid obstacle avoidance algorithm based on A* and optimal control to address the path planning problem of the UAV in dynamic building environments. First, the above references didn’t consider the performance constraints of the UAV comprehensively. The generated path usually contains many sharp turns which will be hard for the UAV to track, ad it increases the energy consumption. Second, only the flight direction and attitudes will be adjusted to avoid the threat of the dynamic obstacles in the most of the existing literatures. The constant speed assumption of them will restrict the flexibility and effectiveness of the UAV to avoid the dynamic obstacles’ influence. Moreover, most of the existing literatures didn’t classify the dynamic obstacles according to their different characteristics. A universal planning strategy is utilized to avoid all kinds of dynamic obstacles’ threat. However, single strategy cannot guarantee the effectiveness of online planning in terms of computation time, obstacle avoidance and energy consumption.

Inspired by the existing results, this paper proposes a novel bi-level structure algorithm to address the dynamic path planning problem for the UAV in low-altitude complex environment. First, an improved HGS is utilized to generate an offline optimized path. Then, utilizing the offline optimized path’s information and the airborne sensors’ information, three kinds of online planning strategies are proposed to avoid the threat of dynamic obstacles by updating the path or speed of the UAV in real-time. The main contributions are described as follows:

1. An improved HGS based on dynamic grouping strategy (noted as DGSHGS) is proposed to generate an offline optimized path for the UAV under its performance constraints and the static obstacles’ constraints. A multiple-groups division strategy, a dynamic grouping strategy, and an improved quantum-behaved foraging behavior are developed in DGSHGS to solve the premature convergence problem and the imbalance problem between exploration and exploitation ability in the traditional HGS.

2. To improve the dynamic obstacle avoidance efficiency, the dynamic obstacles are classified into three categories: newly added no-fly zone, known and unknown dynamic obstacles. Then, utilizing the information of the offline optimized path which is generated by DGSHGS and the airborne sensors’ information, three kinds of online planning strategies—an improved rapid-exploring random tree (RRT), a changing speed strategy, and a novel three-dimensional rolling windows—are proposed to achieve dynamic obstacle avoidance. After avoiding the threat of no-fly zone and unknown dynamic obstacles, the generated local optimal path will return to the offline optimized path as soon as possible to reduce both the modification degree of the offline path and the calculation time. More especially, the planning time and path angles of local path to avoid the no-fly-zone’s threat are improved by 11.3% and 56.8%, respectively, compared with the variation of RRT in Wu and Hu (2020) under the improved RRT based online planning strategy.

Path planning model of the UAV

During the flight process in low-altitude urban areas, the security of the UAV is influenced by its performance constraints, the external environment constraints and the airborne sensors’ limited detection range. Therefore, it is crucial to plan a collision-free flight path with lower energy consumption under the related constraints in complex dynamic environments. In this section, the planning environment, the related constraints, and the cost function of path planning are modeled as follows.

Environment modeling

Environment modeling is the foundation before path planning of the UAV. In this section, the grid-based method is utilized to build the physical model of flying space. First, a three-dimensional rectangular coordinate system is established, with the vertex in the bottom-left of the map as the origin point. Then, the space of three-dimensional path planning is obtained by taking the maximum length of the map along the coordinate axes. As shown in Figure 1, m planes ${\mathrm{\Pi }}_p$ $(p=1,2,\dots ,m)$ are obtained by dividing the space equally along the edge AD with the grid length $\mathrm{\Delta }d$ . For any of the above planes, this section divides it into m equal parts along edge AB and n equal parts along edge $\mathrm{AA}\prime$ . The path planning space is divided into grids of equal size, and the generated path points are located on the endpoint of grids. According to this structure, the definition of path points can be expressed as the serial number coordinates $P(o,p,q)(o=1,2,...,l;p=1,2,..,m;q=1,2,..,n)$ and the position coordinates $P(x,y,z)$ .

Figure 1.

Three-dimensional space division.

Problem formulation

UAV performance constraints

The generated path points must meet the related constraints to ensure the security of the UAV. To generate a more feasible path, the related angles or attitudes constraints must be taken into account. The mathematical equations of the $j\mathrm{th}$ path point’s maximum turning angle ( ${\psi }^{\max }$ ) and maximum climbing angle ( ${\varphi }^{\max }$ ) constraints are shown as equation (1) (Phung and Ha, 2021) and equation (2)

$$0\le {\psi }_{j-1,j,j+1}=\left|\arccos \left(\frac{{\mathit{P}}_{j-1}^{\prime }{\mathit{P}}_j^{\prime }·{\mathit{P}}_j^{\prime }{\mathit{P}}_{j+1}^{\prime }}{\Vert {\mathit{P}}_{j-1}^{\prime }{\mathit{P}}_j^{\prime }\Vert ·\Vert {\mathit{P}}_j^{\prime }{\mathit{P}}_{j+1}^{\prime }\Vert }\right)\right|\le {\psi }^{\max }$$

(1)

$$0\le {\varphi }_{j-1,j,j+1}=|\frac{\pi }{2}-\arcsin \left(\frac{\mathbf{\tau }·\mathit{n}}{\Vert \mathbf{\tau }\Vert ·\Vert \mathit{n}\Vert }\right)|\le {\varphi }^{\max }$$

(2)

where the coordinate of the $j\mathrm{th}$ path point is represented as $P_j(x_j,y_j,z_j)(j=1,2,...,N_{\max }-1)$ . The turning angle ${\psi }_{j-1,j,j+1}$ is the angle between projection vectors of consecutive path segment ${\mathit{P}}_{j-1}{\mathit{P}}_j$ and ${\mathit{P}}_j{\mathit{P}}_{j+1}$ on the horizontal plane in Figure 2. The projection vectors are noted as ${\mathit{P}}_{j-1}^{\prime }{\mathit{P}}_j^{\prime }$ and ${\mathit{P}}_j^{\prime }{\mathit{P}}_{j+1}^{\prime }$ . The climbing angle ${\varphi }_{j-1,j,j+1}$ is the angle between path segment ${\mathit{P}}_j{\mathit{P}}_{j+1}$ and its projection vector ${\mathit{P}}_j^{\prime }{\mathit{P}}_{j+1}^{\prime }$ . $N_{\max }$ represents the maximum index of path points. $\mathbf{\tau }=\frac{{\mathit{P}}_j{\mathit{P}}_{j+1}}{\Vert {\mathit{P}}_j{\mathit{P}}_{j+1}\Vert }$ , and $\mathit{n}=[0,0,1]$ .

Figure 2.

Turning and climbing angle calculation.

Moreover, the minimum distance that UAV can fly is set as the minimum length of path segment $l_{\min }$ . The above constraint is described as

$$\forall l_h\ge l_{\min },h=1,2,...,N_{\max }-1$$

(3)

where, $l_h$ represents the length of the $h\mathrm{th}$ path segment.

Environment constraints

Since the UAV is not a mass point but rather a rigid body that occupies space, it is important to ensure that the generated path points and the path segments keep some certain distances from the obstacles. Inspired by Wu and Hu (2020), the position constraints of the path points and the path segments are modeled as

$$\sqrt{(x_j-\overline{x})+(y_j-\overline{y})+(z_j-\overline{z})}\ge d_{\mathrm{safe}}$$

(4)

$$\frac{\Vert \mathbf{\tau }\times {\mathit{P}}_j{\mathit{P}}_o\Vert }{\Vert \mathbf{\tau }\Vert }\ge d_{\mathrm{safe}}$$

(5)

where $P_o=(\overline{x},\overline{y},\overline{z})$ represents points randomly selected in the space which is occupied by the obstacles. Any points located in the space which is occupied by the obstacles need to keep a safe distance from the generated path points. The value of $d_{\mathrm{safe}}$ is set as 2 $\times \mathrm{\Delta }d$ .

Overall cost function

In this paper, planning strategies are proposed to update the UAV’s speed or path to avoid the obstacles’ threat. It is assumed that the UAV’s speed remains constant in changing path strategies. Therefore, the path length can be considered as the evaluation index of the energy consumption. Inspired by Qu et al. (2020), the number of the path points and the angles are also considered as the evaluation indexes of the energy consumption. Based on the related constraints above, a novel path quality cost function J is shown as

$$\begin{array}{c}\multicolumn{1}{c}{J=k_1·\sum _{j=1}^{N_{\max }-1}\parallel P_{j+1}-P_j\parallel }\\ \multicolumn{1}{c}{+k_2·\mathrm{turns}+k_3·\sum _{j=2}^{N_{\max }-1}({\psi }_{j-1,j,j+1}+A)}\\ \multicolumn{1}{c}{+k_4·\sum _{j=3}^{N_{\max }-1}(|{\varphi }_{j-1,j,j+1}-{\varphi }_{j-2,j-1,j}|+A)}\\ \multicolumn{1}{c}{+k_5·\sum _{j=1}^{N_{\max }-1}B}\end{array}$$

(6)

$$\mathrm{s}.\mathrm{t}.\{\begin{array}{l}0\le \left|\arccos \left(\frac{P_{j-1}^{\prime }{P}_j^{\prime }\cdot P_j^{\prime }{P}_{j+1}^{\prime }}{\Vert P_{j-1}^{\prime }{P}_j^{\prime }\Vert \cdot \Vert P_j^{\prime }{P}_{j+1}^{\prime }\Vert }\right)\right|\le {\psi }^{max}\\ 0\le \Vert \frac{\pi }{2}-\arcsin \left(\frac{\tau \cdot n}{\parallel \tau \parallel \cdot \parallel n\parallel }\right)\Vert \le {\varphi }^{max},\tau =\frac{P_j{P}_{j+1}}{\Vert P_j{P}_{j+1}\Vert },n=\left[0,0,1\right]\\ \forall l_h\ge l_{\min },h=1,2,\cdots ,N_{max}-1\\ \sqrt{\left(x_j-\overline{x}\right)+\left(y_j-\overline{y}\right)+\left(z_j-\overline{z}\right)}\ge d_{safe},j=1,2,\dots ,N_{\max }\\ \frac{\Vert \tau \times P_j{P}_o\Vert }{\Vert \tau \Vert }\ge d_{safe},P_o=\left(\overline{x},\overline{y},\overline{z}\right)\end{array}$$

(7)

where turns is the number of generated path’s turning points. If the angle constraints are not satisfied during the mission execution process, A will be set as a relatively large number. Otherwise, A will be set as a relatively small number to balance the proportion of the angle part. B represents the threat of the obstacles, and it will be set as a relatively large number while the distance between the $j\mathrm{th}$ path segment and the obstacle is smaller than $d_{\mathrm{safe}}$ . $k_1$ , $k_2$ , $k_3$ , $k_4$ and $k_5$ refer to the weights of different index, respectively.

Proposed method

An offline and online bi-level structure algorithm is proposed in this paper to solve the dynamic planning problem of the UAV in low-altitude complex urban areas. First, an improved HGS is proposed to generate an offline optimized path under the static obstacles’ constraints. And then three online planning strategies are proposed to avoid three kinds of dynamic obstacles under the offline path’s constraints, the UAV’s performance constraints, and the sensors’ information.

Offline path planning strategy

Standard hunger games search algorithm

With the path quality evaluation function defined in equation (7), the offline planning problem becomes a complex global optimization problem which aims to find the path that minimizes this function. In this section, HGS is adopted to realize the offline planning due to its well performance in solving the complex optimization problems. Inspired by the activities driven by hunger and behavioral choice of social animals, the search games are categorized into two types in HGS. The criteria of classification are whether social animals are involved in the cooperation or not during foraging (Clutton-Brock, 2009). For a d-dimensional optimization problem, the HGS algorithm will create a population which is composed by several random individuals for initialization. Then, it constantly finds the optimal candidate solution according to the evaluation function value during the iteration course. The mathematical models of HGS are proposed in Yang et al. (2021):

1. Approaching food. The mathematical updating model of social animals’ foraging position is expressed as

$$X_i^d(t+1)=\{\begin{array}{c}\mathrm{Gam}{\mathrm{e}}_1:X_i^d(t)·(1+\mathrm{randn}(1)),\hfill \\ \multicolumn{1}{c}{r_1<u}\\ \multicolumn{1}{c}{\mathrm{Gam}{\mathrm{e}}_2:W_{1i}^d·X_b^d(t)+R_i^d·W_{2i}^d·|X_b^d(t)-X_i^d(t)|,}\\ \multicolumn{1}{c}{r_1>u,r_2>E}\\ \multicolumn{1}{c}{\mathrm{Gam}{\mathrm{e}}_3:W_{1i}^d·X_b^d(t)-R_i^d·W_{2i}^d·|X_b^d(t)-X_i^d(t)|,}\\ \multicolumn{1}{c}{r_1>u,r_2<E}\end{array}$$

(8)

where $X_i^d(t)(i=1,2,\dots ,N;d=1,2,\dots ,D)$ is the $d\mathrm{th}$ dimension of the $i\mathrm{th}$ individual’s position vector at the current iteration, and $X_b^d(t)$ is the $d\mathrm{th}$ dimension of the best individual’s position vector. $\mathrm{randn}(1)$ is a random number satisfying normal distribution. $R_i^d=2·\mathrm{shrink}·\mathrm{rand}-\mathrm{shrink}$ , $\mathrm{shrink}=2·(1-t/T_{\max })$ and $T_{\max }$ is the maximum iterations number. $r_1$ , $r_2$ and rand are random numbers in $[0,1]$ , and u is a constant number that improves the algorithm’s performance. E is a variation control parameter that acts on all positions. It is expressed as $E=\mathrm{sech}(|J(i)-\mathrm{BF}|)$ . $J(i)$ is the evaluation function value of the $i\mathrm{th}$ individual; BF is the best evaluation function value in the current iteration process (so far), and $\mathrm{sech}(x)=2/(e^x+e^{-x}).$

As it can be seen in equation (8), the search pattern is switched randomly with the parameters u and E to enrich the flexibility of HGS. The first pattern expresses the activities of a few social animals, which search for food independently. ${\mathit{X}}_i(t)·(1+\mathrm{randn}(1))$ represents that the $i\mathrm{th}$ individual sets the current position as the center point. ${\mathit{X}}_i(t)=(X_i^1(t),X_i^2(t),\dots ,X_i^D(t))$ is the $i\mathrm{th}$ individual’s position vector. Then, it floats in the range of ${\mathit{X}}_i(t)·\mathrm{randn}(1)$ when searching for food. The second pattern indicates how an agent can cooperate with others during the foraging course. According to the notifications from peers, the $i\mathrm{th}$ individual of the population arrived at ${\mathit{W}}_{1i}\odot {\mathit{X}}_b(t)$ , and then searched for food in a floating range of $\pm {\mathit{R}}_i\odot {\mathit{W}}_{2i}\odot |{\mathit{X}}_b(t)-{\mathit{X}}_i(t)|$ . The position vector of the best individual is ${\mathit{X}}_b(t)=(X_b^1(t),X_b^2(t),\cdots ,X_b^D(t))$ . The operational symbol ⊙ represents the Hadamard product. $|{\mathit{X}}_b(t)-{\mathit{X}}_i(t)|$ simulates the activity scope of the $i\mathrm{th}$ individual at the current iteration. The weight coefficient vectors are expressed as ${\mathit{W}}_{1i}=(W_{1i}^1,W_{1i}^2,\dots ,W_{1i}^D)$ and ${\mathit{W}}_{2i}=(W_{2i}^1,W_{2i}^2,\dots ,W_{2i}^D)$ , respectively. ${\mathit{R}}_i$ is introduced to control the activity range of the $i^{\mathrm{th}}$ individual. It is expressed as ${\mathit{R}}_i=(R_i^1,R_i^2,\dots ,R_i^D)$ . The variation range of ${\mathit{R}}_i$ declines in a linear speed, which ensures the ability of the algorithm to keep balance on global exploration and local exploitation.

1. Hunger role. In equations (9) and (10), ${\mathit{W}}_{1i}$ is the deviation in obtaining the actual food location of the $i^{\mathrm{th}}$ individual; ${\mathit{W}}_{2i}$ simulates the negative and positive influence of hunger on the $i^{\mathrm{th}}$ individual’s activity scope. These two coefficient vectors act as perturbations to prevent HGS from falling into local extremum. Each dimension of the ${\mathit{W}}_{1i}$ and ${\mathit{W}}_{2i}$ can be presented in the form of

   $$W_{1i}^d=\{\begin{array}{c}\mathrm{hungr}{y}_i·\frac{N}{\mathrm{SH}\mathrm{ungry}}·r_4,r_3<u\hfill \\ \multicolumn{1}{c}{1,r_3>u}\end{array}$$

   (9)
   $$W_{2i}^d=(1-\exp (-|h{g}_i-\mathrm{SH}\mathrm{ungry}|))·r_5·2$$

   (10)

where SHungry is the total hunger extent of the population; $r_3$ , $r_4$ , and $r_5$ are random numbers in the range between 0 and 1. $\mathrm{hungr}{y}_i$ simulates the hunger extent of the $i^{\mathrm{th}}$ individual, which is demonstrated as

$$\mathrm{hungr}{y}_i=\{\begin{array}{c}0,\mathrm{AllFitnes}{s}_i=\mathrm{BF}\hfill \\ \multicolumn{1}{c}{\mathrm{hungr}{y}_i+H,\mathrm{AllFitnessi}!=\mathrm{BF}}\end{array}$$

(11)

In equation (11), $\mathrm{AllFitnes}{s}_i$ presents the $i\mathrm{th}$ individual’s evaluation function value at the current iteration. During the evolutionary process, the hunger extent of optimal candidate solution is set as $0$ . For the rest of the individuals, a new variable H is defined as

$$H=\{\begin{array}{c}\mathrm{LH}·(1+r),\mathrm{TH}<\mathrm{LH}\hfill \\ \multicolumn{1}{c}{\mathrm{TH},\mathrm{TH}\ge \mathrm{LH}}\end{array}$$

(12)

$$\mathrm{TH}=\frac{J(i)-\mathrm{BF}}{\mathrm{WF}-\mathrm{BF}}·2·r_6·(\mathrm{UB}-\mathrm{LB})$$

(13)

where WF refers the worst evaluation function value in the current iteration process (so far); $\frac{J\left(i\right)-\mathrm{BF}}{\mathrm{WF}-\mathrm{BF}}$ indicates the ratio of hunger; r and $r_6$ are random numbers in the range of $[0,1]$ . In addition, $2·r_6$ refers the negative or positive environmental influence on individual’s hunger; LH is a constant number; LB and UB represent the lower and upper bounds of the solution space. The introduce of H in equation (12) takes the different function value of individuals into consideration. It can enhance the accuracy of the algorithm to simulate the foraging behavior of social animals.

Proposed improved HGS algorithm (dynamic grouping strategy based HGS algorithm)

It is known that keeping balance on the abilities of global exploration and local exploitation is of great importance to the SIOA. Global exploration represents the ability to explore more feasible solution space. It maintains the diversity of the candidate solutions and avoids the algorithm trapping into the local optimum. Local exploitation is closely related to algorithm’s convergence speed. It indicates the abilities of individuals to search elaborately around local areas. HGS has the problem of premature convergence and the imbalance problem between exploration and exploitation ability. Therefore, an improved dynamic grouping strategy based hunger games search algorithm called DGSHGS is proposed in this section to improve the performance of HGS:

• 1. Dynamic grouping strategy. The whole population is randomly divided into H subgroups with K individuals each in DGSHGS to enrich the diversity of the population. Each subgroup utilizes its own individuals’ information to obtain more promising foraging positions. In addition, the dynamic grouping strategy is proposed to promote information exchanges between different groups. That aims to avoid the algorithm from trapping into local optimums. The dynamic grouping strategy occurs according to a fixed probability P. The triggering condition is $r_7<P$ , and $r_7$ represents a random number in the range of $[0,1]$ . Then, two random arrays ${\mathit{a}}_1$ and ${\mathit{a}}_2$ indicating the location of the individuals who exchange information in different groups are generated as

$${\mathit{a}}_1=[m_1,n_1];{\mathit{a}}_2=[m_2,n_2]$$

(14)

with $m_1=\mathrm{randi}(H),m_2=\mathrm{randi}(H),n_1=\mathrm{randi}(K),n_2=\mathrm{randi}(K),m_1\ne m_2$ . These operations increase the population diversity of some groups to improve the exploration ability. Furthermore, they reduce the individuals’ gap of some groups to increase convergence accuracy.

• 2. Nonlinear control parameter. In the early search stage, $R_i^d$ is greater than one. It ensures the broadly exploration of the feasible search space, which increases the ability to find global optimal solution. In the later stage, $R_i^d$ is less than one. The population’s search scope will be narrowed, and the individuals exploit finely in local areas to enhance the convergence accuracy. Whereas, for different optimization problems, the linear control parameter shrink is difficult to adapt to the actual optimization search process. Inspired by the setting of inertia weight in the PSO (Poli et al., 2007), a nonlinear control parameter based on anti-sinusoidal function is modeled as

  $$\mathrm{shrin}{k}_1=2-\frac{4}{\pi }·\arcsin \left({\left(\frac{t}{T_{\max }}\right)}^2\right)$$

  (15)

The attenuation speed of $\mathrm{shrin}{k}_1$ is slow in the early stage, which can increase the exploration ability of the algorithm. The decreasing speed of $s{k}_1$ is fast in the later stage, that can improve convergence speed of the algorithm. Correspondingly, $R_i^d$ may still be greater than one in the later stage, which ensures the ability of global search. Thus, the proposed control parameter $\mathrm{shrin}{k}_1$ can balance the ability of exploration and exploitation of the algorithm more effectively:

• 3. Novel foraging operations. The optimal individual’s information in corresponding group will be considered into the individual’s position updating formula. That helps to improve the convergence accuracy and the ability to avoid the local optimum of the algorithm. According to the definition of E in equation (8), it is noted that the value of E is large when the individuals’ gap is small. Thus, the algorithm will be inclined to execute the $\mathrm{Gam}{e}_3$ when the individuals becoming more similar. In order to avoid the algorithm from trapping into the local optimum, $\mathrm{Gam}{e}_2$ and $\mathrm{Gam}{e}_3$ are combined into one formulation. Inspired by Zou et al. (2014), the quantum-behaved cooperation operation is adopted to enhance the algorithm’s performance. Furthermore, E is replaced by the adaptive probability parameter $\sqrt{t/T_{\max }}$ to increase the flexibility of the algorithm. A novel updating formula of individuals’ foraging position is proposed as

$$X_e^{\mathrm{wd}}(t+1)=\{\begin{array}{c}\mathrm{Gam}{e}_1:X_e^{\mathrm{wd}}(t)·(1+\mathrm{randn}\left(1\right)),r_1<u\hfill \\ \multicolumn{1}{c}{\mathrm{Gam}{e}_2:\{\begin{array}{c}a(t)-\beta ·|\mathrm{Groupte}{n}^{\mathrm{wd}}(t)-X_e^{\mathrm{wd}}(t)|\hfill \\ \multicolumn{1}{c}{·\ln \left(\frac{1}{Q}\right),}\\ \multicolumn{1}{c}{r_1>u,r_2>\sqrt{t/T_{\max }},f<0.5}\\ \multicolumn{1}{c}{a(t)+\beta ·|\mathrm{Groupte}{n}^{\mathrm{wd}}(t)-X_e^{\mathrm{wd}}(t)|}\\ \multicolumn{1}{c}{·\ln \left(\frac{1}{Q}\right),}\\ \multicolumn{1}{c}{r_1>u,r_2>\sqrt{t/T_{\max }},f>0.5}\end{array}}\\ \multicolumn{1}{c}{\mathrm{Gam}{e}_3:W_{1e}^{\mathrm{wd}}·\mathrm{Groupbes}{t}^{\mathrm{wd}}(t)}\\ \multicolumn{1}{c}{+R_e^{\mathrm{wd}}·c·(1-W_{2e}^{\mathrm{wd}})·|p_{\mathrm{gb}}^d(t)-X_e^{\mathrm{wd}}(t)|,}\\ \multicolumn{1}{c}{r_1>u,r_2<\sqrt{t/T_{\max }}}\end{array}$$

(16)

$$a(t)=\phi ·\mathrm{Groupbes}{t}^{\mathrm{wd}}(t)+(1-\phi )·p_{\mathrm{gb}}^d(t)$$

(17)

$$\mathrm{Groupte}{n}^{\mathrm{wd}}(t)=\{\begin{array}{c}{O}_{\frac{(K+1)}{2}}^{\mathrm{wd}}(t),\mathrm{Kis}\mathrm{odd}\hfill \\ \multicolumn{1}{c}{\frac{O_{\frac{K}{2}}^{\mathrm{wd}}(t)+O_{\frac{(K+1)}{2}}^{\mathrm{wd}}(t)}{2},\mathrm{otherwise}}\end{array}$$

(18)

It should be noted that $\mathrm{Groupte}{n}^{\mathrm{wd}}(t)(w=1,2,\dots ,H)$ is inspired by the setting of cultural tendency in the coyote optimization algorithm (COA) (Pierezan et al., 2019). The superscript w represents the $w\mathrm{th}$ group; e represents the $e\mathrm{th}$ individual of group $(e=1,2,\cdots ,K)$ ; $O^{\mathrm{wd}}(t)$ represents the $d\mathrm{th}$ dimension of the median for the $w\mathrm{th}$ group’s individuals; $\beta$ , Q, f, and $\phi$ are random numbers between 0 and 1. In the proposed novel term $\mathrm{Gam}{e}_3$ , c is a weight coefficient. $\mathrm{Groupbes}{t}^{\mathrm{wd}}(t)$ and $p_{\mathrm{gb}}^d(t)$ refer to the $d\mathrm{th}$ dimension of the optimal position of the $w\mathrm{th}$ group and the population at the current iteration course, respectively. $X_e^{\mathrm{wd}}(t)$ is the $d\mathrm{th}$ dimension of the $e\mathrm{th}$ individual’s position vector of the $w\mathrm{th}$ subgroups at the current iteration. $R_e^{\mathrm{wd}}=2·\mathrm{shrin}{k}_1·\mathrm{rand}-\mathrm{shrin}{k}_1$ , and the calculation methods of $W_{1e}^{\mathrm{wd}}$ and $W_{2e}^{\mathrm{wd}}$ are the same as equations (9) and (10), respectively. The value of $1-W_{2e}^{\mathrm{wd}}$ is randomly selected in the range of $[-1,1]$ , which provides more promising foraging positions for the social animals. From equation (16), the value of $\sqrt{t/T_{\max }}$ is increasing through the iteration process: it is smaller earlier, and quantum-behaved operation is frequently implemented to increase the information exchanges between individuals. It becomes larger in the later stage, and all individuals in the population approach to the global optimal solution to improve the convergence speed. These operations are of benefit to keep balance between global exploration and local exploitation. Inspired by COA, the birth and death operations are introduced into DGSHGS to improve the quality of candidate solutions

$$\mathrm{new}\_\mathrm{ind}{}^{\mathrm{wd}}(t)=\{\begin{array}{c}{X}_{e_1}^{\mathrm{wd}}\left(t\right),\mathrm{rand}<P_s\mathrm{or}d=1\hfill \\ \multicolumn{1}{c}{X_{e_2}^{\mathrm{wd}}\left(t\right),\mathrm{rand}\ge P_s+P_a\mathrm{or}d=2}\\ \multicolumn{1}{c}{R_d,\mathrm{otherwise}}\end{array}$$

(19)

Equation (19) denotes that the birth of a new individual that is affected by the environment and the randomly selected parents in the $w\mathrm{th}$ group. The ratios of them are represented as the scatter probability $(P_s)$ and the association probability $(P_a)$ . $e_1$ and $e_2$ denote the index of random parent in the $w\mathrm{th}$ group; rand is random number in the range of $[0,1]$ ; $R_d$ is a random number inside the solution space. In order to keep the number of individuals being unchanged, the individual with the worst evaluation function value in each group will die. If there are more than one individual meeting the condition, their ages will be compared. The pseudo-code of DGSHGS is demonstrated in Algorithm 1.

Algorithm 1: Pseudo-code of DGSHGS

1 Initialize the parameters $u,T_{\max },N,D,\mathrm{SH}\mathrm{ungry}$ 2 Initialize the positions of Individuals ${\mathit{X}}_i(t)(i=1,2,...,N)$ 3 white $t\le T_{\max }$ do 4 Update ${\mathit{p}}_{\mathit{gb}}(t)$ 5 for eachGroup do 6  Calculate the evaluation value of all Individuals 7  Update $\mathrm{BF},\mathrm{WF},\mathit{Groupbest}(t),\mathit{Groupten}(t)$ 8  for eachIndividual do 9   Calculate the hungry by equation (11) 10   Calculate the ${\mathit{W}}_1$ by equation (9) 11   Calculate the ${\mathit{W}}_2$ by equation (10) 12   Update $\mathit{R}$ 13   Update positions by equation (16) 14  end 15 end 16 Birth and death by equation (19) 17 Transition between packs by equation (14) 18  $t=t+1$ 19 end

Online path planning strategies

During the flight missions of the UAV in low-altitude complex environment, various dynamic obstacles will threat its security and other performances. In order to increase the flexibility and efficiency of realizing dynamic obstacle avoidance, obstacles encountered by the UAV will be divided into three categories: dynamic obstacle with known trajectory, fixed no-fly zone and dynamic obstacle with unknown trajectory. It is assumed that the airborne sensors are utilized to detect the dynamic changes during the flying process. The range of detection is defined as a sphere which sets the UAV as the center. In this section, R refers to the detection range of the airborne sensors, and $P(T)$ denotes the current position of the UAV. The sampling time of the flying course is defined as $T_0$ .

The obstacle with known trajectory

Since the motion law of the dynamic obstacles with known trajectory is known, the changing speed or path strategy can be selected. Considering the lower modification degree of the offline optimized path, a changing speed strategy is adopted in this paper (Wu and Hu, 2020). There are three speed expressed as high speed $(V_H)$ , middle speed $(V_M)$ , and low speed $(V_L)$ for the UAV to select. When the distance between the known obstacle and the UAV will be less than $d_{\mathrm{safe}}$ after a short period of time, the changing speed strategy will be activated. After changing the speed, the speed of the UAV will maintain constant until the next speed change. If the supplied speed cannot guarantee the security of the UAV, its velocity will be reduced to zero. And, it will hover in the original position and wait for the next sampling time to judge again. If there exists more than one speed that can satisfy the dynamic obstacle avoidance’s requirement, the lower one will be selected to reduce energy consumption. The threat of the obstacle will be avoided successfully by the UAV until the distance between them is always no less than $d_{\mathrm{safe}}$ .

The newly added no-fly zone

Due to the characteristic of RRT to generate feasible path in a short period of time, it is very suitable for path planning in real-time. The standard form of RRT is suitable for solving planning problem in continuous space. In Wu and Hu (2020), an improved RRT (IRRT) which is suitable for solving path planning problem in discrete space is proposed, whereas it does not consider the UAV’s performance constraints. Therefore, an improved RRT with angle constraints (noted as CIRRT) and a novel evaluation function of path are proposed. Furthermore, the replanned path will return to the offline path quickly to reduce the calculation complexity. And the distance between the replanned path and the newly added no-fly zone should be always no less than $d_{\mathrm{safe}}$ . This strategy will be executed at once if the distance between the UAV and the obstacles will be shorter than $d_{\mathrm{safe}}$ after $T_0$ . The critical collision path point represents as $(P_c)$ .

In the standard RRT, the obtained candidate nodes can be added to the growing-tree only when there is no obstacle between the current node and the candidate node. As shown in Figure 3, this constraint is automatically satisfied by the node expansion definition of the proposed CIRRT. Only the nodes around the current node which are not an obstacle can be selected as the candidate nodes. Meanwhile, the climbing and turning angle constraints of the UAV need to be considered. The unfeasible directions need to be removed to reduce unnecessary node expansion. One of the angle constraints situation is shown in Figure 3. $P_j$ , $P_{j-1}$ and $P_{j+1}$ defined in equations (1) and (2) must satisfy the constraints of turning and climbing angle. The vectors which are composed of the three consecutive path points are expressed as ${\mathit{P}}_{j-1}{\mathit{P}}_j$ and ${\mathit{P}}_j{\mathit{P}}_{j+1}$ .

Figure 3.

Diagram of CIRRT.

For the turning angle, the dot product of the projection vectors of the two vectors on the horizontal plane must be greater than $\sqrt{2}/2$ . For the climbing angle, the dot product of ${\mathit{P}}_j{\mathit{P}}_{j+1}$ and its projection vector must be no less than 0. The red nodes in Figure 4 represent the feasible candidate nodes, and the red arrows indicate the directions that cannot be selected. It is noted that if the distance between $P_j$ and the local goal $(P_{\lg })$ is smaller than $(2\times \sqrt{5}\times \mathrm{\Delta }d)$ , the local goal also must be considered into the angle constraints. And the angle constraints of $P_c$ should be considered when finding the first replanned path point. When the threat of the fixed no-fly zone is successfully evaded by the UAV, it will return to the offline optimized path as soon as possible. The pseudo-code of CIRRT is shown as Algorithm 2.

Figure 4.

The candidate nodes which consider angle constraints.

Algorithm 2: Pseudo-code of CIRRT

1 while $\parallel P_{j+1}-P_{lg}\parallel >\sqrt{3}\times \mathrm{\Delta }d$ do 2  $\mathrm{mark}$ = 0’ 3 while $mark=0$ do 4  for $n=1:\mathrm{length}(T_U)$ do 5    $d(n)=\parallel T_U(n)-T_R\parallel$ do 6  end 7   $p_j=\arg min\left(d\left(n\right)\right)$ 8  if $P_j=T_c$ then 9   Select the feasible candidate nodes by equations (1) and (2) (Consider the angle constraints of the former path point of the $P_c$ , $P_c$ and $P_{j+1})$ 10   else if $\parallel P_j-P_{lg}\parallel <2\times \sqrt{5}\times \mathrm{\Delta }d$ then 11    Select the feasible candidate nodes by equations (1) and (2) (Taking the angle constraints of ${\mathit{P}}_j{\mathit{P}}_{j+1}$ and ${\mathit{P}}_{j+1}{\mathit{P}}_{lg}$ ; the angle constraints of ${\mathit{P}}_{j+1}{\mathit{P}}_{\lg }$ and ${\mathit{P}}_{\lg }{\mathit{P}}_{\lg +1}$ into consideration. $(P_{\lg +1}$ : the latter path point of the $P_{\lg })$ ) 12    else 13     Select the feasible candidate nodes by equations (1) and (2) 14    end 15   end 16  end 17  if $T_{\mathrm{cp}}\ne \Phi$ then 17   mark = 1 18  end 19 end 20 for $s=1:length(T_{cp})$ do 21   $r(s)=\parallel T_{\mathrm{cp}}(s)-T_R\parallel$ 22 end 23  $P_{j+1}=$ arg min $(r(s))$ 24 end

CIRRT: improved RRT with angle constraints.

In Algorithm 2, $T_U$ represents a set of nodes which is composed of all nodes in the growing-tree. $T_R$ is a point randomly selected in the searching space, and $T_{\mathrm{cp}}$ is constituted by the feasible candidate nodes of the current node. $P_g$ is the goal point. arg min represents the value of the variable that minimizes the value of the objective function. So, arg min $(d(n))$ represents n that minimizes d. arg min $(r(s))$ represents s that minimizes r. Therefore, $P_j$ is the node on the growing-tree which is closest to the generated random point. $P_{j+1}$ is the node in the candidate nodes’ set which is closest to the random point. The candidate nodes are selected from the eight directions around the current node. As shown in Figure 3, the distance from these red candidate nodes to the blue random point will be calculated. And the node with the shortest distance will be appended to the growing-tree. If there are no feasible nodes, the random point will be regenerated. And the above process will be repeated until reaching the goal. Finally, feasible paths from the start point to the goal point will be generated. A novel path quality’s evaluation function is proposed to select the optimal path, which is in the form of

$$\begin{array}{c}{J}_1=l_1·\frac{\sum _{j=1}^{N_{\mathrm{maxs}}-1}\parallel P_{j+1}-P_j\parallel }{\sum _{s=1}^M\sum _{j=1}^{N_{\mathrm{maxs}}-1}\parallel P_{j+1}-P_j\parallel }\hfill \\ +l_2·\frac{(N_{\mathrm{maxs}}-2)}{\sum _{s=1}^M(N_{\mathrm{maxs}}-2)}\hfill \\ +l_3·\frac{\sum _{j=2}^{N_{\mathrm{maxs}}-1}{\psi }_{j-1,j,j+1}}{\sum _{s=1}^M\sum _{j=2}^{N_{\mathrm{maxs}}-1}{\psi }_{j-1,j,j+1}}\hfill \\ +l_4·\frac{\sum _{j=3}^{N_{\mathrm{maxs}}-1}|{\varphi }_{j-1,j,j+1}-{\varphi }_{j-2,j-1,j}|}{\sum _{s=1}^M\sum _{j=3}^{N_{\mathrm{maxs}}-1}|{\varphi }_{j-1,j,j+1}-{\varphi }_{j-2,j-1,j}|}\hfill \end{array}$$

(20)

The function contains parts of path length, number of path points, and angle, which aims to reduce the energy consumption. The measurement of each index in the function is different; therefore, the normalization operation is carried out to avoid unfair advantage over others due to the larger value of one of them. M denotes the number of generated paths, and $N_{\mathrm{maxs}}$ denotes the path point number of the $s\mathrm{th}$ path. In addition, $l_1$ , $l_2$ , $l_3$ , and $l_4$ indicate weights of different index.

The obstacle with unknown trajectory

Inspired by Xi and Zhang (2003), a rolling optimization principle combined with three-dimensional rolling windows method is proposed to avoid the threat of the unknown dynamic obstacles. The rolling optimization principle has been studied in 2D environment. And this principle is extended to three-dimensional environment for the UAV in this section. The sensors’ detection area of the UAV is known as the optimization window, and it will move forward with each $T_0$ . The environment information in the window is constantly updated. Then, the information is utilized to judge whether there is a collision. If there is no collision, the UAV will fly along the original path. Otherwise, it will plan a local optimal path. The UAV will fly along the new path until the next time of online planning. The sensors will be always detecting the global goal in the current optimization window. If the global goal is detected, it indicates that it is the last time of online planning.

At first, the model of obstacles with unknown trajectory is proposed. The path of them during a rolling period $(T_r)$ is difficult to predict due to the uncertainty of its moving direction. Therefore, the possible influence range of the obstacle must be expanded to predict its moving range during $T_r$ . The possible influence range of obstacle is defined as a sphere which sets its current point $(O=(x_k,y_k,z_k))$ as the center point, and the threat radius of the obstacle is defined as $R_k$ . The speed of obstacle is in the range of $[0,v_{\mathrm{ok}}]$ , in which $k=1,2,...,v$ indicates the $k\mathrm{th}$ unknown dynamic obstacle detected by the sensors. The obstacle’s possible influence range during $T_r$ is defined as

$$(x-x_k{)}^2+(y-y_k{)}^2+(z-z_k{)}^2={(R_k+v_{\mathrm{ok}}\times T_r)}^2$$

(21)

Then, it is necessary to predict whether the expansion region of the obstacle intersects with the global path or not. If not, there will be no collision during the next rolling period. On the contrary, a further prediction is required. Then, it is necessary to predict the position of the UAV and the obstacle during the next $T_0$ and the next two $T_0$ . The possible situations that the UAV will encounter are shown in Figure 5. The fan-shaped area represents partial detection scope of the UAV. $P(T)$ and $V_a$ denote the current position and speed of the UAV; $P(T+T_0)$ and $D{O}_k(T+T_0)$ are the maximum motion position of the UAV and the obstacle during the next $T_0$ , respectively. $D{O}_k(T+2T_0)$ refers to the maximum motion range of the obstacle during the next two $T_0$ . (1) As shown in Figure 5(a), it indicates that there will be an intersection between the motion range of the UAV and the unknown dynamic obstacle during the next $T_0$ , which is called dangerous state. (2) In the quasi dangerous state in Figure 5(b), the collision will not happen during the next $T_0$ , but it will happen during the next two $T_0$ . If the UAV continues to fly along the offline path, the dangerous state will occur after the next $T_0$ . In order to avoid the influence of these two situations, obstacle avoidance measures must be taken immediately at the current time. The corresponding expansion time is $T_0$ . (3) In the potential dangerous state shown in Figure 5(c), the collision will not happen during the next two $T_0$ . If obstacle avoidance measure is taken immediately, the effectiveness of the replanned local path for some time to come cannot be guaranteed. (4) In Figure 5(d), it indicates that the moving direction of the obstacle is opposite to that of the UAV. Thus, the UAV will continue to flight along the offline optimized path without collision.

Figure 5.

The possible position of the UAV and the obstacles with unknown trajectory (a) dangerous state, (b) quasi dangerous state,(c) potentially hazardous state and (d) secure state.

The definition of the local goal during the rolling optimization process is proposed. First, the intersection part of the obstacle’s expansion range and the current rolling window will be recognized as the static obstacle in the grid-based environment. Then, the global path point within the rolling window which is closest to the intersection point will generally be selected as the local goal point. If the local goal is located on the maximum arrival range of the obstacle, the latter global path point will be set as the new local goal. After determining the local goal, the feasible candidate nodes are selected similarly to the situation of fixed no-fly zones. Different from the situation of fixed no-fly zones, only one optimal path will be generated. A novel cost function as shown in equation (22) is proposed to select the optimal path point

$$J_3=d_{\mathrm{jc}}+\Vert P_{j+1}-P_j\Vert +\Vert P_{\lg }-P_{j+1}\Vert$$

(22)

$d_{\mathrm{jc}}$ is the actual distance from the local starting point $P_c$ to the current path point $P_j$ . Point $P_{j+1}$ with the lowest function value will be selected as the next path point of the UAV. Finally, the UAV will perform the predicting, avoiding and flying process constantly to ensure itself security when advancing along the replanned path.

Simulation results

Benchmark function tests and comparisons

In order to verify the performance of DGSHGS, it is used to find the optimal values of the standard optimization functions which is commonly utilized in Rashedi et al. (2009). There are 23 benchmark functions, in which $f_1-f_7$ are seven uni-modal functions, $f_8-f_{13}$ are six multi-modal functions, and $f_{14}-f_{23}$ are multi-modal functions with fixed dimension. The dimension of uni-modal functions are set as 50 and the specific information of all benchmark functions are shown in Table 1. Some existing algorithms and their variations are compared with the proposed DGSHGS, such as HGS (Yang et al., 2021), COA (Pierezan et al., 2019), FPA (Yang, 2012), Grey Wolf Optimization (GWO) (Mirjalili et al., 2014), Selective Opposition Grey Wolf Optimization (SOGWO) (Dhargupta et al., 2020), Particle Swarm Optimization Based Grey Wolf Optimization (PSOGWO) (Gul et al., 2021), Multi-Strategy Ensemble Grey Wolf Optimizer (MEGWO) (Tu et al., 2019), and Grey Wolf Optimization Algorithm Based on Adaptive Normal Cloud Model (CGWO) (Zhang et al., 2020). In order to realize a fair comparison, the total times of utilizing the evaluation function during the running process of each algorithm are set as 6000. The population’s size of all algorithms is $30$ . In order to reduce the random factors’ influence, the times of independent tests for each algorithm is $25$ . The parameters’ detailed definitions in all algorithm are shown in Table 2.

Table 1.

Description of the 23 benchmark functions.

ID	Function equation	Dim	Range	${\mathrm{f}}_{\min }$
F1	$f_1(x)=\sum _{i=1}^n{x}_i^2$	50	[–100,100]	0
F2	$f_2(x)=\sum _{i=1}^n\left|x_i\right|+\underset{i=1}{\overset{n}{\Pi }}\left|x_i\right|$	50	[–10,10]	0
F3	$f_3(x)=\sum _{i=1}^n{\left(\sum _{j=1}^i{x}_j\right)}^2$	50	[–100,100]	0
F4	$f_4(x)=\underset{i}{max}\{\left|x_i\right|,1\le i\le n\}$	50	[–100,100]	0
F5	$f_5(x)=\sum _{i=1}^{n-1}[100{(x_{i+1}-x_i^2)}^2+{(x_i-1)}^2]$	50	[–30,30]	0
F6	$f_6(x)=\sum _{i=1}^n{\left([x_i+0.5]\right)}^2$	50	[–100,100]	0
F7	$f_7(x)=\sum _{i=1}^n{\mathrm{ix}}_i^4+\mathrm{random}[0,1)$	50	[–1.28,1.28]	0
F8	$f_8(x)=\sum _{i=1}^n-x_i\sin \left(\sqrt{\left|x_i\right|}\right)$	50	[–500,500]	–418.982
F9	$f_9(x)=\sum _{i=1}^n[x_i^2-10\cos \left(2\pi x_i\right)+10]$	50	[–5.12,5.12]	0
F10	$f_{10}(x)=-20\exp \{-0.2\sqrt{\frac{1}{n}\sum _{i=1}^n{x}_i}\}-\exp \left\{\frac{1}{n}\sum _{i=1}^n\cos \left(2\pi x_i\right)\right\}+20+e$	50	[–32,32]	0
F11	$f_{11}(x)=\frac{1}{4000}\sum _{i=1}^n{x}_i^2-\underset{i=1}{\overset{n}{\Pi }}\cos \left(\frac{x_i}{\sqrt{i}}\right)+1$	50	[–600,600]	0
F12	$f_{12}(x)=\frac{\pi }{n}\{10\sin \left(a{y}_1\right)+\sum _{i=1}^{n-1}{(y_i-1)}^2[1+10\stackrel{2}{\sin }\left(\pi y_{i+1}\right)]+{(y_n-1)}^2+\sum _{i=1}^n\mu (x_i,10,100,4)\}$	50	[–50,50]	0
F13	$f_{13}(x)=0.1\{\stackrel{2}{\sin }\left(3\pi x_i\right)+\sum _{i=1}^n{(x_i-1)}^2[1+\stackrel{2}{\sin }(3\pi x_i+1)]+{(x_n-1)}^2[1+\stackrel{2}{\sin }\left(2\pi x_n\right)]+\sum _{i=1}^n\mu (x_i,5,100,4)$	50	[–50,50]	0
F14	$f_{14}(x)={(\frac{1}{500}+\sum _{j=1}^{25}\frac{1}{j+\sum _{i=1}^2{(x_i+a_{\mathrm{ij}})}^6})}^{-1}$	2	[–65,65]	1
F15	$f_{15}(x)=\sum _{i=1}^{11}{[a_i-\frac{x_1(b_i^2+b_i{x}_2)}{b_i^2+b_i{x}_3+x_4}]}^2$	4	[–5,5]	0.00030
F16	$f_{16}(x)=4x_1^2-2·1x_1^4+\frac{1}{2}{x}_1^6+x_1{x}_2-4x_2^2+4x_2^4$	2	[–5,5]	–1.0316
F17	$f_{17}(x)={(x_2-\frac{5.1}{4{\pi }^2}{x}_1^2+\frac{5}{\pi }{x}_1-6)}^2+10(1-\frac{1}{8\pi })\cos x_1+10$	2	[–5,5]	0.398
F18	$\begin{array}{c}{f}_{18}(x)=[1+{(x_1+x_2+1)}^2(19-14x_1+3x_1^2-14x_2+6x_1{x}_2+3x_2^2)]\times \hfill \\ \multicolumn{1}{c}{[30+{(2x_1-3x_2)}^2\times (18-32x_1+12x_1^2+48x_2-36x_1{x}_2+27x_2^2)]}\end{array}$	2	[–2,2]	3
F19	$f_{19}(x)=-\sum _{i=1}^4{c}_i\exp (-\sum _{j=1}^3{a}_{\mathrm{ij}}{(x_j-p_{\mathrm{ij}})}^2)$	3	[1,3]	–3.86
F20	$f_{20}(x)=-\sum _{i=1}^4{c}_i\exp (-\sum _{j=1}^6{a}_{\mathrm{ij}}{(x_j-p_{\mathrm{ij}})}^2)$	6	[0,1]	–3.32
F21	$f_{21}(x)=-\sum _{i=1}^5{[(X-a_i){(X-a_i)}^T+c_i]}^{-1}$	4	[0,10]	–10.1532
F22	$f_{22}(x)=-\sum _{i=1}^7{[(X-a_i){(X-a_i)}^T+c_i]}^{-1}$	4	[0,10]	–10.4028
F23	$f_{23}(x)=-\sum _{i=1}^{10}{[(X-a_i){(X-a_i)}^T+c_i]}^{-1}$	4	[0,10]	–10.5363

Table 2.

Parametric details of the compared algorithms.

Algorithm	Parameter settings
FPA	$\mathrm{UN}(0,{\sigma }^2);\mathrm{VN}(0,1);p=0.5$
GWO	$a=[2,0]$
CGWO	$\mathrm{nfeva}{l}_{max}=60,000;a=[2,0]$
MEGWO	$\mathrm{nfeva}{l}_{max}=60,000;a=[2,0]$
SOGWO	$a=[2,0],\mathrm{threshold}=a$
PSOGWO	$c_1=c_2=c_3=[2,0];{\omega }_{\max }=1;{\omega }_{\min }=0.5$
HGS	$u=0.03;\mathrm{LH}=100$
COA	$\mathrm{nfeva}{l}_{max}=60,000,n_p=5,n_c=6$

The comparison results on benchmark functions are shown in Tables 3 and 4. In Tables 3 and 4, the optimal value of the best, the average (Mean) and the standard deviation (std) are taken in bold. The mean value and the standard deviation are inversely proportional to the optimization ability and the stability of the algorithm, respectively. It can be seen from the related numerical results in Tables 3 and 4 that DGSHGS has better performance on 13 functions, similar performance on 7 functions and worse performance on 3 functions compared with the traditional HGS from the aspect of the value of Mean. Among the seven functions with similar performance, DGSHGS and HGS obtain the optimal values on six functions. Compared with HGS, DGSHGS obtains better performance on 13 functions, similar performance on 9 functions and worse performance on 1 function in terms of the value of Std. It is noted that DGSHGS and HGS obtain the optimal values on the all 9 functions with similar performance. Compared with the SOGWO, DGSHGS achieves the better value on the 69.6% of functions in terms of the value of Mean and the better standard deviation value on the 91.3% of functions. DGSHGS achieves the better mean value on the 73.9% of functions and the better standard deviation value on the 91.3% of functions compared with CGWO. Compared with GWO, COA, and PSOGWO, the mean and standard deviation value of DGSHGS on most functions are better than them. The DGSHGS is also superior to FPA and MEGWO. Meanwhile, it is noted that more than half of the functions with similar performance in terms of the value of Mean achieve the optimal function value by them. According to the simulation results and analysis above, it can be concluded that the performances of DGSHGS on the benchmark functions are superior to the other eight algorithms.

Table 3.

Results of DGSHGS, FPA, GWO, COA, and HGS on benchmark functions.

Function	FPA		GWO		COA		HGS		DGSHGS
	Mean	Std	Mean	Std	Mean	Std	Mean	Std	Mean	Std
F1	7.1745E–01	2.5751E–01	1.7197E–91	1.9929E–91	3.8084E–01	1.4131E–01	0.0000E+00	0.0000E+00	0.0000E+00	0.0000E+00
F2	7.1964E–01	9.6921E–02	5.8455E–54	7.0212E–54	1.6217E–02	7.7824E–03	8.1861E–236	0.0000E+00	0.0000E+00	0.0000E+00
F3	1.9667E+04	4.8286E+03	3.5060E–20	3.3445E–20	1.3670E+04	3.6373E+03	0.0000E+00	0.0000E+00	8.6239E–280	0.0000E+00
F4	1.5111E+01	2.0653E+00	6.3231E–19	1.0680E–18	2.7665E+01	6.4622E+00	1.3219E–201	0.0000E+00	0.0000E+00	0.0000E+00
F5	2.3110E+02	2.9401E+01	4.7077E+01	6.1743E–01	4.1713E+02	1.8081E+02	4.7006E+01	1.3113E+00	4.4693E+01	1.3044E–01
F6	5.2046E–01	1.3524E–01	1.9176E+00	3.8094E–01	1.7635E–01	1.1323E–01	1.8159E+00	1.6159E–01	4.1424E–05	1.5467E–05
F7	2.0785E–01	7.7283E–02	4.8195E–04	1.8434E–04	2.1799E–01	5.9731E–02	3.5571E–04	3.9118E–04	1.0758E–04	4.9251E–05
F8	–1.2357E+04	3.3670E+02	–9.1625E+03	1.0813E+03	–2.0829E+04	8.1650E+01	–1.6859E+04	5.4546E+02	–2.0909E+04	7.8955E+01
F9	1.6387E+02	8.8230E+00	0.0000E+00	0.0000E+00	1.6471E+01	3.0620E+00	0.0000E+00	0.0000E+00	0.0000E+00	0.0000E+00
F10	4.9915E+00	2.6985E+00	1.6520E–14	3.1776E–15	2.9572E+00	5.3982E–01	8.8818E–16	0.0000E+00	8.8818E–16	0.0000E+00
F11	5.1916E–01	1.5158E–01	3.4388E–04	1.7194E–03	1.8622E–01	1.2341E–01	0.0000E+00	0.0000E+00	0.0000E+00	0.0000E+00
F12	4.1503E+00	9.2350E–01	8.0233E–02	3.0672E–02	4.0751E–01	4.4978E–01	2.2262E–02	1.0278E–02	7.9439E–07	7.1421E–07
F13	2.3927E+01	1.0957E+01	1.7498E+00	3.7165E–01	7.2171E+00	6.5220E+00	7.5876E–01	2.6919E–01	3.0103E–02	3.6803E–02
F14	9.9800E–01	1.0279E–16	4.5280E+00	3.9828E+00	9.9800E–01	0.0000E+00	1.6490E+00	2.5214E+00	1.4614E+00	7.3762E–01
F15	4.1420E–04	1.0104E–04	2.9816E–03	7.0570E–03	4.0021E–04	2.5847E–04	3.2080E–03	6.9746E–03	7.2413E–04	2.6145E–04
F16	–1.0316E+00	0.0000E+00	–1.0316E+00	1.0273E–09	–1.0316E+00	0.0000E+00	–1.0316E+00	0.0000E+00	–1.0316E+00	0.0000E+00
F17	3.9789E–01	0.0000E+00	3.9789E–01	4.4047E–08	3.9789E–01	0.0000E+00	3.9789E–01	0.0000E+00	3.9789E–01	0.0000E+00
F18	3.0000E+00	6.5455E–16	5.7000E+00	1.4789E+01	3.0000E+00	1.8254E–15	3.0000E+00	3.0367E–15	3.0000E+00	1.5973E–14
F19	–3.8628E+00	2.7101E–15	–3.8614E+00	2.6991E–03	–3.8628E+00	2.7101E–15	–3.8612E+00	3.2065E–03	–3.8628E+00	2.2662E–15
F20	–3.3220E+00	4.8793E–12	–3.2778E+00	5.9029E–02	–3.2744E+00	5.9241E–02	–3.2432E+00	1.0323E–01	–3.2784E+00	5.8273E–02
F21	–1.01304E+01	1.2466E–01	–9.8148E+00	1.2876E+00	–9.9493E+00	9.2244E–01	–9.2310E+00	2.4230E+00	–9.9833E+00	9.3076E–01
F22	–1.04029E+01	1.8514E–04	–1.04028E+01	8.0283E–05	–1.04029E+01	1.1629E–15	–8.2619E+00	3.9953E+00	–1.04029E+01	1.8467E–05
F23	–1.05364E+01	2.6360E–06	–1.05363E+01	3.6840E–05	–1.0090E+01	1.7302E+00	–8.4296E+00	3.6373E+00	–1.05364E+01	4.9634E–05
Lose	12	15	16	21	15	16	13	13	–	–
Win	5	6	4	1	3	5	3	1	–	–
Draw	6	2	3	1	5	2	7	9	–	–

Table 4.

Results of DGSHGS, SOGWO, PSOGWO, MEGWO, and CGWO on benchmark functions.

Function	SOGWO		PSOGWO		MEGWO		CGWO		DGSHGS
	Mean	Std	Mean	Std	Mean	Std	Mean	Std	Mean	Std
F1	2.2619E–91	1.8212E–91	1.3155E–33	2.2079E–33	7.8724E–44	1.1459E–43	2.6922E–91	2.8747E–91	0.0000E+00	0.0000E+00
F2	4.2539E–54	2.9258E–54	8.9738E–07	1.5543E–06	8.5346E–31	5.3917E–31	1.1592E–52	8.1200E–53	0.0000E+00	0.0000E+00
F3	3.0178E–15	5.0605E–15	9.9753E+02	1.7278E+03	2.0017E+03	1.1551E+03	6.7526E–11	7.8908E–11	8.6239E–280	0.0000E+00
F4	1.7638E–19	1.9903E–19	1.3408E+01	1.8707E+01	2.7405E–07	1.6710E–07	4.0662E–20	3.7528E–20	0.0000E+00	0.0000E+00
F5	4.7043E+01	8.2088E–01	4.4791E+07	7.7530E+07	4.5021E+01	4.8207E–01	4.6550E+01	1.1038E+00	4.4693E+01	1.3044E–01
F6	2.3286E+00	6.3110E–01	4.3861E–01	3.5175E–01	7.3345E–20	7.3640E–20	1.4063E–02	7.6224E–04	4.1424E–05	1.5467E–05
F7	5.5591E–04	2.6138E–04	1.0309E+00	1.6619E+00	4.6844E–03	1.6695E–03	4.1426E–04	2.0091E–04	1.0758E–04	4.9251E–05
F8	–9.6437E+03	1.1459E+03	–9.4738E+03	3.4257E+03	–2.0949E+04	3.8348E–12	–7.9692E+03	1.2402E+03	–2.0909E+04	7.8955E+01
F9	0.0000E+00	0.0000E+00	1.3708E+02	1.8554E+02	0.0000E+00	0.0000E+00	0.0000E+00	0.0000E+00	0.0000E+00	0.0000E+00
F10	1.5810E–14	1.5888E–15	3.4671E+00	7.7236E+00	6.5725E–15	1.9459E–15	1.5099E–14	0.0000E+00	8.8818E–16	0.0000E+00
F11	2.0305E–03	4.7618E–03	1.9443E+01	9.5147E+01	0.0000E+00	0.0000E+00	1.7872E–03	5.2066E–03	0.0000E+00	0.0000E+00
F12	7.7612E–02	2.8976E–02	1.4201E+05	5.4999E+05	4.8600E–14	8.7869E–14	2.9912E–04	3.6329E–05	7.9439E–07	7.1421E–07
F13	1.7188E+00	3.0290E–01	1.2011E+08	3.3538E+08	2.5541E–13	4.8514E–13	3.1589E–02	7.8223E–02	3.0103E–02	3.6803E–02
F14	4.1158E+00	5.0866E+00	5.7389E+00	5.6464E+00	3.3325E+00	4.8329E+00	3.6073E+00	4.1454E+00	1.4614E+00	7.3762E+00
F15	2.9894E–03	7.0539E–03	5.6557E–03	9.1803E–03	4.1737E–04	3.0370E–04	1.6447E–03	5.1784E–03	7.2413E–04	2.6145E–04
F16	–1.0316E+00	2.5630E–09	–1.0316E+00	2.6936E–05	–1.0316E+00	9.8203E–15	–1.0316E+00	1.8592E–07	–1.0316E+00	0.0000E+00
F17	3.9789E–01	4.7174E–08	3.9789E–01	1.0055E–05	3.9789E–01	4.5475E–14	3.9789E–01	3.9110E–06	3.9789E–01	0.0000E+00
F18	5.7000E+00	1.4789E+01	3.0136E+00	4.5988E–02	3.0000E+00	5.8185E–08	3.0000E+00	3.2368E–06	3.0000E+00	1.5973E–14
F19	–3.8615E+00	2.7593E–03	–3.8604E+00	3.6756E–03	–3.8628E+00	2.6823E–15	–3.8627E+00	1.1224E–04	–3.8628E+00	2.2662E–15
F20	–3.3220E+00	7.3972E–02	–3.1035E+00	3.2083E–01	–3.2744E+00	2.3707E–09	–3.2364E+00	6.4291E–02	–3.2784E+00	5.8273E–02
F21	–9.3080E+00	1.9220E+00	–7.2683E+00	3.0228E+00	–6.5602E+00	3.3248E+00	–8.9531E+00	2.1657E+00	–9.9833E+00	9.3076E–01
F22	–9.7435E+00	1.8647E+00	–8.6112E+00	3.3265E+00	–7.4636E+00	3.2361E+00	–9.7166E+00	1.8550E+00	–1.04029E+01	1.8467E–05
F23	–1.05363E+01	5.7077E–05	–7.2592E+00	3.2209E+00	–6.9029E+00	3.6152E+00	–9.6078E+00	2.4237E+00	–1.05364E+01	4.9634E–05
Lose	16	21	19	23	13	16	17	21	–	–
Win	4	1	2	0	4	5	2	0	–	–
Draw	3	1	2	0	6	2	4	2	–	–

Offline path planning results and comparisons

At first, two $40\times 30\times 40$ three-dimensional grid-based environments with some static obstacles are established. The first and the second scenarios have the different environment, and the start point and the goal point of them are different. In the first scenario, the positions of the start and the goal point are located on $(3,3,15)$ and $(28,27,28)$ , respectively. In the second scenario, the coordinates of them are on $(6,7,12)$ and $(28,26,22)$ , respectively. The value of ${\psi }_{\max }$ is set as $45°$ , and the value of ${\varphi }_{\max }$ is set as $90°$ . In order to keep fair, the times of utilizing the evaluation function of each algorithm are set as $20,000$ . $k_1=0.6,k_2=0.1,k_3=0.15,k_4=0.15,k_5=1.$ Figure 6 shows the comparison results between DGSHGS, GWO, PSOGWO, COA, and HGS when being applied into the problem of three-dimensional offline path planning.

Figure 6.

3D path planning results comparison: (a) scenario 1 and (b) scenario 2.

As shown in Table 5, the path generated by DGSHGS has fewer number of turns and shorter length than other four compared algorithms. Meanwhile, the climbing and turning angles are also smaller than other four algorithms. As shown in Figure 6, the generated paths keep a certain distance from the obstacles, and the angle constraints have been satisfied. The convergence curves in Figure 7 denote that the convergence value of DGSHGS is always the lowest at the same iteration during the evolutionary process. Moreover, the number of iterations when reaching the same evaluation function value is lesser than the others. It can be concluded that DGSHGS has excellent global exploration ability and high efficiency during iterations course. Thus, the above results verify the effectiveness and superiority of DGSHGS for the UAV path planning in low-altitude complex urban environment.

Table 5.

Simulation values of the offline planning.

		The objective function value	Path length ( $\mathrm{\Delta }d$ )	Turns	Turning angle (rad)	Climbing angle (rad)
Scenario 1	GWO	65.6000	45.2876	17	5.9088	6.8802
	PSOGWO	57.4467	46.1783	11	3.3903	4.1732
	COA	56.7700	44.6438	13	4.0689	4.1135
	HGS	55.1240	43.1545	10	2.4981	3.7541
	DGSHGS	50.9834	41.3550	4	1.4289	1.1668
Scenario 2	GWO	55.4219	36.0274	15	2.8198	7.4806
	PSOGWO	60.7245	40.0302	10	2.8929	8.6279
	COA	53.7465	35.1078	13	1.8925	3.0554
	HGS	54.2367	36.0766	5	1.7507	2.2373
	DGSHGS	52.0446	34.9581	5	1.4289	1.7908

Figure 7.

Convergence curves of DGSHGS, GWO, PSOGWO, COA and HGS for 3D path planning: (a) scenario 1 and (b) scenario 2.

Online path planning results and comparisons

The serial number coordinates of the newly added fixed no-fly area is set as $(10,11,27)$ , and the initial serial number coordinate of the unknown dynamic obstacles are set as $(23,15,28)$ and $(26,22,28)$ , respectively. The initial serial number coordinate of the known dynamic obstacles are set as $(16,26,28)$ and $(27,12,28)$ . And the center of the unknown and known dynamic obstacles are represented as cyan and purple-red dot, respectively. $l_1$  = 0.1, $l_2$  = 0.3, $l_3$  = 0.3, and $l_4$  = 0.3 The related settings of parameters are shown in Table 6, and the comparison results of IRRT (Wu and Hu, 2020) and CIRRT-based online planning strategies to avoid the influence of fixed no-fly zone are shown in Figure 8. A local path is generated between the critical collision path point and the global goal in Figure 8(a). The maximum turning angle of the path which is planned by IRRT is greater than ${45}^{°}$ . It is hard for the UAV to perform. And the generated local path deviates from the offline optimized path. As shown in Figure 8(b), the UAV abandons the original path and replans a local optimal path by proposed CIRRT when confronting with the threat of fixed no-fly zone. The online planning algorithm starts immediately if the UAV skips over the critical collision path point at the next sampling time. After avoiding the threat of the obstacle, the UAV returns to the offline optimized path quickly. That reduces the computation time of the algorithm obviously.

Table 6.

Parametric details of the online path planning algorithms.

Parameter	Value	Parameter	Value
$T_0$	$0.1$ $\mathrm{s}$	$\mathrm{\Delta }d$	$10$ $\mathrm{m}$
$d_{\mathrm{safe}}$	$20$ $\mathrm{m}$	R	$50$ $\mathrm{m}$
$V_H$	$3$ $\mathrm{m}/\mathrm{s}$	$R_k$	$20$ $\mathrm{m}$
$V_M$	$6$ $\mathrm{m}/\mathrm{s}$	$V_{\mathrm{ok}}$	$1.7$ $\mathrm{m}/\mathrm{s}$
$V_L$	$9$ $\mathrm{m}/\mathrm{s}$	M	$5$

Figure 8.

Planning results comparison of IRRT and CIRRT based online planning strategies: (a) IRRT and (b) CIRRT.

Therefore, the effectiveness of CIRRT based online planning strategy to avoid the newly added no-fly zone is verified. When encountering the newly added no-fly zone, the UAV’s flight speed remains constant. Therefore, the length of the path can be considered as an evaluation index of the energy consumption. The turning number and the smoothness of the UAV’s path also have great influence on energy consumption. As shown in Table 7, the performance of the proposed CIRRT based online planning strategy is better than that based on IRRT in terms of the path planning time, path length and angles when working in the same environment. The planning time and path angles are improved by 11.3% and 56.8%, respectively. And the length of the path planned by proposed CIRRT-based online planning strategy is also smaller than that based on IRRT. In conclusion, the energy consumption of the path generated by the proposed CIRRT is better than that generated by the IRRT.

Table 7.

Comparison of IRRT and CIRRT based online strategies.

	IRRT	CIRRT
Planning time	$3.2174\times {10}^1$ $\mathrm{ms}$	$2.8551\times {10}^1$ $\mathrm{ms}$
Path length	$3.1780\times {10}^2$ $\mathrm{m}$	$3.0509\times {10}^2$ $\mathrm{m}$
Climbing angle	$423.{1727}^{°}$	$215.{2644}^{°}$
Turning angle	${180}^{°}$	${45}^{°}$

IRRT: improved RRT; CIRRT: improved RRT with angle constraints.

After avoiding the threat of no-fly zone, the first known dynamic obstacle is being detected by sensors. As shown in Table 8, the position coordinates of this obstacle which is detected by sensors have been recorded. Meanwhile, the maximum motion velocity $(v_{o1})$ of the first unknown dynamic obstacle is calculated as $1.7$ m/s. The center of the unknown dynamic obstacles are represented as cyan dots. In Figure 9(a), the path replanned by the proposed 3D rolling windows method successfully avoids this unknown dynamic obstacle and satisfies the UAV’s performance constraints. And the distance between the UAV and this unknown dynamic obstacle is always bigger than $d_{\mathrm{safe}}$ during the execution of the proposed rolling windows method. When encountering the first known dynamic obstacle, the UAV avoids it by changing its speed. The speed of the known dynamic obstacle is $1.5$ m/s, and the speed of the UAV is set as $6$ m/s first. And the center of the known dynamic obstacle is represented as a purple-red dot. The velocity change curve of the UAV is shown in Figure 11. When the distance between the known obstacle and the UAV is less than $d_{\mathrm{safe}}$ after a short period of time, the changing speed strategy will be activated. Whereas, there are no feasible speed for the UAV to avoid the obstacle successfully. Therefore, the UAV begins suddenly to slow down when t = 49.4 s, and the deceleration is set as 6 m/ ${\mathrm{s}}^2$ . As shown in Figure 9(b), the UAV hovers in the original position when t = 50.4 s and waits for the next sampling time to judge again. This process continues until $t=54.2$ s, and the UAV adjusts the speed to 3 m/s. After that, the UAV remains the current speed (3 m/s) until the next speed changing. The position coordinates of the second unknown obstacle is shown in Table 9. The maximum speed of it is calculated as $1.6$ m/s. In Figure 10(a), the replanned path successfully avoids this unknown dynamic obstacle and satisfies the related constraints. The distance between the UAV and the second unknown dynamic obstacle is always bigger than $d_{\mathrm{safe}}$ after the implementation of proposed method. Finally, the UAV is influenced by the second known dynamic obstacle. The speed of it is set as $1.5$ m/s. As shown in Figure 11(a), the speed of the UAV is $3$ m/s through the last speed changing. The UAV begins suddenly to slow down when $t=66.4$ s. As shown in Figure 10(b), the UAV hovers in the original position when $t=66.9$ s and waits for the next sampling time to judge again. The above process continues until $t=70.3$ s, and the UAV adjusts the speed to 3 m/s. As shown in Figure 11(b), the UAV successfully avoids all kinds of obstacles and arrives at the goal. It can be concluded that the proposed bi-level structure algorithm can solve the dynamic planning problem with fine performances for the UAV in low-altitude complex urban environment effectively.

Table 8.

Position coordinates of the first unknown dynamic obstacle which are detected by the airborne sensors.

Time	Position coordinate	Time	Position coordinate	Time	Position coordinate	Time	Position coordinate
t = 35.7 s	(16.9140,15,28)	t = 35.8 s	(16.9039,15,28)	t = 35.9 s	(16.8943,15,28)	t = 36.0 s	(16.8781,15,28)
t = 36.1 s	(16.8650,15,28)	t = 36.2 s	(16.8571,15,28)	t = 36.3 s	(16.8438,15,28)	t = 36.4 s	(16.8349,15,28)
t = 36.5 s	(16.8201,15,28)	t = 36.6 s	(16.8044,15,28)	t = 36.7 s	(16.7941,15,28)	t = 36.8 s	(16.7858,15,28)
t = 36.9 s	(16.7711,15,28)	t = 37.0 s	(16.7609,15,28)	t = 37.1 s	(16.7541,15,28)	t = 37.2 s	(16.7424,15,28)
t = 37.3 s	(16.7347,15,28)	t = 37.4 s	(16.7222,15,28)	t = 37.5 s	(16.7111,15,28)	t = 37.6 s	(16.6966,15,28)
t = 37.7 s	(16.6848,15,28)	t = 37.8 s	(16.6739,15,28)	t = 37.9 s	(16.6663,15,28)	t = 38.0 s	(16.6572,15,28)
t = 38.1 s	(16.6447,15,28)	t = 38.2 s	(16.6336,15,28)	t = 38.3 s	(16.6176,15,28)	t = 38.4 s	(16.6101,15,28)
t = 38.5 s	(16.5986,15,28)	t = 38.6 s	(16.5862,15,28)	t = 38.7 s	(16.5779,15,28)	t = 38.8 s	(16.5655,15,28)
t = 38.9 s	(16.5449,15,28)	t = 39.0 s	(16.5409,15,28)	t = 39.1 s	(16.5324,15,28)	t = 39.2 s	(16.5248,15,28)
t = 39.3 s	(16.5112,15,28)	t = 39.4 s	(16.5040,15,28)	t = 39.5 s	(16.4941,15,28)

Figure 9.

Online path planning in the case of the first unknown and known dynamic obstacle: (a) first unknown dynamic obstacle and (b) first known dynamic obstacle.

Figure 10.

Online path planning in the case of the second unknown and known dynamic obstacle: (a) second unknown dynamic obstacle and (b) second known dynamic obstacle.

Table 9.

Position coordinates of the second unknown dynamic obstacle which are detected by the airborne sensors.

Time	Position coordinate	Time	Position coordinate	Time	Position coordinate	Time	Position coordinate
t = 59.2 s	(18.9845,22,28)	t = 59.3 s	(18.9762,22,28)	t = 59.4 s	(18.9612,22,28)	t = 59.5 s	(18.9494,22,28)
t = 59.6 s	(18.9418,22,28)	t = 59.7 s	(18.9282,22,28)	t = 59.8 s	(18.9206,22,28)	t = 59.9 s	(18.9057,22,28)
t = 60.0 s	(18.8924,22,28)	t = 60.1 s	(18.8812,22,28)	t = 60.2 s	(18.8728,22,28)	t = 60.3 s	(18.8596,22,28)
t = 60.4 s	(18.8443,22,28)	t = 60.5 s	(18.8362,22,28)	t = 60.6 s	(18.8217,22,28)	t = 60.7 s	(18.8125,22,28)
t = 60.8 s	(18.7986,22,28)	t = 60.9 s	(18.7852,22,28)	t = 61.0 s	(18.7762,22,28)	t = 61.1 s	(18.7638,22,28)
t = 61.2 s	(18.7521,22,28)	t = 61.3 s	(18.7390,22,28)	t = 61.4 s	(18.7301,22,28)	t = 61.5 s	(18.7153,22,28)
t = 61.6 s	(18.7052,22,28)	t = 61.7 s	(18.6974,22,28)	t = 61.8 s	(18.6843,22,28)	t = 61.9 s	(18.6747,22,28)
t = 62.0 s	(18.6606,22,28)	t = 62.1 s	(18.6442,22,28)	t = 62.2 s	(18.6300,22,28)	t = 62.3 s	(18.6175,22,28)
t = 62.4 s	(18.6016,22,28)	t = 62.5 s	(18.5871,22,28)	t = 62.6 s	(18.5753,22,28)	t = 62.7 s	(18.5625,22,28)

Figure 11.

The velocity curve and path of the UAV: (a) velocity curve of the UAV and (b) the UAV successfully reaches the goal point.

Conclusion

In order to solve dynamic path planning problem for the UAV in low-altitude complex urban environment, this work proposes a novel bi-level structure algorithm which combines the offline and online planning. Before path planning of the UAV, a three-dimensional grid-based environment and the related flight constraints are established. The offline planning is a complex global optimization problem with multiple constraints, and it focus on the optimality and stability of the path. Considering the effectiveness of SIOAs to solve the complex global optimization problem, an improved DGSHGS is developed to generate the offline optimized path. In order to improve the efficiency of solving dynamic obstacle avoidance problem, the dynamic obstacles during the mission are classified into three categories according to their different characteristics. Then, on the basis of the offline path, the UAV’s performance constraints and the sensors’ information, three online path planning strategies: CIRRT, a changing speed strategy, and an novel three-dimensional rolling windows method are introduced to update the speed or path of the UAV in real-time. Simulation results shown that DGSHGS outperforms many compared algorithms on the $23$ well-known benchmark functions. The offline optimized path generated by DGSHGS has fewer number of turns and shorter length than other compared algorithms. The online path planning strategies can effectively avoid different kinds of dynamic obstacles under the constraints of offline path. More specifically, the planning time and angles of local path generated by CIRRT to avoid the influence of no-fly zone are improved by $11.3\mathrm{\%}$ and $56.8\mathrm{\%}$ respectively compared with the variation of RRT in Wu and Hu (2020).