Three-Dimensional Path Planning of Constant Thrust Unmanned Aerial Vehicle Based on Artificial Fluid Method

In this paper, a three-dimensional path planning problem of an unmanned aerial vehicle under constant thrust is studied based on the artificial fluid method.'e effect of obstacles on the original fluid field is quantified by the perturbationmatrix, the streamlines can be regarded as the planned path for the unmanned aerial vehicle, and the tangential vector and the disturbance matrix of the artificial fluid method are improved. In particular, this paper addresses a novel algorithm of constant thrust fitting which is proposed through the impulse compensation, and then the constant thrust switching control scheme based on the isochronous interpolation method is given. It is proved that the planned path can avoid all obstacles smoothly and swiftly and reach the destination eventually. Simulation results demonstrate the effectiveness of this method.


Introduction
Unmanned aerial vehicles have the characteristics of high flexibility, low cost, high safety, and strong concealment. ese unique superior performances have enabled the rapid development of unmanned aerial vehicle technology and gradually become the representative technology of the world's cutting-edge technology. ey are widely used in various fields such as military reconnaissance, urban express delivery, terrain exploration, and environmental monitoring. e purpose of the unmanned aerial vehicle three-dimensional path planning is to find the optimal path between the initial position and the target position under the constraints of the unmanned aerial vehicle and the environment.
However, unmanned aerial vehicles still face a challenge; that is, it is difficult to ensure the safety and reliability of flight paths because unmanned aerial vehicles encounter obstacles and threats, such as high-rise buildings and enemy air-defense systems. To address the challenge, a path planning problem has been actively studied, which derives a flight path for an unmanned aerial vehicle from a point to another with respect to navigating a region of interest safely. Various operational constraints such as the maximum energy level of an unmanned aerial vehicle and safety distance from an object (e.g., a building or another unmanned aerial vehicle) are often imposed into the problem.
Following the importance of the path planning in unmanned aerial vehicle deployment, various approaches based on exact and heuristic approaches have been proposed. Model predicted control (MPC) method [1], Voronoi method [2], intelligent algorithms (e.g., genetic algorithm and particle swarm optimization) [3][4][5], rapidly exploring random tree (RRT) method [6], and artificial potential field (APF) [7] are some of the typical algorithms. e improved rapidly exploring random tree (RRT) method produces a time parameterized set of control inputs to make the robot move from the initial point to the destination, which proves to be efficient for 3D path planning. Artificial potential field (APF) method has the advantages of simple principle and real-time computation [8,9], but there exists local minimum when the robot enters into a concave area. Besides, it is hard to obtain a feasible path sometimes even if the magnitude of the attractive or repulsive force is regulated. e intelligent algorithms, such as particle swarm optimization [10], evolutionary algorithms [5], and ant colony algorithm [11], are also widely used in 3D path planning.
ese methods can be easily employed in different environments, but it is possible to trap in a local optimum. However, the abovementioned drawback can be relieved when the intelligent methods are improved or combined with other methods [12][13][14][15][16]. ese traditional approaches are improved to solve the three-dimensional path planning problem. However, the calculation of these algorithms tends to increase exponentially if the planning space enlarges. Besides, the planned path may be not smooth enough for the robot to track. As a result, extra strategy of the smoothing path is usually needed. Several research studies targeted 3D path planning in order to plan a feasible and smooth path. A novel algorithm based on the disturbed fluid and trajectory propagation is developed to solve the 3D path planning problem of an unmanned aerial vehicle in static environment [17]. e core path graph algorithm [18] calculates the core path graph where arcs are minimum-length trajectories satisfying geometrical constraints and searches the optimal trajectory between two arbitrary nodes of the graph. However, multiple quadratics should be resolved, resulting in low computational efficiency. In addition, constant thrust collision avoidance maneuver in path planning is studied in our previous studies [19,20].
In this paper, a novel algorithm of constant thrust fitting is proposed through the impulse compensation for constant thrust maneuver of an unmanned aerial vehicle, and the tangential vector and the disturbance matrix of the artificial fluid method is improved by combing the interfered fluid dynamical system. Although the physical characteristics of the modified streamlines are broadened, they still conform to the basic properties of fluid flow, i.e., smoothness, impenetrability, and accessibility. e rest of the paper is organized as follows. Section 2 focuses on the calculation of the shortest distance between the unmanned aerial vehicle and the obstacle. Section 3 explains 3D path planning based on the improved artificial fluid method. Section 4 describes constant thrust collision avoidance maneuver. e simulation results are given in Section 5. Section 6 concludes the paper.

