Trajectory planning of unmanned aerial vehicle based on the improved biogeography-based optimization algorithm

Abstract

This paper is concerned with trajectory planning for unmanned aerial vehicle in a three-dimensional complex workspace. Biogeography-based optimization algorithm is widely used in solving practical problems because of its fewer parameters, fast convergence rate, and good global optimization ability. In this paper, some improvements, modifying migration and mutation operations are made on the biogeography-based optimization algorithm to make it suitable for solving the trajectory planning problem. The optimal trajectory obtained by the improved algorithm can be used to generate the reference trajectory. Then, a control scheme of unmanned aerial vehicle based on the Lyapunov theory and radial basis function neural network is formed to track the reference trajectory. The improved trajectory algorithm generates the shortest trajectory and the time consumption is the lowest. Finally, the designed control scheme makes the unmanned aerial vehicle track the different trajectories quite well, the effectiveness of it can be illustrated by algorithm accuracy and electricity consumption.

Introduction

Unmanned aerial vehicle (UAV) is often seen for reconnaissance, search, and rescue missions in dangerous areas because of its high mobility, low risk, and low cost.^1,2 In order to ensure the efficiency and safety of UAV in the task completion, it is usually necessary to estimate an optimal flight trajectory before a mission execution. The trajectory planning for UAV in complex environments to avoid a series of obstacles is really a challenging problem. For these reasons, UAV trajectory planning has become a research hotspot.^1–6

Samaniego et al. ⁷ regarded the space of UAV flight as a discrete adaptive grid with minimum step size, and used cost response algorithm and recursive calculation to conduct trajectory planning, which reduced the response time, Lupascu et al. ⁸ extended the square in two dimensions to the cube in three dimensions for trajectory planning, and demonstrated the feasibility of the algorithm. Contreras-Cruz et al. ⁹ trajectory planning in a static environment using artificial bee colony algorithm was proposed. Firstly, artificial bee colony algorithm was used to explore the local environment, and then evolutionary algorithm was used to obtain the optimal trajectory. Li et al. ¹⁰ used the balance strategy to make an improvement on the artificial bee colony algorithm, and collected the convergence information in the iterative process to control the exploration accuracy, so as to achieve a balance between local and global information. Particle swarm optimization algorithm with dynamic adjustment of inertia weight (WLIPSO) used in Lu et al. ¹¹ initialized a swarm of particles with a random solution. The cooperation and competition ideas were introduced to update its speed and it found the optimal solution by searching for the personal best and the global best. However, in the late stage of the convergence, the convergence speed becomes slow, all particles tend to be the same and lose diversity. When the algorithm converges to a certain degree, it cannot be further optimized and the accuracy it can achieve is not high. Biogeography-based optimization with covariance matrix based migration (CMM-BBO) ¹² is a novel evolutionary algorithm. The covariance matrix based migration operator combined with the original migration randomly, thus the population diversity is increased. CMM-BBO can get the optimal solution of the shortest trajectory, but it is inefficient because it traverses many nodes and performs many repeated operations. The above algorithms have some disadvantages, with the expansion of the planning area, their computational complexities will increase exponentially, resulting in a “combinatorial explosion,” and they are prone to fall into local optimal solutions, which needs to be avoided. Compared with these algorithms, the improved BBO algorithm can deal with the uncertainty in trajectory planning more effectively and improve the reliability of trajectory planning and the stability of the system.

Biogeography-based optimization (BBO) algorithm, proposed in Simon, ¹³ is a global optimization algorithm that studies the distribution of organisms in geographical areas. It uses migration operation to realize information sharing among solutions (habitats) and uses mutation operation to enrich the population diversity in habitats. ¹⁴ It describes the migration, formation, and extinction of species by the mathematical model based on biogeography. ¹⁵ The improved BBO algorithm is similar to particle swarm optimization (PSO) algorithm, ant colony optimization algorithm, and genetic algorithm, but it has fewer designed parameters, low computational complexity, and strong optimization performance.^3,13,16,17

According to the grid-based workspace modeling algorithm in Zhang et al., ¹⁸ Islam et al., ¹⁹ and He et al., ²⁰ the UAV’s workspace is divided into a series of small blocks by using grids through the three-dimensional space, and the vertices of each block are selected as trajectory nodes. Mark these nodes as collision and collision-free nodes by determining whether they are in obstacles or not.^21,22 Simply select a series of available trajectory nodes (collision-free nodes) in a certain order and put them into a trajectory node list, that list consists of some available trajectories. Each available trajectory can be regarded as a habitat in the biogeography algorithm, and its nodes can be regarded as suitability index variables (SIVs) ²³ of the habitat. The habitat suitability index (HSI) ²⁴ of each habitat can be obtained by taking its nodes as parameters into the cost function calculations. According to the obtained HSIs, the immigration rate, emigration rate, and mutation probability corresponding to each habitat can be calculated, which can be used for the migration operation and mutation operation.^25,26 Finally, the optimal trajectory is generated. In particular, an improved migration and mutation operation is proposed in section 2.D for UAV trajectory planning. In addition, WLIPSO and CMM-BBO algorithm are used to complete the task of the UAV trajectory planning in the same workspace, and then compared with the proposed algorithm.

The optimal trajectory generated by the proposed biogeography algorithm can be used to generate the flight reference trajectory for UAV. Then, a fixed-wing UAV mathematical model under external interference is constructed and selected as the research object. ²⁷ Finally, according to the generated reference trajectory, a robust controller designed based on Lyapunov stability theory and radial basis function neural network (RBFNN)^28,29 is used for the fixed-wing UAV tracking the reference trajectory.

Our main contributions are summarized as follows:

• In this paper, the algorithm of setting collision nodes is used in the 3D workspace to model the UAV operation environment, which effectively improves the operation efficiency of the trajectory planning algorithm in the 3D workspace.