Calculation of the Shortest Distance between
Unmanned Aerial Vehicle and Obstacle e purpose of an unmanned aerial vehicle is to avoid obstacles and reach the destination. In this paper, the obstacle is described approximately as a 3D space surface. We define the relative motion coordinate system o − xyz as the path planning space, where the origin is the center of the obstacle, and R � (x, y, z) is taken as the position of the unmanned aerial vehicle relative to the obstacle. Suppose the parametric equations of the obstacle's space curved surface as follows: where f u and f v are the partial derivatives of f(u, v) on u and v. Obviously, in order to get the minimum distance between R � R(u, v) and Q(u, v), RQ ��→ should be parallel to the normal vector n q , that is to say, Equation (3) can be written as follows: where f 1 , f 2 , and f 3 are the nonlinear equations. en, equation (4) can be written as follows: e derivative matrix of equation (5) can be written as the following form: Equation (6) is a first-order partial differential equation, and we can directly find its analytical solution or use the toolbox in MATLAB to get its analytical solution. We assume that X * (u * , v * ) is the analytical solution of equation (6), then the following results can be obtained by using the Taylor formula of multivariate functions: us, the least square solution of X * can be obtained as follows: erefore, the point Q * (u, v) on the obstacle's space curved surface minimizes the distance between the unmanned aerial vehicle and the obstacle. e shortest distance between the unmanned aerial vehicle to the obstacle is where |·| represents the modular of the vector from R(u, v) to Q * (u, v).

3D Path Planning Based on the Artificial Fluid Method
Based on the description in Section 2, the procedure for 3D path planning is as follows. First, the perturbation matrix P(u, v) is calculated [17]. Next, we calculate the disturbed fluid velocity v d by modifying the target velocity v T . en, the planned path is obtained by the recursive integration of v d . Finally, constant thrust collision avoidance maneuver is studied and the switching control scheme based on the isochronous interpolation method is given. To describe the influence of the obstacle on the original flow, the perturbation matrix P(u, v) is defined as follows: where I is a 2 × 2 identity matrix, n q is a column vector given by equation (2), f u is a tangential vector (the partial derivative of f(u, v) on u) at the point Q(u, v), λ(u, v) is a saturation function defining the orientation of tangential velocity, ‖·‖ is the 2-norm of a vector or a matrix, and ρ(u, v) and σ(u, v) are defined as the weight of n q and t q , respectively: where ρ 0 is the repulsive parameter, σ 0 is the tangential parameter, D 0 is the maximum radius of the unmanned aerial vehicle, and λ 0 is a small positive threshold of the saturation function λ(u, v). en, the disturbed fluid velocity v d can be calculated by

e Planned Path Can Avoid Obstacles Safely.
To avoid possible collisions, an unmanned aerial vehicle cannot approach obstacles indefinitely, so we introduce the maximum radius of the unmanned aerial vehicle D 0 , i.e., the distance between the boundary of the unmanned aerial vehicle and the obstacle should be greater than D 0 . Suppose that D * (u, v) � D 0 + δ(u, v) and δ(u, v) is a monotonically decreasing function: It can be inferred that Because vectors n q and f u are perpendicular exactly, i.e., n T q f u � 0, and the equation n T q f u � 0 means that n T q v d � 0, the path is outside of the minimum permitted distance and there is no collision.

e Planned Path Can Reach the Destination Eventually.
Because the goal of the path planning is to make the unmanned aerial vehicle reach the destination safely, the velocity of the unmanned aerial vehicle should have a component in the direction of the target velocity, i.e., velocity v T and v d should satisfy v T T v d ≥ 0, and the planned path will converge to the target point. Besides, v d ≈ v T should be satisfied if the unmanned aerial vehicle is near to the destination ξ T � (x T , y T , z T ):

Discrete Dynamics in Nature and Society
where (·) is the inner product of vectors, where 〈v T , n q 〉 denotes the angle between v T and n q . It is obvious that when the unmanned aerial vehicle approaches the destination, As D * (u, v) ρ(u,v) ⟶ + ∞ and cos 2 〈v T , n q 〉 ≤ 1 hold, we

Analysis of the Disturbed Fluid Velocity v d .
e modified velocity v d defined by equation (15) can be expressed as It can be seen from equation (20) that v d consists three parts: v T can be called the target velocity; ‖f u ‖‖n q ‖)f u can be called the tangential velocity. Similarly, the perturbation matrix P(u, v) can be divided into three parts: attractive matrix I, repulsive matrix n q n T q /D * (u, v) ρ(u,v) ‖n q ‖ 2 , and tangential matrix . It can be analyzed that the magnitudes of repulsive and tangential velocities increase with ρ(u, v) and σ(u, v), respectively. erefore, we can readjust the shape of the path by changing parameters ρ(u, v) and σ (u, v). is method is similar to the virtual force method or then the artificial potential field method to some degree. However, the perturbation matrix by this method can describe the effect of obstacles on path more objectively, considering the shape of obstacles and the position of the unmanned aerial vehicle.

Constant Thrust Collision Avoidance Maneuver
Suppose that the time of the unmanned aerial vehicle's collision avoidance maneuver is T and the shortest switching time interval is ΔT. ere are M shortest switching time intervals and N target maneuver positions, and T i represents the time of the i-th thrust arc. e process of collision avoidance maneuver can be considered as the system state variables change from a nonzero initial state x(0) to a desired state x(T) � 0: Suppose that the actual constant thrusts of the unmanned aerial vehicle is F, the maximum thrusts is F, and the theoretical continuous thrusts is F * . e thrusters can provide different sizes of constant thrust to meet different thrust requirements. ere are N different sizes of constant thrust which can be denoted as follows: e size of the constant thrust is calculated as follows: there are N + 1 thrust levels which can be selected; N), and the level of the constant thrust can be calculated as follows, taking the i-th thrust arc as example: where [] means bracket function and |F * (t)| means the absolute value of F * (t).

Constant rust Fitting through the Impulse Compensation.
e constant thrust fitting should be discussed in several categories; for convenience, let us take the i-th thrust arc as example.   1, 2, . . . , N), and the level of the constant thrust can be calculated as follows: where [] means bracket function and |F * (t)| means the absolute value of F * (t). en, the impulse error can be calculated as follows: Suppose that the value of the impulse compensation threshold is a positive constant c > 0: (1) If the impulse error ΔI i satisfies the following condition: and the actual constant thrust of the unmanned aerial vehicle can be calculated as follows: then the unmanned aerial vehicle will not carry out impulse compensation. (2) If the impulse error ΔI i satisfies the following condition: then the unmanned aerial vehicle should carry out impulse compensation, and the size of the constant thrust impulse compensation can be calculated as follows: Case 4. If the theoretical working time in the i-th thrust arc t * � M i ΔT, then the impulse error in the i-th thrust arc ΔI i can be calculated as follows: Discrete Dynamics in Nature and Society Furthermore, if there exist n 1 shortest switching time intervals satisfying the following conditions, without loss of generality, we suppose that these time intervals are the first n 1 shortest switching time intervals, taking the j-th shortest switching time interval as an example: then the size of the impulse compensation can be calculated as follows: (1) If the impulse error ΔI i satisfies the following condition: and the actual constant thrust of the unmanned aerial vehicle can be calculated as follows, taking the j-th shortest switching time interval as an example: then the unmanned aerial vehicle will not carry out impulse compensation.
(2) Suppose that Furthermore, if the impulse error ΔI zi satisfies the following condition: and the actual constant thrust of the unmanned aerial vehicle can be calculated as follows: 6 Discrete Dynamics in Nature and Society then the unmanned aerial vehicle will not carry out impulse compensation.
(3) If the impulse error ΔI zi satisfies the following condition: then the unmanned aerial vehicle should carry out impulse compensation, and the size of the constant thrust impulse compensation can be calculated as follows: Case 5. If the theoretical working time in the i-th thrust arc t * � m 1 ΔT, 1 < M 1 < M i and t * can be any M 1 shortest switching time interval in the i-th thrust arc; without loss of generality, we suppose that t * is the first m 1 shortest switching time interval, and the impulse error in the i-th thrust arc ΔI i can be calculated as follows: Furthermore, if there exist m 1 shortest switching time intervals satisfying the following conditions, without loss of generality, we suppose that these time intervals are the first m 1 shortest switching time interval, taking the j-th shortest switching time interval as an example: en, the size of the impulse compensation can be calculated as follows: (1) If the impulse error ΔI i satisfies the following condition: and the actual constant thrust of the unmanned aerial vehicle can be calculated as follows, taking the j-th shortest switching time interval as an example: Discrete Dynamics in Nature and Society 7 then the unmanned aerial vehicle will not carry out impulse compensation.
(2) Suppose that Furthermore, if the impulse error ΔI i satisfies the following condition: and the actual constant thrust of the unmanned aerial vehicle can be calculated as follows: then the unmanned aerial vehicle will not carry out impulse compensation. (3) If the impulse error ΔI i satisfies the following condition: then the unmanned aerial vehicle should carry out impulse compensation, and the size of the constant thrust impulse compensation can be calculated as follows: (48)