• The novel trajectory planning algorithm based on biogeography algorithm improves the migration and mutation operations compared with the original biogeography algorithm, making it more suitable for UAV trajectory planning.

• According to the trajectory generated by the improved BBO algorithm proposed in this paper, an interpolation algorithm based on time step is designed, which can convert the discrete points in the workspace to a relatively continuous trajectory, and can effectively realize the UAV’s accurate tracking of the reference trajectory after combining with the designed UAV control scheme.

The arrangement of this paper is made as follows: In section 2, the improved trajectory planning algorithm based on BBO in the 3D workspace is presented. We give the mathematical model of the UAV in section 3. The control scheme of UAV with the planned trajectory is proposed in section 4. Section 5 presents the experimental details and shows the effectiveness of the improved BBO algorithm and the control scheme. We give a conclusion in section 6.

The improved trajectory planning algorithm based on BBO in three-dimensional environment

For ease of expression, the following notations are exploited in this paper:

1. $\mathit{R}$ denotes a real constant space, ${\mathit{R}}^n$ denotes an n-dimensional real vector space, ${\mathit{R}}^{m\times n}$ represents an $m\times n$ -dimensional real matrix pace.

2. ${•}^T$ represents the transpose operation of a matrix, ‖•‖ represents the Euclidean norm of a matrix and can be expressed as $\sqrt{\mathrm{tr}({•}^T•)}$ , and $\mathrm{tr}(•)$ represents the trace of a matrix.

3. There exists a pseudo inverse in a non-full rank matrix $\mathit{M}\in {\mathit{R}}^{m\times l}$ , $m\ne l$ , such that $\mathit{M}{\mathit{M}}^{+}={\mathit{I}}_m$ . ³⁰

BBO is a biology evolutionary algorithm using biogeographic dynamic models and ideas that describe the movement of species from habitat to habitat in an ecosystem and how species are created or lost. ³¹ In BBO, each habitat resembles an island in the sea, is regarded as a possible solution and has a HSI to test the suitability and the priority of the habitat. Other aspects, for example, plant and animal diversity, rainfall capacity, etc., which characterize HSI are called SIVs.^31,32

In the following proposed BBO based trajectory optimizing process, the habitat is considered as the trajectory in the trajectory planning, a SIV of each habitat as a route node, the HIS of each habitat as the cost function.

The construction of the UAV’s workspace

The first step of the vehicle’s trajectory planning in three-dimensional space is to construct the vehicle’s work environment in a meaningful way. According to He et al. ²⁰ and Garrido et al., ³³ the grid-based workspace modeling algorithm is introduced to describe the vehicle’s workspace. As shown in Figures 1 and 2, this algorithm cuts the UAV’s workspace into n_x × n_y × n_z blocks (n_x, n_y, and n_z are the number of cutting layers on the x, y, and z axes, respectively), the vertexes of each block are selected as the trajectory node, and the nodes in obstacles are defined as collision nodes, ³⁴ which means that UAV cannot pass through them. On the contrary, UAV can plan a trajectory from start to end by constructing a list only containing collision-free nodes.

Figure 1.

The UAV’s workspace in three-dimensional workspace.

Figure 2.

The UAV’s constructed workspace in three-dimensional space (red stars and blue dots represent collision nodes and collision-free nodes, respectively).

A trajectory generation algorithm based on the constructed workspace

Before operating BBO algorithm, some initial trajectories need to be generated randomly so that they can be used as initial habitats in BBO to facilitate the formation of optimal trajectories by immigration, mutation, and other operations.

As shown in Figure 3, according to the construction algorithm of UAV’s workspace described above, the node where the UAV is currently located can be marked as node(i), the i-th node in UAV’s constructed workspace, $i\in [1,2,...,n_x*n_y*n_z]$ . Where n_x, n_x n_x are the number of UAV’s workspace cutting layers on x-axis, y-axis, and z-axis in the inertial coordinate frame, respectively. The adjacent nodes on its x, y, and z axes can be marked as node(i − 1), node(i + 1), node(i − n_y), node(i + n_y), node(i − n_x × n_y), and node(i + n_x × n_y), respectively.

Figure 3.

The adjacent nodes of node(i) in x, y, and z axis.

The steps of the trajectory generation algorithm are as follows:

• Step 1. A trajectory node list PL is initialized as an empty list; the iteration times IR is initialized as zero, the maximum iteration times is set as IR_max; the current node and goal node is initialized as node(i) and node(j), $i,j\in [1,2,..,n_x*n_y*n_z]$ , respectively; an available node set VS is initialized as an empty set.

• Step 2. Firstly, the node(i) is selected into PL; Then, judging whether node(i) or node(x) is a collision node, if it is, return error information, else, obtain the three-dimensional position of node(i) and node(j), which is $(x_i,y_i,z_i)$ and $(x_j,y_j,z_j)$ , respectively. Here, $x$ , $y$ , and $z$ are x-axis position, y-axis position, and z-axis position in the inertial coordinate frame, respectively. If $x_j>x_i$ , node(i + 1) is selected into VS, else if $x_j<x_i$ , node(i − 1) is selected into VS, else, there is no node on the x-axis of node(i) can be selected. If $y_j>y_i$ , node(i + n_y) is selected into VS, else if $y_j<y_i$ , node(i − n_y) is selected into VS, else, there is no node on the y-axis of node(i) can be selected. If $z_j>z_i$ , node(i + n_x × n_y) is selected into VS, else if $z_j<z_i$ , node(i − n_x × n_y) is selected into VS, else, there is no node on the z-axis of node(i) can be selected.

• Step 3. Firstly, the collision nodes in VS is eliminated; then, one node in VS is randomly selected as the next position of UAV; particularly, an adjacent node of node(i) that is not collision node is randomly selected as the next position of UAV, when there is no node in VS. The randomly selected node is set as node(k), and IR = IR + 1.