Compare Fuel Consumption and Design Switch Control
Laws. e fuel consumption under the theoretical continuous thrust and under the actual constant thrust is compared from the perspective of impulse compensation. We have already calculated the different impulse compensation according to the different conditions in Section 4, then the fuel savings under the actual constant thrust can be calculated as follows. Without loss of generality, taking the fuel savings in Case 5 as an example, we suppose that the mass flow rate of the propellant of the unmanned aerial vehicle's thruster is assumed to be p 0 g/s. Since the impulse error ΔI zi is as follows: 8 Discrete Dynamics in Nature and Society then the fuel savings in the i-th thrust arc ΔP i can be calculated as follows: ere are three types of time intervals in each thrust arc: the accelerating time intervals, the zero-thrust time intervals, and the decelerating time intervals. e task of the unmanned aerial vehicle collision avoidance maneuver is converted into the calculation of the number and sequence of three types of time intervals, respectively. In this section, the fuel consumption under the theoretical continuous thrust and the actual constant thrust is calculated and compared by using the method proposed in this paper. At last, the actual constant thrust switch control laws are obtained through the isochronous interpolation method, without loss of generality and taking Case 5 as an example. If the impulse error ΔI i satisfies the following condition: If the impulse error ΔI i satisfies the following condition: If the impulse error ΔI i satisfies the following condition: At last, the switch control laws for the collision avoidance maneuver can be given. For convenience, let us take the time intervals in the i-th thrust arc for an example: Discrete Dynamics in Nature and Society

Simulations
In order to present some of the results of the algorithms presented, we evaluated them in some nontrivial simulation scenarios. Suppose that ξ � (x, y, z) is anunmanned aerial vehicle position and obstacle centers are ξ 0 � (x 0 , y 0 , z 0 ), with axis lengths a, b, and c and index parametersd, e, and f. en, we can construct the function as follows [17]: where parameters a, b, c, d, e, and f determine the shape and size of the obstacle: if a � b � c and d � e � f � 1, the obstacle is a sphere; if a � b and d � e � 1, 0 < f < 1, the obstacle is regarded as a cone; if a and b are variables meeting the condition a � b � R 1 + (R 2 − R 1 )z/c and d � e � 1, f > 1 holds, the obstacle is a circular truncated cone approximately, where R 1 and R 2 are the radii of the two bases. Suppose that there are seven obstacles in flight space as shown in Table 1 and ρ 0 (ξ) � 1 and σ 0 (ξ) � 1. e starting position of the unmanned aerial vehicle is [0, 0, 0.5], and the target position is [40,40,6].
It can be seen from Table 1 that there are seven obstacles in the environment, of which four obstacles can be simplified as spheres, two obstacles can be simplified as cones, and one obstacle can be simplified as a cylinder.
It can be seen from Figure 1 that there are two planned paths for the unmanned aerial vehicle. e blue color curve is the planned path of the original algorithm, and the path length is 65.98288. e red color curve is the planned path of the improved algorithm, and the path length is 64.38021. Both planned paths of the two algorithms can avoid all obstacles and reach the target position smoothly. For the obstacle 3, the planned path by the original algorithm is to pass from the side of obstacle 3, but the planned path by the improved algorithm is to fly over the obstacle 3; this is because the tangential vector t(ξ) of the improved algorithm is parallel to the plane o − xy, so the length of the planned path is shorter.
It can be seen from Figure 2 that there are two flying height curves for the unmanned aerial vehicle. e blue color curve is the flying height curve of the original algorithm, and the red color curve is the flying height curve of the improved algorithm. e flying height of the improved algorithm is higher than the original algorithm, and the highest point is 6.97 at (24.23, 15.97, 6.97). It shows that the improved algorithm has better climbing performance.  Table 2. ρ 0 (ξ) � 1 and σ 0 (ξ) � 1 for each obstacle.
It can be seen from Figure 3 that there are two flying height curves for an unmanned aerial vehicle. e blue color curve is the flying height curve of the original algorithm, and the red color curve is the flying height curve of the improved algorithm. e planned path by the original algorithm can be stuck in the trap area formed by obstacles 3 and 4, but the unmanned aerial vehicle can fly over the trap area based on the improved algorithm. is  The original method The improved method  is because the improved algorithm takes the unmanned aerial vehicle's climbing ability into account. When the unmanned aerial vehicle encounters a trap area, it does not need to leave the trap area horizontally. is avoids the selection of the virtual target and simplifies the calculation process, resulting in greater flexibility.

Comparison of Different Saturation λ.
In the process of unmanned aerial vehicle 3D path planning, with the change of saturation λ, different paths will be planned. erefore, we need to design a reasonable saturation and give an optimal three-dimensional path. e different saturation λ are given in Table 3.
It can be seen from Figure 4 that there are four planned paths according to different saturation λ, and all four paths can bypass obstacles to reach the destination. It can be found that each path facing obstacles is different in the climbing height and the deflection angle. Further, the number of iterations, the length, and the highest point of each planned path are calculated and given in Table 4.
Path 1 tends to bypass obstacles horizontally. e maximum height of the unmanned aerial vehicle is only 4.71, but the path length is the longest. As the threshold value gradually decreases, the lengths of path 2 and path 3 are gradually shortened, and the maximum height of the corresponding unmanned aerial vehicle gradually increases, reaching 5.41 and 6.43, respectively. In path 4, the threshold is set to 0, which means that path 4 will tend to the shortest distance indefinitely, avoiding all obstacles and reaching the destination, but the rate of the climb drops a lot. Figure 5 shows the flight height curves of the above four planned paths. For path 1, path 2, and path 3, the unmanned aerial vehicle reaches the highest flying height near the obstacle 3. For path 4, the unmanned aerial vehicle reaches the highest flying height near the obstacle 1. erefore, we need to find a suitable threshold to make the unmanned aerial vehicle not only to have the ability to climb, but also to make the path as short as possible to ensures efficiency and safety of the planned path. e results in Figure 6 show the constant thrust fitting of F. According to the proposed criterion of this article, the unmanned aerial vehicle should carry out impulse compensation and the size of the constant thrust impulse compensation is the same, but the time of the constant thrust impulse compensation is different. e switch control laws can be given according to the sizes and the directions of the thrust accelerations of the unmanned aerial vehicle. e switch control law is given as follows:

Conclusions
is paper deals with the three-dimensional path planning problem of the unmanned aerial vehicle under constant thrust. e tangential vector and the disturbance matrix of the artificial fluid method are improved. A novel algorithm of constant thrust fitting is proposed through the impulse compensation, and then the constant thrust switching control scheme based on the isochronous interpolation method is given. It is proved that the planned path can avoid all obstacles smoothly and swiftly and reach the destination eventually.
e simulation results show that the switch control laws can effectively guarantee the unmanned aerial vehicle moving along the planned path.

Data Availability
e data used to support the findings of this study are included within the article.

Conflicts of Interest
e authors declare that they have no conflicts of interest. unmanned aerial vehicle (UAV) on the sea using hybrid differential evolution and quantum-behaved particle swarm optimization," IEEE Transactions on Systems, Man, and Cybernetics: Systems, vol. 43, no. 6, pp. 1451-1465, 2013.