• Step 4. To avoid duplicate nodes in PL, node(i) is set as a collision node, then the selected node in Step 3 is set as the current node by letting i = k. If three-dimensional position of the current node is equal to goal node or IR > IR_max, save the trajectory list PL and return success information, else, skip to Step 2.

To illustrate the trajectory generation algorithm visually, all available trajectories generated from the current node to the goal node are shown in Figure 4.

Figure 4.

The trajectories generated by the proposed trajectory generation algorithm from the current node to the goal node.

Cost function

In the process of UAV trajectory planning, choosing an optimal trajectory is the ultimate goal. Usually, there is a cost function to measure whether the current trajectory is optimal. If the algorithm has a small value of cost function, then the computational complexity is small. Compared with other trajectories, the generated trajectory with the smallest value of cost function will be considered as the optimal trajectory. ⁵ The cost function is defined as follows:

$$F_{\mathrm{cost}}=N_{\mathrm{number}}+P_{\mathrm{length}}+P_{\mathrm{smoothing}}$$

(1)

In the cost function, $N_{\mathrm{number}}$ represents the number of nodes in the trajectory’s node list.

The term related to the length of the planned trajectory is represented as follows:

$$P_{\mathrm{length}}=\sum _{i=1}^{N_{\mathrm{number}}-1}\Vert \mathrm{pos}(i)-\mathrm{pos}(i+1)\Vert$$

(2)

Where $\mathrm{pos}(i)=\left(x_i{y}_i{z}_i\right)$ , $i\in [1,2,...,N_{\mathrm{number}}]$ is three-dimensional position of the node(i) in the trajectory node list.

The term related to the smoothness of the planned trajectory is represented as follows:

$$P_{\mathrm{smoothing}}=\sum _{i=1}^{N_{\mathrm{number}}-1}\frac{1}{\arccos (-q_{i+1}·q_i)}$$

(3)

Where $q_i=(\mathrm{pos}(i)-\mathrm{pos}(i-1))/\Vert \mathrm{pos}(i)-\mathrm{pos}(i-1)\Vert$ , $i\in [1,2,..,N_{\mathrm{number}}]$ .

The migration and mutation operation

Migration

The i-th habitat PL _i (a trajectory node list of the i-th habitat) and j-th habitat PL _j (a trajectory node list of the j-th habitat) are selected for immigration and emigration operation respectively, $i,j\in [1,2,..,\mathrm{Np}]$ , where $\mathrm{Np}$ is the number of habitats or trajectories. Before the migration operation, in order to improve the diversity of species in each habitat, if there is a same node node(k), $k\in [1,2,...,n_x*n_y*n_z]$ , between PL _i and PL _j , all nodes in PL _i after node(k) will be cleared. Then, node(m) is randomly selected from PL _j and use it as the start node to generate a new trajectory to the goal node by using the trajectory generation algorithm proposed in section 2.B. Similarly, node(p) is randomly selected from PL _i and use it as the start node to generate a new trajectory to node(m). In the end, clear the node(p) and all subsequent nodes in PL _i and put the two new generated trajectories into PL _i .

Mutation

The i-th habitat PL _i is selected for mutation operation. Firstly, a node node(m) that is not a collision code is randomly selected from the UAV’s workspace in Figure 2 and use it as the start node to generate a new trajectory to the goal node by using the trajectory generation algorithm, then, similarly, node(p) is randomly selected from PL _i and use it as the start node to generate a new trajectory to the node(m). In the end, clear the node(p) and all subsequent nodes in PL _i and put the two new generated trajectories into PL _i .

An improved BBO algorithm

Each habitat has its specific emigration rate ${\mu }_k$ and immigration rate ${\lambda }_k$ in BBO,^13,32 which are the functions of the HSIs in habitats and can be represented as follows:

$${\mu }_k=E\frac{F_{\mathrm{cost}}^k}{F_{\mathrm{cost}}^{\mathrm{maximum}}},{\lambda }_k=I(1-\frac{F_{\mathrm{cost}}^k}{F_{\mathrm{cost}}^{\mathrm{maximum}}})$$

(4)

Where $E$ and $I$ denote the maximum value of the emigration and immigration rate, respectively. $F_{\mathrm{cost}}^k$ is the k-th habitat’s cost function, and $F_{\mathrm{cost}}^{maximum}=\max [F_{\mathrm{cost}}^1,F_{\mathrm{cost}}^2,\dots ,F_{\mathrm{cost}}^{Np}]$ .

Now, the variable $P_k$ denotes the probability of the habitat’s population migration, and can be updated as $P_k=P_k+{\stackrel{·}{P}}_k$ , and according to Mo and Li ³² and Huang et al., ³⁵ the updating laws of the probability ${\stackrel{·}{P}}_k$ are defined as:

$${\stackrel{·}{P}}_k=\{\begin{array}{c}-({\lambda }_k+{\mu }_k)P_k+{\mu }_{k+1}{P}_{k+1},\hfill \\ \multicolumn{1}{c}{-({\lambda }_k+{\mu }_k)P_k+{\lambda }_{k-1}{P}_{k-1}+{\mu }_{k+1}{P}_{k+1}}\\ \multicolumn{1}{c}{-({\lambda }_k+{\mu }_k)P_k+{\lambda }_{k-1}{P}_{k-1},}\end{array},\begin{array}{c}k=0\hfill \\ \multicolumn{1}{c}{1\le k<\mathrm{Np}}\\ \multicolumn{1}{c}{k=\mathrm{Np}}\end{array}$$

(5)

The mutation probability $m_k$ for each habitat’s population is defined as:

$$m_k=m_k(1-\frac{P_k}{P_{max}})$$

(6)

Where $k\in [1,2,...,\mathrm{Np}]$ , $P_{max}=max[P_1,P_2,..,P_{\mathrm{Np}}]$ .

The steps of the improved BBO algorithm are summarized as follows:

• Step 1. The trajectory generation algorithm proposed in section 2.2 is used to randomly generate Np trajectories, which are used as the initial habitats, between the start node and the goal node. These parameters are initialized: the max emigration rate E, the max immigration rate I, the probability $P_k$ , and the mutation probability $m_k$ of each habitat. In the end, the iteration times NR of the improved BBO algorithm is initialized as zero, the maximum iteration times is set as NR_max.

• Step 2. The cost function $F_{\mathrm{cost}}^k$ , the immigration ${\lambda }_k$ , and emigration rate ${\mu }_k$ , the population migration probability $P_k$ and the mutation probability $m_k$ of each habitat are calculated in this iteration process. The optimal trajectory, that is the trajectory with the lowest cost function, is saved.

• Step 3. If ${\lambda }_i>\mathrm{rand}(0,d)$ , ${\mu }_j>\mathrm{rand}(e,1)$ , $0<d,1>e$ , $i\ne j$ , $i,j\in [1,2,...,\mathrm{Np}]$ are satisfied, the i-th habitat and the j-th habitat are selected for migration operation in section 2.D; especially, the i-th habitat is selected for the immigrant operation and the j-th habitat is selected for the emigrant operation; besides, the upper limit $d$ and the lower limit $e$ are selected to ensure that elites in habitats do not lose. Then, if $m_i>\mathrm{rand}(0,1)$ , the i-th habitat can be selected for mutation operation, and NR = NR + 1.

• Step 4. If NR > NR_max, save the optimal trajectory and return the success information, else, skip to Step 2.

The mathematical model of the fixed-wing UAV

Usually, the UAV’s motion in three-dimensional space needs to be transformed among three coordinate frames: body coordinate frame, inertial coordinate frame, and airflow coordinate frame, expressed in letters n, b, w, respectively.^36–38 The relationships between three coordinate frames are shown in Figure 5.

Figure 5.

The fixed-UAV in three-dimensional space.

Firstly, some definitions are given as: $\mathit{P}={(x,y,z)}^T$ and $\stackrel{·}{\mathit{P}}={(\stackrel{·}{x},\stackrel{·}{y},\stackrel{·}{z})}^T$ represent the position vector and velocity vector in the inertial coordinate frame, respectively. $\mathit{E}\mathbf{=}{(\phi ,\theta ,\psi )}^T$ and $\stackrel{·}{\mathit{E}}={(\stackrel{·}{\phi },\stackrel{·}{\theta },\stackrel{·}{\psi })}^T$ are the attitude angle and attitude angular velocity vectors in the inertial coordinate frame, respectively. $\phi ,\theta ,\psi$ are the roll angle, pitch angle, and yaw angle in the inertial coordinate frame, respectively. ${\mathit{V}}^b={(u,v,w)}^T$ and ${\mathit{\Omega }}^b={(p,q,r)}^T$ are the linear velocity and attitude angle velocity vectors in the body coordinate frame, respectively. $u,v,w$ are linear velocities in the body coordinate frame of x-axis, y-axis, and z-axis, respectively. The relationships between them are given as:

$${\mathit{V}}^b={\mathbf{R}}_n^b\stackrel{·}{\mathit{P}},{\mathit{\Omega }}^b={\mathbf{C}}_n^b\stackrel{·}{\mathit{E}}$$

(7)

The coordinate transform matrices ${\mathit{R}}_n^b$ and ${\mathit{C}}_n^b$ are given in Appendix 1. According to Newtonian mechanics,^27,39 the UAV’s mathematical model can be defined as:

$${\stackrel{·}{\mathit{V}}}^b=\frac{({\mathit{G}}^b+{\mathit{F}}_{\mathrm{engine}}^b+{\mathbf{C}}_w^b{\mathit{F}}_{\mathrm{aerodynamic}}^w)}{d}-{\mathit{\Omega }}^b\times {\mathit{V}}^b$$

(8)

$${\stackrel{·}{\mathit{\Omega }}}^b=({\mathit{I}}^b{)}^{-1}({\mathit{M}}^b-{\mathit{\Omega }}^b\times ({\mathit{I}}^b·{\mathit{\Omega }}^b))$$

(9)

The coordinate transform matrix ${\mathbf{C}}_w^b$ is given in Appendix 1. The forces ${\mathit{G}}^b$ , ${\mathit{F}}_{\mathrm{engine}}^b$ , and ${\mathit{F}}_{\mathrm{aerodynamic}}^w$ acting on the UAV are given in Appendix 2. The moment ${\mathit{M}}^{\mathbf{b}}$ and the inertial matrix ${\mathit{I}}^{\mathbf{b}}$ are given in Appendix 3.

Let ${\mathit{u}}_v={(F_u,F_v,F_w)}^T$ , ${\mathit{u}}_{\mathrm{\Omega }}={({\delta }_{a1},{\delta }_{a2},{\delta }_{e1},{\delta }_{e2},{\delta }_r)}^T$ , $F_u,F_v,F_w$ are the UAV’s engine force on x-axis, y-axis, and z-axis in the body coordinate frame, respectively. ${\delta }_{a1}$ is the controller of the left aileron, ${\delta }_{a2}$ is the controller of the right aileron, ${\delta }_{e1}$ is the controller of the left elevator, ${\delta }_{e2}$ is the controller of the right elevator, and ${\delta }_r$ is the controller of the rudder. The UAV’s mathematical model is simplified as:

$${\stackrel{·}{\mathit{V}}}^b=\mathit{A}{\mathit{V}}^b+\mathit{B}{\mathit{u}}_{\mathrm{v}}+{\mathit{Q}}_v$$

(10)

$${\stackrel{·}{\mathit{\Omega }}}^b=\mathit{F}{\mathit{\Omega }}^b+\mathit{G}{\mathit{u}}_{\mathrm{\Omega }}+{\mathit{Q}}_{\mathit{\Omega }}$$

(11)

The complete form of $\mathit{A}$ , $\mathit{B}$ , ${\mathit{Q}}_v$ , $\mathit{F}$ , $\mathbf{G}$ , and ${\mathit{Q}}_{\mathit{\Omega }}$ are given in Appendix 4.

The external disturbances and model uncertainty issues are considered, equations (10) and (11) are rewritten as:

$${\stackrel{·}{\mathit{V}}}^b=\mathit{A}{\mathit{V}}^b+\mathit{B}{\mathit{u}}_{\mathrm{v}}+{\mathit{Q}}_v+{\mathit{f}}_v$$

(12)

$${\stackrel{·}{\mathit{\Omega }}}^b=\mathit{F}{\mathit{\Omega }}^b+\mathit{G}{\mathit{u}}_{\mathit{\Omega }}+{\mathit{Q}}_{\mathit{\Omega }}+{\mathit{f}}_{\mathrm{\Omega }}$$

(13)

Where ${\mathit{f}}_v={(f_{v1},f_{v2},f_{v3})}^T$ and ${\mathit{f}}_{\mathit{\Omega }}={(f_{\mathrm{\Omega }1},f_{\mathrm{\Omega }2},f_{\mathrm{\Omega }3})}^T$ are the unknown external disturbances or model uncertainties.

The control scheme of UAV with the planned trajectory

According to the trajectory optimization algorithm proposed in section 2.E, an optimal trajectory PL connecting the initial node and the goal node is obtained. In particular, this trajectory consists of a series of trajectory nodes.

The reference trajectory and its derivative required by UAV is as follows:

• Step 1. The start node in PL is initialized as the current node node(i), the air speed of UAV is set as $v_s$ , and the time step is set as h. The time t is initialized as zero. The time T is initialized as zero. According to equation (2), the total length of the trajectory PL is $P_{\mathrm{length}}$ . The time the UAV spends on this trajectory is $T_{\mathrm{all}}=P_{\mathrm{length}}/v_s$ .

• Step 2. The time it takes to fly from the current node(i) to the following node(i + 1) is designed as $T_s=\mathrm{norm}(\mathrm{pos}(i+1)-\mathrm{pos}(i))/v_s$ , $\mathrm{pos}(i)=\left(x_i{y}_i{z}_i\right)$ , and $\mathrm{pos}(i+1)=\left(x_{i+1}{y}_{i+1}{z}_{i+1}\right)$ are three-dimensional positions of the node(i) and node(i + 1) respectively. Let $x_s=0$ , $y_s=0$ , $z_s=0$ , while $T\le t<T+T_s$ , do the following:

$$\begin{array}{c}{\stackrel{·}{x}}_r(t)=(x_{i+1}-x_i)/T_s,{\stackrel{·}{y}}_r(t)=(y_{i+1}-y_i)/T_s,\hfill \end{array}$$

$$\begin{array}{c}{\stackrel{·}{z}}_r(t)=(z_{i+1}-z_i)/T_s\hfill \end{array}$$

$$x_s=x_s+h*{\stackrel{·}{x}}_r,y_s=y_s+h*{\stackrel{·}{y}}_r,z_s=z_s+h*{\stackrel{·}{z}}_r,$$

$$x_r(t)=x_i+x_s,y_r(t)=y_i+y_s,z_r(t)=z_i+z_s,t=t+h$$

• Step 3. The time T is updated by $T=T+T_s$ , the current node is updated by letting $i=i+1$ , that is transferred from the i-th node to the i + 1-th node in PL. If node(i) is equal to the goal node in PL, return the reference trajectory $x_r(t)$ , $y_r(t)$ , $z_r(t)$ and its derivatives ${\stackrel{·}{x}}_r(t)$ , ${\stackrel{·}{y}}_r(t)$ , ${\stackrel{·}{z}}_r(t)$ between $t=0$ and $t=T_{\mathrm{all}}$ , else, skip to Step 2.

Lemma 1^40,42: Let $R(t)\in \mathit{R}$ be a continuous positive function with a bounded initial $R(0)$ . If the inequality holds $\stackrel{·}{R}(t)\le -\beta R(t)+\eta$ , where $\beta$ and $\eta$ are positive constants, then the following inequality is obtained:

$$R(t)\le R(0)e^{-\beta t}+\frac{\eta }{\beta }(1-e^{-\beta t})$$

Lemma2^28,42 : $x(t),y(t)\in {\mathit{R}}^m$ are bounded functions with continuous values, the following formula can be obtained by Young’s inequality and Cauchy inequality:

$$\begin{array}{c}\multicolumn{1}{c}{x(t)y(t)\le \Vert x(t)\Vert \Vert y(t)\Vert \le \frac{{\Vert x(t)\Vert }^2}{\lambda }}\\ +\frac{{\Vert y(t)\Vert }^2}{\gamma },\lambda ,\gamma >0,\frac{1}{\lambda }+\frac{1}{\gamma }=1\hfill \end{array}$$

The control scheme for UAV’s attitude and position is designed as follows:

Let ${\mathit{P}}_r={(x_r,y_r,z_r)}^T$ , ${\stackrel{·}{\mathit{P}}}_r={({\stackrel{·}{x}}_r,{\stackrel{·}{y}}_r,{\stackrel{·}{z}}_r)}^T$ , ${\mathit{V}}_r^b={\mathit{R}}_n^b{\stackrel{·}{\mathit{P}}}_r$ , and $\stackrel{·}{\mathit{P}}={\left({\mathit{R}}_n^b\right)}^{-1}{\mathit{V}}^b$ . The error dynamic between the agent and the reference is represented as:

$$\mathit{e}=\mathit{P}-{\mathit{P}}_r+{\mathit{V}}^b-{\mathit{V}}_r^b$$

(14)

A Lyapunov function is given as follows:

$${\mathit{V}}_1=\frac{1}{2}{\mathit{e}}^T\mathit{e}$$

(15)

According to equation (12), the derivative of ${\mathit{V}}_1$ is obtained:

$${\stackrel{·}{\mathit{V}}}_1={\mathit{e}}^T\stackrel{·}{\mathit{e}}={\mathit{e}}^T(\mathbf{A}{\mathit{V}}^b+\mathbf{B}{\mathit{u}}_{\mathrm{v}}+{\mathbf{Q}}_v+{\mathit{f}}_v+{\mathit{f}}_l)$$

(16)

Where ${\mathit{f}}_l=\mathit{P}-{\mathit{P}}_r-{\mathit{V}}_r^b$ is a smooth function with bounded area, and satisfies $\Vert {\mathit{f}}_l\Vert \le a$ . Here, $a$ is a positive value. According to Lemma 2, the following one is obtained:

$${\mathit{e}}^T{\mathit{f}}_l\le a{\mathit{e}}^T\le \frac{a^2}{2}+\frac{{\Vert \mathit{e}\Vert }^2}{2}$$

(17)

Use equation (17), equation (16) is rewritten as:

$${\stackrel{·}{\mathit{V}}}_1\le {\mathit{e}}^T\mathbf{A}{\mathit{V}}^b+{\mathit{e}}^T\mathbf{B}{\mathit{u}}_{\mathrm{v}}+{\mathit{e}}^T{\mathbf{Q}}_v+{\mathit{e}}^T{\mathit{f}}_v+\frac{{\mathit{e}}^T\mathit{e}}{2}+\frac{a^2}{2}$$

(18)

To complete the UAV’s task, the desired controller ${\mathit{u}}_{\mathrm{v}}$ is designed as:

$${\mathit{u}}_{\mathrm{v}}=-\mathit{ke}-{\mathbf{B}}^{+}(\mathbf{A}{\mathit{V}}^b+{\mathbf{Q}}_v+{\mathit{f}}_v)$$

(19)

Let $k=\gamma (1/2+\raisebox{1ex}{$1$}\!\left/ \!\raisebox{-1ex}{$2\gamma $}\right.)B^{+}$ , $\gamma >0$ , and substitute it into equation (16), yields:

$${\stackrel{·}{\mathit{V}}}_1\le -\gamma {\mathit{V}}_1+\frac{a^2}{2}$$

(20)

Use Lemma 1, one can obtain that:

$${\mathit{V}}_1(t)\le \frac{a^2}{2\gamma }+({\mathit{V}}_1(0)-\frac{a^2}{2\gamma })e^{-\gamma t}$$

(21)

${\mathit{f}}_v$ is unknown, thus the controller in equation (19) cannot be introduced to equation (12) immediately. According to Luan et al. ⁴¹ and Tang et al., ⁴² the RBFNN is often used to estimate the unknown function.

$${\mathit{f}}_v={\mathit{w}}^{*T}\mathit{\eta }(\mathit{z})+\mathit{\epsilon }(\mathit{z}),{\hat{\mathit{f}}}_v={\hat{\mathit{w}}}^T\mathit{\eta }(\mathit{z})$$

(22)

where ${\mathit{w}}^{*}\in {\mathit{R}}^{N\times p}$ is the optimal weight, $\mathit{\epsilon }(\mathit{z})$ is an estimated error, and $\Vert \mathit{\epsilon }(\mathit{z})\Vert \le \delta$ , $\hat{\mathit{w}}\in {\mathit{R}}^{N\times p}$ is used to estimate the optimal weight, and the error between them can be defined as $\tilde{\mathit{w}}=\hat{\mathit{w}}-{\mathit{w}}^{*}$ . ${\mathit{\eta }}_i(\mathit{z})=\exp (-\Vert \mathit{z}-{\mathbf{\xi }}_i{\Vert }^2/{\zeta }_i^2)$ and $\mathit{\eta }(\mathit{z})={[{\mathit{\eta }}_1(\mathit{z}),...,{\mathit{\eta }}_N(\mathit{z})]}^T\in {\mathit{R}}^N$ . ${\mathbf{\xi }}_i$ is the center in the receptive field and ${\zeta }_i$ is the width of the receptive field. $\mathit{z}=\mathit{e}\in {\mathit{R}}^m$ is the input unknown disturbance.

Then, a redesigned Lyapunov function is generated as:

$$\mathit{V}={\mathit{V}}_1+\frac{1}{2\eta }\mathrm{tr}({\tilde{\mathit{w}}}^T\tilde{\mathit{w}})$$

(23)

By combining equations (12), (17), and (22), the derivative of equation (23) is obtained as:

$$\begin{array}{c}\multicolumn{1}{c}{\stackrel{·}{\mathit{V}}\le {\mathit{e}}^T(\mathbf{A}{\mathit{V}}^b+\mathbf{B}{\mathit{u}}_{\mathrm{v}}+{\mathbf{Q}}_v+{\mathit{w}}^{*T}\mathit{\eta }(\mathit{z})+\mathit{\epsilon }(\mathit{z}))}\\ +\frac{{\mathit{e}}^T\mathit{e}}{2}+\frac{a^2}{2}+\frac{1}{\eta }\mathrm{tr}({\tilde{\mathit{w}}}^T\stackrel{·}{\stackrel{~}{\mathit{w}}})\hfill \end{array}$$

(24)

According to Lemma 2, the following one is obtained:

$${\mathit{e}}^T\mathit{\epsilon }(\mathit{z})\le \delta {\mathit{e}}^T\le \frac{{\delta }^2}{2}+\frac{{\Vert \mathit{e}\Vert }^2}{2}$$

(25)

Use equation (25), equation (24) is rewritten as:

$$\begin{array}{c}\multicolumn{1}{c}{\stackrel{·}{\mathit{V}}\le {\mathit{e}}^T(\mathbf{A}{\mathit{V}}^b+\mathbf{B}{\mathit{u}}_{\mathrm{v}}+{\mathbf{Q}}_v+{\mathit{w}}^{*T}\mathit{\eta }(\mathit{z}))}\\ +{\mathit{e}}^T\mathit{e}+\frac{a^2+{\delta }^2}{2}+\frac{1}{\eta }\mathrm{tr}({\tilde{\mathit{w}}}^T\stackrel{·}{\stackrel{~}{\mathit{w}}})\hfill \end{array}$$

(26)

To finish the control task, the controller is redesigned as:

$${\mathit{u}}_{\mathrm{v}}=-\mathit{ke}-{\mathbf{B}}^{+}(\mathbf{A}{\mathit{V}}^b+{\mathbf{Q}}_v+{\hat{\mathit{w}}}^T\mathit{\eta }(\mathit{z}))$$

(27)

Where $\mathit{k}=\gamma (1/2+1/\gamma ){\mathbf{B}}^{+}$ , $\gamma >0$ , and the weight adaptive law of the neural network is designed as:

$$\stackrel{·}{\hat{\mathit{w}}}=\stackrel{·}{\stackrel{~}{\mathit{w}}}=\eta (\mathit{\eta }(\mathit{z})e^T-\sigma \hat{\mathit{w}})$$

(28)

Substitute equations (27) and (28) into equation (26), yields:

$$\begin{array}{c}\multicolumn{1}{c}{\stackrel{·}{\mathit{V}}\le -\frac{\gamma }{2}{\mathit{e}}^T\mathit{e}-{\mathit{e}}^T{\tilde{\mathit{w}}}^T\mathit{\eta }(\mathit{z})+\mathrm{tr}({\tilde{\mathit{w}}}^T\mathit{\eta }(\mathit{z})e^T-\sigma {\tilde{\mathit{w}}}^T\hat{\mathit{w}})}\\ +\frac{a^2+{\delta }^2}{2}\hfill \\ =-\frac{\gamma }{2}{\mathit{e}}^T\mathit{e}-\sigma \mathrm{tr}\left({\tilde{\mathit{w}}}^T\hat{\mathit{w}}\right)+\frac{a^2+{\delta }^2}{2}\hfill \end{array}$$

(29)

Based on the fact that ²² :

$${\tilde{\mathit{w}}}^{\mathrm{T}}\widehat{\mathit{w}}={\tilde{\mathit{w}}}^T\tilde{\mathit{w}}/2+{\widehat{\mathit{w}}}^T\widehat{\mathit{w}}/2-{\mathit{w}}^{\ast T}{\mathit{w}}^{\ast }/2$$

(30)

The following result is obtained:

$$-\sigma \mathrm{tr}({\tilde{\mathit{w}}}^T\hat{\mathit{w}})\le -\frac{1}{2}\sigma \mathrm{tr}({\tilde{\mathit{w}}}^T\tilde{\mathit{w}})+\frac{1}{2}\sigma \mathrm{tr}({\mathit{w}}^{*T}{\mathit{w}}^{*})$$

(31)

So, equation (29) is rewritten as:

$$\begin{array}{c}\multicolumn{1}{c}{\stackrel{·}{\mathit{V}}\le -\frac{\gamma }{2}{\mathit{e}}^T\mathit{e}-\sigma \eta \frac{1}{2\eta }\mathrm{tr}({\tilde{\mathit{w}}}^T\tilde{\mathit{w}})+\frac{1}{2}\sigma \mathrm{tr}({\mathit{w}}^{*T}{\mathit{w}}^{*})}\\ +\frac{a^2+{\delta }^2}{2}\hfill \\ \le -\overline{\beta }\mathit{V}+\overline{\eta }\hfill \end{array}$$

(32)

Where $\overline{\beta }=min[\gamma ,\sigma \eta ]$ , $\overline{\eta }=\frac{1}{2}\sigma \mathrm{tr}({\mathit{w}}^{*T}{\mathit{w}}^{*})+\frac{a^2+{\delta }^2}{2}$ .

Use Lemma 1, the following one can be obtained:

$$\mathit{V}(t)\le \frac{\overline{\eta }}{\overline{\beta }}+(\mathit{V}(0)-\frac{\overline{\eta }}{\overline{\beta }})e^{-\overline{\beta }t}$$

(33)

From it one can conclude that the trajectory tracking error can become very small in the situation of selecting the suitable parameters and choosing RBFNN to estimate the unknown disturbance.

Similarly, the controller ${\mathit{u}}_{\mathrm{\Omega }}$ in equation (13) is designed as:

$${\mathit{u}}_{\mathrm{\Omega }}=-\mathit{k}{\mathit{e}}_{\mathrm{\Omega }}-{\mathbf{G}}^{+}(\mathbf{F}{\mathit{\Omega }}^b+{\mathbf{Q}}_{\mathrm{\Omega }}+{\hat{\mathit{w}}}^T\mathit{\eta }(\mathit{z}))$$

(34)

Where $\mathit{k}=\gamma (1/2+1/\gamma ){\mathbf{G}}^{+}$ , $\stackrel{·}{\hat{\mathit{w}}}=\eta (\mathit{\eta }(\mathit{z}){\mathit{e}}_{\mathrm{\Omega }}^T-\sigma \hat{\mathit{w}})$ , $\mathit{z}={\mathit{e}}_{\mathrm{\Omega }}$ , ${\mathit{e}}_{\mathrm{\Omega }}=\mathit{E}-{\mathit{E}}_r+{\mathit{\Omega }}^b-{\mathit{\Omega }}_r^b$ , ${\mathit{\Omega }}_r^b={\mathbf{C}}_n^b{\stackrel{·}{\mathit{E}}}_r$ . The elements of the reference attitude vector ${\mathit{E}}_r={({\phi }_r,{\theta }_r,{\psi }_r)}^T$ related to the inertial coordinate frame are follows: ${\phi }_r=\arcsin (u/V_T)$ , ${\theta }_r=\arcsin (w/V_T)$ , and ${\psi }_r=\arcsin (v/V_T)$ .

To verify the RBFNN’s performance of estimating the external disturbance, the contrasted controllers without RBFNN of equations (27) and (34) are set as:

$${\mathit{u}}_{\mathrm{v}}=-\mathit{ke}-{\mathit{B}}^{+}(\mathit{A}{\mathit{V}}^b+{\mathit{Q}}_v)$$

(35)

$${\mathit{u}}_{\mathrm{\Omega }}=-\mathit{k}{\mathit{e}}_{\mathrm{\Omega }}-{\mathit{G}}^{+}(\mathit{F}{\mathit{\Omega }}^b+{\mathit{Q}}_{\mathrm{\Omega }})$$

(36)

Simulation results

The experiment is conducted on a personal computer of Inter Core I7 eight-core processor, 16 GB memory, 1T hard disk, and MATLAB R2018a is used for simulation analysis. The effectiveness of the proposed trajectory planning algorithm and UAV control scheme is verified in this section. The specific simulation process is represented as follows: Firstly, the map in Figure 1 is selected as the UAV’s workspace, the coordinates $\mathrm{pos}(\mathrm{start})=(5,47,0)$ and $\mathrm{pos}(\mathrm{goal})=(1,1,0)$ are selected as the of the start node and the goal node, respectively. Then the proposed trajectory planning algorithm is compared with the state-of-the-art competitive algorithms, such as WLIPSO and CMM-BBO. Secondly, a reference trajectory is generated by the proposed BBO based trajectory planning algorithm. Finally, the proposed control scheme is used to track the reference trajectory to complete the task of transportation. The detailed default parameters can be found in Appendix 5, Appendix 6, and Appendix 7. The specific simulation results are as follows:

The simulation results of trajectory planning using WLIPSO, CMM-BBO, and the improved BBO algorithm can be obtained as follows. The optimal trajectory length obtained by WLIPSO is 108 m, and it takes 7.1 s to solve the trajectory planning problem. The CMM-BBO algorithm is 99 m and it takes 6.7 s. The improved BBO is 94 m and it takes 4.8 s. After 30 iterations, the values of the cost function of WLIPSO, CMM-BBO, and improved BBO are 93.5, 90.3, and 88.2, respectively. The improved BBO achieves the lowest value of cost function, the computational complexity is the lowest compared with other two algorithms. After obtaining the optimal path quickly, the UAV follow the optimal path to complete the specific task. It can be concluded that the improved BBO based trajectory planning algorithm has good performance, and the trajectories planned by these algorithms in three-dimensional workspace are shown in Figure 6.

Figure 6.

Trajectory planning using different algorithms.

As shown in Figures 7 and 8, the improved BBO based trajectory planning algorithm is used to find an optimal trajectory connecting the start node and the goal node in the complex working environment. The trajectory generated by the improved BBO is selected as the reference trajectory. From the two-dimensional plane and the three-dimensional plane, it can be seen that the actual flight trajectory of the UAV almost coincides with the reference trajectory. From Figure 9, it can be seen that the control scheme without RFBNN is bad in following the reference trajectory. The reason why the proposed control scheme has good tracking effect is that it can better compensate for the compound interference terms and eliminate the influence of various interferences in the flight process.

Figure 7.

Trajectory planning and following in three-dimensional workspace.

Figure 8.

The vertical view of UAV trajectory tracking.

Figure 9.

The vertical view of UAV trajectory tracking without RBFNN.

The trajectory tracking effects of different methods can be seen from Figures 10 to 15. By using the proposed control scheme, the UAV can track the trajectory quite well in three different situations. The position deviation and the attitude deviation are small in all tracking process. The algorithm accuracies of positions generated by WLIPSO, CMM-BBO, and improved BBO are 96.30%, 96.86%, and 96.48%, respectively. The average algorithm accuracy is 96.55%. The electricity consumption is a crucial problem for UAV tracking. The time consumption of UAV in tracking the trajectories generated by WLIPSO, CMM-BBO, and improved BBO are 20.40, 18.80, and 18.20 s, respectively. Among the three algorithms, the trajectory length and time consumption of the proposed algorithm is the lowest, and the electricity consumption of it is the lowest.

Figure 10.

UAV’s positions tracking of improved BBO.

Figure 11.

UAV’s attitudes tracking of improved BBO.

Figure 12.

UAV’s positions tracking of WLIPSO.

Figure 13.

UAV’s attitudes tracking of WLIPSO.

Figure 14.

UAV’s positions tracking of CMM-BBO.

Figure 15.

UAV’s attitudes tracking of CMM-BBO.

To evaluate the good tracking performance of the designed controller, controllers with and without RFBNN are compared. The tracking effects of positions and attitudes in three directions are illustrated from Figures 10 to 15. The controllers with neural network proposed in equations (27) and (34) are used to track the different trajectories, the control effects are better than controllers without neural network proposed in equations (35) and (36). That is because the neural network used in this paper can accurately estimate the unknown external disturbance. From Figures 16 to 17, it can be concluded that the estimated values are close to the actual values, the neural network has good estimation ability.

Figure 16.

The disturbances in equation (12) and its estimates.

Figure 17.

The disturbances in equation (13) and its estimates.

Conclusions

This paper proposes an off-line trajectory planning algorithm and a robust trajectory tracking control scheme for fixed-wing UAV. In this paper, an improved trajectory planning algorithm based on BBO is proposed to generate an optimal trajectory connecting the start node and the goal node for obstacle collision avoidance. By comparing with WLIPSO and CMM-BBO algorithm, it is concluded that this algorithm has certain advantages in UAV trajectory planning. A robust controller with RBFNN to estimate unknown disturbances is proposed to track the trajectory generated by three algorithms. Compared with other two algorithms, the improved BBO algorithm and the designed control scheme generate the optimal trajectory in the subjective and objective aspects. In the following work, the reduction of the algorithm complexity will be considered.