Next Article in Journal
Study on Driver Cross-Subject Emotion Recognition Based on Raw Multi-Channels EEG Data
Previous Article in Journal
A Joint-Learning-Based Dynamic Graph Learning Framework for Structured Prediction
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

3D Smooth Trajectory Planning for UAVs under Navigation Relayed by Multiple Stations Using Bézier Curves

1
School of Automation, Beijing Institute of Technology, Beijing 100081, China
2
National Key Laboratory of Autonomous Intelligent Unmanned Systems, Beijing 100081, China
*
Author to whom correspondence should be addressed.
Electronics 2023, 12(11), 2358; https://doi.org/10.3390/electronics12112358
Submission received: 13 April 2023 / Revised: 19 May 2023 / Accepted: 22 May 2023 / Published: 23 May 2023

Abstract

:
Navigation relayed by multiple stations (NRMS) is a promising technique that can significantly extend the operational range of unmanned aerial vehicles (UAVs) and hence facilitate the execution of long-range tasks. However, NRMS employs multiple external stations in sequence to guide a UAV to its destination, introducing additional variables and constraints for UAV trajectory planning. This paper investigates the trajectory planning problem for a UAV under NRMS from its initial location to a pre-determined destination while maintaining a connection with one of the stations for safety reasons. Instead of line segments used in prior studies, a piecewise Bézier curve is applied to represent a smooth trajectory in three-dimensional (3D) continuous space, which brings both benefits and complexity. This problem is a bi-level optimisation problem consisting of upper-level station routing and lower-level UAV trajectory planning. A station sequence must be obtained first to construct a flight corridor for UAV trajectory planning while the planned trajectory evaluates it. To tackle this challenging bi-level optimisation problem, a novel efficient decoupling framework is proposed. First, the upper-level sub-problem is solved by leveraging techniques from graph theory to obtain an approximate station sequence. Then, an alternative minimisation-based algorithm is presented to address the non-linear and non-convex UAV trajectory planning sub-problem by optimising the spatial and temporal parameters of the piecewise Bézier curve iteratively. Computational experiments demonstrate the efficiency of the proposed decoupling framework and the quality of the obtained approximate station sequence. Additionally, the alternative minimisation-based algorithm is shown to outperform other non-linear optimisation methods in finding a better trajectory for the UAV within the given computational time.

1. Introduction

Recently, unmanned aerial vehicles (UAVs), or drones, have gained significant popularity for diverse applications, including aerial photography and videography, precision agriculture, infrastructure inspection, security and surveillance, and cargo delivery [1,2,3,4,5]. As the scope of applications extends and technology advances, the demand for UAVs is increasing rapidly [6]. However, most UAV operations are conducted within the visual line of sight (VLOS) range [7] for safety and airspace integration considerations, which greatly limits their full potential and necessitates beyond visual line of sight (BVLOS) flights [8]. For BVLOS, more than one ground control station (GCS) might be involved in the navigation of a single UAV, which is called navigation relayed by multiple stations (NRMS). NRMS enables a UAV to cover larger areas, which facilitates the execution of long-range tasks [9,10,11].
Despite the promising benefits, NRMS employs multiple external stations in sequence to guide a UAV to its destination, which raises an additional station routing problem. During the flight, the UAV needs to be connected with one of the stations in order to receive the uplink data (e.g., control commands or navigation information) or transmit the downlink data (e.g., flight status or sensor data) in real time [12,13]. This introduces new characteristics and constraints for UAV trajectory planning.
Researchers have investigated the communication-aware trajectory planning problem for cellular-connected UAVs, which is an application of interest in NRMS. Zhang et al. [14] first addressed a connectivity-constrained trajectory planning problem and optimised the horizontal line-segment path by assuming that UAVs fly at a fixed minimum height and a constant maximum speed. Two methods were proposed to find high-quality approximate trajectories using graph theory and convex optimisation. Rather than optimising the horizontal path, Fonseca et al. [15] presented a reinforcement learning-based approach to optimise the altitude of a UAV to improve connectivity performance. Zhang et al. [16] further studied a radio map-based three-dimensional (3D) path planning problem for cellular-connected UAVs and derived the optimal UAV path by solving an equivalent shortest path problem (SPP). In addition to the connectivity constraint, disconnectivity constraints and outage-related constraints have also been explored in [17,18,19,20,21]. However, all these studies focused mainly on enhancing the communication performance by optimising the UAV’s trajectory with some simplifications on the manoeuvrability, ignoring the flyability of the planned trajectory. Specifically, they degraded a 3D trajectory planning problem to a two-dimensional (2D) path planning problem by assuming a constant height and speed for the UAV. Moreover, they used a line-segment path with path discretisation or time discretisation, which may need further path smoothing.
To generate smooth trajectories for UAVs, piecewise polynomial curves [22,23,24,25,26,27,28] have been commonly used in trajectory planning due to their simplicity and ease of differentiation, which is useful for velocity and acceleration planning. When using piecewise polynomial curves, the trajectory, with spatial and temporal parameters, has high non-linearity and non-convexity. However, the temporal parameters could be optimised efficiently using convex optimisation when the timing is fixed. How to find the optimal time allocation is still an open problem. Gradient descent with finite difference was used in [22,23], which is computationally expensive, especially for a high-dimensional time allocation vector. To decrease the computational cost, researchers have used heuristics to estimate the segment times [24,25,26] or mapped the segment duration to a velocity field induced by the Euclidean signed distance field of the map to achieve better time allocation [27]. Instead of using a traditional monomial polynomial basis, Gao et al. [27] used a Bernstein polynomial basis to represent the piecewise trajectory, as the safety and dynamical feasibility constraints can be enforced on the entire trajectory directly. Wang et al. [28] conducted alternating minimisation between boundary conditions and time durations for trajectory pieces, which can generate an optimal trajectory with the best time allocation.
In this paper, the 3D smooth trajectory planning problem for a quadrotor under NRMS is studied. The UAV has a mission (e.g., parcel delivery) that requires flying from an initial location to a destination, guided by multiple stations along the trajectory. The novelties and contributions are summarised as follows:
  • To the best of our knowledge, this is the first attempt to address the 3D smooth trajectory planning problem for a UAV under NRMS. Instead of line segments used in prior studies, a piecewise Bézier curve is applied to represent a continuous and smooth UAV trajectory in 3D space, which transforms an infinite number of trajectory variables into limited ones without loss of smoothness. This problem involves two levels of variables: the upper-level discrete station sequence and the lower-level continuous spatial and temporal parameters of the UAV trajectory. It is formulated as a bi-level, mixed-variable, constrained optimisation problem with the objective function being the mission completion time, which is subject to various constraints related to the NRMS and the UAV’s manoeuvrability. This problem is difficult to solve due to the interdependence of the two levels, non-convexity, and the presence of mixed variables.
  • To tackle this challenging problem, an efficient decoupling framework is proposed, where a high-quality approximate station sequence is obtained by first leveraging techniques from graph theory and then performing UAV trajectory planning based on the obtained station sequence. Compared with the enumeration method, the proposed decoupling approach can obtain a feasible and high-quality solution in a short computational time and shows excellent scalability, making it suitable for large-scale scenarios.
  • For lower-level UAV trajectory planning, a feasible initial trajectory is carefully designed, which provides a good start for later iterations. As the sub-problem with the fixed time allocation is a second-order cone programming (SOCP) problem and can be solved by off-the-shelf solvers efficiently, alternative minimisation is conducted between the time durations for trajectory segments and the control points’ locations of the curves. To update the temporal parameters, the entire trajectory is divided into several segments by fixing the boundary conditions at the endpoints, and then the time durations can be optimised separately. Compared with other non-linear optimisation algorithms, the proposed alternative minimisation-based method can obtain a better trajectory within the given computational time.
The rest of this paper is organised as follows. In Section 2, the optimisation model is given. Then, a decoupling approach is described in Section 3. Section 4 first demonstrates the effectiveness of the decoupling approach and then presents the comparison results between the alternative minimisation-based trajectory planning method and other optimisation algorithms. Finally, Section 5 concludes this paper.

2. Problem Formulation

The symbols and their descriptions are presented in Table 1. As shown in Figure 1, this paper considers M stations, which are scattered in different locations, O m = [ a m , b m , h m ] ( h m 0 , m = 1 , , M ) , over a large area and equipped with omnidirectional sensors, whose navigation range can be represented as hemispheres with radii R = [ R 1 , R 2 , , R M ] . An UAV has a mission to fly from a designated initial location W 0 = [ x 0 , y 0 , z 0 ] to its destination W F = [ x F , y F , z F ] . As  W 0 and W F are located within the range of two different stations, the UAV needs to be guided by multiple stations sequentially along its trajectory.

2.1. Decision Variables

The trajectory planning problem for a UAV under NRMS involves two levels of variables, including the upper-level station sequence and the lower-level UAV trajectory.

2.1.1. Station Sequence

Suppose there are M available navigation stations numbered as 1 , , M , a navigation station sequence S consisting of N stations can be represented by an N dimensional vector
S = [ S 1 , , S N ] ,
where S i { 1 , , M } for i = 1 , , N . As N denotes the number of stations chosen from the total pool of M stations, we have 1 < N M . Note that the choice of S 1 or S N may not be unique, as the UAV might start or end within the navigation range of multiple stations.

2.1.2. UAV Trajectory

Given a navigation station sequence S = [ S 1 , , S N ] , a piecewise trajectory u ( t ) consisting of N segments { u 1 ( t ) , , u N ( t ) } is defined as follows:
u ( t ) = u 1 ( t ) , t [ t 0 , t 1 ] , u 1 ( t 0 ) = W 0 , u 1 ( t 1 ) = W 1 u 2 ( t ) , t [ t 1 , t 2 ] , u 2 ( t 1 ) = W 1 , u 2 ( t 2 ) = W 2 u N ( t ) , t [ t N 1 , t N ] , u N ( t N 1 ) = W N 1 , u N ( t N ) = W F
where u i ( t ) , i = 1 , , N is the ith trajectory segment with respect to the time instant t [ t i 1 , t i ] . Here, t i 1 and t i represent the times when the ith navigation station S i starts and ends, guiding the UAV, respectively. [ W 1 , , W N 1 ] denotes the locations of a sequence of handover waypoints at which the control of the UAV is passed from one station to another. To ensure these trajectory segments are connected, we have u i ( t i ) = u i + 1 ( t i ) = W i for i = 1 , , N 1 .
To generate a smooth trajectory for the UAV, a piecewise Bézier curve is used to represent the UAV trajectory, which transforms an infinite number of trajectory variables into limited ones without sacrificing smoothness. The general form of a Bézier curve C ( u ) is given as follows [27,29]:
C ( u ) = j = 0 n B n , j ( u ) · P j ,
where n is the order of the curve and P j ( j = 0 , , n ) represents the location of the ( j + 1 ) th control point. u is a parameter that varies within the range [ 0 , 1 ] . B n , j ( u ) is the jth-order Bernstein polynomial basis, defined as [27,29]
B n , j ( u ) = n j u j ( 1 u ) n j ,
where n j is a binomial coefficient.
In this paper, we consider the smoothness of acceleration, which requires that n 4 . For simplicity, all trajectory segments are represented by Bézier curves of the same degree. As the parameter u [ 0 , 1 ] , the ith Bézier curve-based trajectory segment with respect to the time instant can be expressed as follows:
u i ( t ) = j = 0 n B n , j t t i 1 T i · P i , j , t [ t i 1 , t i ] ,
where P i , j ( i = 1 , , N , j = 0 , , n ) denotes the location of the ( j + 1 ) th control point of the ith trajectory segment. T i ( i = 1 , , N ) is the time duration for the ith trajectory segment, and  T i = t i t i 1 . As  T i appears in the denominator term of Equation (5), we have T i R + .
Based on the above, a UAV trajectory can be represented by the time durations for trajectory segments T , given by an N-dimensional vector
T = [ T 1 , , T N ] ,
where T i R + for i = 1 , , N , and  N ( n + 1 ) control points of the trajectory P , given by
P = P i , 0 P 1 , 1 P 1 , n P 2 , 0 P 2 , 1 P 2 , n P N , 0 P N , 1 P N , n
where P i , j R 3 for i = 1 , , N , j = 0 , , n .
Compared with other polynomial curves, Bézier curves have several special properties [29], which are useful for enforcing various constraints in this paper:
Endpoint interpolation property: a Bézier curve always starts at its first control point and ends at its last control point.
Trajectory segments can be connected easily.
Convex hull property: a Bézier curve is entirely confined within the convex hull defined by its control points.
A trajectory segment is constrained within a station’s navigation range (i.e., a hemisphere, which is convex) as long as all its control points are located within this convex area.
Hodograph property: the derivative of a Bézier curve is also a Bézier curve.
The velocity and acceleration along the trajectory can also be easily constrained.
An example of a piecewise Bézier curve U ( t ) of three segments is illustrated in Figure 2. t 0 and t 3 denote the times when the UAV starts and ends its mission. t 1 and t 2 denote the time when the transfer of control happens. T 1 , T 2 , and T 3 are the time durations for these three segments. The first segment is entirely constrained within the convex hull defined by its control points { P 1 , 0 , P 1 , 1 , P 1 , 2 , P 1 , 3 } . The first and last control points (i.e., P 1 , 0 and P 1 , 3 ) are located at the initial location W 0 of the UAV and the first handover waypoint W 1 to ensure this segment passes through these two key waypoints. The first control point P 2 , 0 of the second segment has the same location as the last control point P 1 , 3 to make these two segments connected.

2.2. Constraints

For the trajectory planning problem, a number of constraints must be taken into account to ensure the smoothness, dynamical feasibility, and some problem-specific constraints. Thanks to the Hodograph property, the velocity u i ˙ ( t ) and the acceleration u i ¨ ( t ) are also piecewise Bézier curves. Specifically, the ith velocity segment can be written as
u i ˙ ( t ) = j = 0 n 1 B n 1 , j t t i 1 T i · Q i , j ,
where Q i , j is the ( j + 1 ) th control point’s location of the ith velocity Bézier curve and can be derived as
Q i , j = n ( P i , j + 1 P i , j ) T i ,
where i = 1 , , N and j = 0 , , n 1 .
Additionally, the ith acceleration segment can be written as
u i ¨ ( t ) = j = 0 n 2 B n 2 , j t t i 1 T i · K i , j ,
where K i , j is the ( j + 1 ) th control point’s location of the ith acceleration Bézier curve and can be derived as
K i , j = n ( n 1 ) ( P i , j + 2 2 P i , j + 1 + P i , j ) T i 2 ,
where i = 1 , , N and j = 0 , , n 2 .

2.2.1. Constraints on Initial and Final States

As the initial and final locations of the UAV are predetermined, the trajectory is subject to the initial and final state constraints, which can be expressed as
u ( 0 ) = W 0 , u ˙ ( 0 ) = [ 0 , 0 , 0 ] , u ¨ ( 0 ) = [ 0 , 0 , 0 ] ,
u ( t N ) = W F , u ˙ ( t N ) = [ 0 , 0 , 0 ] , u ¨ ( t N ) = [ 0 , 0 , 0 ] ,
where t N denotes the mission completion time.
Due to the endpoint interpolation property of Bézier curves, the constraints on the initial and final states can be enforced on the first and last control points of the piecewise trajectory, velocity, and acceleration Bézier curves. Equation (12) can be rewritten as
P 1 , 0 = W 0 , Q 1 , 0 = [ 0 , 0 , 0 ] , K 1 , 0 = [ 0 , 0 , 0 ] ,
P N , n = W F , Q N , n 1 = [ 0 , 0 , 0 ] , K N , n 2 = [ 0 , 0 , 0 ] .

2.2.2. Constraints on Continuity

To ensure the smoothness of the piecewise trajectory, all the derivatives at the handover waypoints must be continuous. At the handover time T i , the continuity constraint can be expressed as
u i ( T i ) = u i + 1 ( T i ) , u ˙ i ( T i ) = u ˙ i + 1 ( T i ) , u ¨ i ( T i ) = u ¨ i + 1 ( T i ) , i = 1 , , N 1 .
Additionally, due to the endpoint interpolation property of Bézier curves, Equation (14) can be rewritten as
P i , n = P i + 1 , 0 , Q i , n 1 = Q i + 1 , 0 , K i , n 2 = K i + 1 , 0 , i = 1 , , N 1 .

2.2.3. Constraints on Connectivity

For safety concerns, the UAV must maintain a connection with one of the stations to be monitored or controlled (if necessary), which means that each UAV trajectory segment must be constrained within the navigation range of the corresponding station. For the ith trajectory segment, the connectivity constraint can be expressed as
| | u ( t ) O S i | | 2 R S i , t [ t i 1 , t i ] , i = 1 , , N ,
where | | · | | 2 denotes the Euclidean norm, and O S i and R S i are the location and radius of the ith station, respectively.
Thanks to the convex hull property of Bézier curves, the connectivity constraint is satisfied as long as all the control points of the Bézier curves are located within the corresponding stations’ range. Equation (16) can be rewritten as
| | P i , j O S i | | 2 R S i , i = 1 , , N , j = 0 , , n .

2.2.4. Constraints on Altitude

For UAV operations, the minimum and maximum allowable altitudes are specified by aviation authorities to ensure safety and prevent collisions with other aircraft or objects on the ground. The altitude constraint can be expressed as
z min u z ( t ) z max , 0 t t N ,
where u z ( t ) is the z coordinate of u ( t ) , and z min and z max are the minimum and maximum allowable altitudes, respectively.
According to the convex hull property of Bézier curves, Equation (18) can be rewritten as
z min P i , j z z max , i = 1 , , N , j = 0 , , n ,
where P i , j z is the z coordinate of the control point P i , j .

2.2.5. Constraints on Dynamic Feasibility

Due to the limitations of UAVs’ manoeuvrability and regulations on flight operations, the velocity and acceleration of a UAV must be less than a given threshold. The dynamic feasibility constraints can be expressed as
| | u ˙ ( t ) | | 2 v max , | | u ¨ ( t ) | | 2 a max , 0 t t N ,
where v max and a max are the maximum velocity and acceleration, respectively.
In addition, thanks to the convex hull property, the dynamic feasibility constraints can be easily enforced on the velocity and acceleration control points. Equation (20) can be rewritten as
| | Q i , j | | 2 v max , i = 1 , , N , j = 0 , , n 1 ,
| | K i , j | | 2 a max , i = 1 , , N , j = 0 , , n 2 .

2.3. Optimisation Model

This paper considers a commonly used objective function, i.e., the mission completion time [14,16,17,18,19,20]. By collecting the above-mentioned constraints, the optimisation problem is formulated as
( P 1 ) min S , T , P f = i = 1 N T i s . t . ( 13 ) , ( 15 ) , ( 17 ) , ( 19 ) , ( 21 ) .
 Remark 1. 
This model involves two coupling levels of variables: the upper-level discrete station sequence and S the lower-level continuous spatial and temporal parameters of the UAV trajectory { T , P } . The station sequence represents a constrained flight corridor in Equation (17) for UAV trajectory planning. Moreover, the dimension of S = [ S 1 , , S N ] determines the dimensions of T = [ T 1 , , T N ] and P = [ P i , 0 , , P i , n ] , i = 1 , , N . The quality of the station sequence is evaluated by the planned trajectory. The formulated problem is difficult to solve due to the interdependence of the two levels, non-convexity, and the presence of mixed variables.

3. Algorithm Design

To tackle this challenging bi-level problem, a naive approach is enumeration (i.e., list all possible station sequences, conduct trajectory planning based on all the sequences, and choose the best trajectory as the final solution), which is simple and straightforward but computationally expensive.
To solve this problem efficiently, a decoupling approach is proposed, as shown in Figure 3. First, the upper-level station sequence planning problem is solved by leveraging techniques from graph theory, and an approximate station sequence S is obtained. Then, the trajectory { T , P } is initialised within the constrained flight corridor based on S and refined iteratively via alternative minimisation.

3.1. Graph Theory-Based Station Routing Algorithm

The upper-level station routing sub-problem aims to build a constrained flight corridor for the UAV in the 3D continuous space, which is composed of a sequence of station navigation ranges. Considering Constraints (17) and (19), the flight corridor S can be expressed as
S = i = 1 N S i ,
where S i denotes the navigation range of the ith station S i and is given by
S i = { ( x , y , z ) | | | ( x , y , z ) O i | | 2 R i , z min z z max , ( x , y , z ) R 3 } .
Given M stations located in O m ( a m , b m , h m ) , m = 1 , , M with radii R m , an undirected graph G = ( V , E ) is defined as a vertex set V = { V O 1 , , V O M , V W 0 , V W F } and an edge set E. To judge the adjacency of two stations, we have the following proposition.
 Proposition 1 
(The adjacency between two station vertexes). Two station vertexes V O m and V O m , m m are adjacent as long as the distance between their projections on the horizontal plane is greater than or equal to the sum of the radii of S m and S m on the horizontal plane z = z min .
 Proof. 
As shown in Figure 4, two station vertexes are adjacent, indicating that they can provide continuous guidance for the UAV. In other words, the intersection of their navigation range cannot be an empty set, i.e.,  S m S m . Denote the radii of S m and S m on the horizontal plane z = z min as r m and r m , respectively. According to geometric characteristics, we have r m = R m 2 ( z min h m ) 2 and r m = R m 2 ( z min h m ) 2 . It is easy to draw the conclusion that V O m and V O m , m m are adjacent as long as | | ( a m , b m ) ( a m , b m ) | | 2 r m + r m holds.    □
Then, the edge set E is given by
E = { ( V W 0 , V O m ) : | | W 0 O m | | 2 R m , m = 1 , , M } { ( V O m , V O m ) : | | ( a m , b m ) ( a m , b m ) | | 2 r m + r m , m , m = 1 , , M , m m } { ( V W F , V O m ) : | | W F O m | | 2 R m , m = 1 , , M } .
As the objective of (P1) is to minimise the mission completion time, the Euclidean distance between vertex locations is used to approximately evaluate the cost. The weight of each edge is given by
w ( V W 0 , V O m ) = 0 ,
w ( V O m , V O m ) = | | ( a m , b m ) ( a m , b m ) | | 2 ,
w ( V W F , V O m ) = 0 .
Based on the above, the upper-level station routing problem can be transformed into an SPP on the undirected graph and efficiently solved by Dijkstra’s algorithm [30]. An approximate station sequence is thereby obtained, as illustrated in Figure 5.

3.2. Alternative Minimisation-Based UAV Trajectory Planning

The lower-level sub-problem aims to generate a smooth UAV trajectory within the flight corridor while satisfying various constraints. Given the navigation station sequence S = [ S 1 , , S N ] , (P1) can be rewritten as
( P 2 ) min T , P f = j = 1 N T i s . t . ( 13 ) , ( 15 ) , ( 17 ) , ( 19 ) , ( 21 ) .
 Remark 2. 
In (P2), the flight corridor is built according to Constraint (17) and the dimension of the lower-level trajectory planning sub-problem is determined. This sub-problem includes two parts of variables: the time durations for trajectory segments T and the control points’ locations for Bézier curves P , which have high non-linearity and non-convexity coming from their constraints. Specifically, T appears in the denominator term of Constraints (13), (15), and (21).
Note that when T is fixed, (P2) only has convex constraints. Specifically, Constraints (13), (15), and (19) are linear equality constraints, and Constraints (17) and (21) are second-order cone constraints. (P2) with constant T can be reformulated as a constraint satisfaction problem:
( P 3 ) s . t . P ( 13 ) , ( 15 ) , ( 17 ) , ( 19 ) , ( 21 ) .
Alternative minimisation (AM) is a popular optimisation technique that alternately minimises the objective with respect to each of its variables while fixing the others and has been successfully applied to solve many problems, including trajectory planning [28], traffic control [31], radar beamforming [32], and image enhancement [33]. As  P can be efficiently optimised with a constant T in (P3), alternative minimisation is considered for use in solving (P2). To optimise T , boundary conditions of trajectory segments D are introduced as intermediate variables. By fixing the boundary conditions of trajectory segments, the time durations can be minimised separately while maintaining continuity of the piecewise trajectory. In this paper, boundary conditions refer to location, velocity, and acceleration and are defined as
D 0 = P 1 , 0 Q 1 , 0 K 1 , 0 , D i = P i , n Q i , n 1 K i , n 2 , i = 1 , , N .
Based on the above analysis, an AM-based UAV trajectory planning algorithm is proposed, as shown in Algorithm 1. AM starts with a carefully designed, feasible initial trajectory { T ( 0 ) , P ( 0 ) } and then optimises P and T alternatively. In this paper, we consider two stopping criteria: the maximum runtime r u n t i m e max and the minimum difference threshold of two consecutive trajectories on the objective function value (i.e., mission completion time) δ T . The iterative process terminates and returns to the planned trajectory when either criterion is satisfied.
Algorithm 1 Alternative Minimisation-Based UAV Trajectory Planning Algorithm.
 Input: 
O , R , [ S 1 , , S N ] , W 0 , W F , z min , z max , V max , a max .
 Output: 
{ T * , P * }
 1:
Set i t e r = 0
 2:
Generate a feasible initial trajectory { T ( 0 ) , P ( 0 ) } via Algorithm 2
 3:
while   r u n t i m e < r u n t i m e max or i = 1 N T i ( i t e r 1 ) i = 1 N T i ( i t e r ) > δ T  do
 4:
    i t e r i t e r + 1
 5:
   Update P ( i t e r ) arg min P J ( T ( i t e r 1 ) , P ) based on Equation (37)
 6:
   Derive D ( i t e r ) from { T ( i t e r 1 ) , P ( i t e r ) } based on Equation (29)
 7:
   Minimise T i ( i t e r ) , i = 1 , , N separately via Algorithm 3;
 8:
end while
 9:
{ T * , P * } { T ( i t e r ) , P ( i t e r ) }
 10:
return { T * , P * }
Algorithm 2 Initial Trajectory Generation.
 Input: 
O , R , [ S 1 , , S N ] , W 0 , W F , z min , z max .
 Output: 
{ T ( 0 ) , P ( 0 ) }
 1:
W * arg min f p a t h ( W ) based on Equation (30)
 2:
Derive P ( 0 ) from [ W 0 , W * , W F ] based on Equation (32)
 3:
Derive T i , min , i = 1 , , N separately based on Equation (35)
 4:
T ( 0 ) [ T 1 , min , , T N , min ]
 5:
return { T ( 0 ) , P ( 0 ) }
Algorithm 3 Binary Search.
 Input: 
T i
 Output: 
T i
 1:
Set T l b = 0 and T u b = T i
 2:
while  T u b T l b > δ   do
 3:
    T t e m p = T l b + ( T u b T l b ) / 2
 4:
   Derive [ P i , 0 , , P i , n ] according to Equation (29) with T t e m p , D i 1 , and  D i
 5:
   Check the feasibility of { T t e m p , [ P i , 0 , , P i , n ] } in Constraints (17), (19), and (21)
 6:
   if  T t e m p is feasible then
 7:
      T u b = T t e m p
 8:
   else
 9:
       T l b = T t e m p
 10:
   end if
 11:
end while
 12:
T i = T u b

3.2.1. Initial Trajectory Generation

To start alternative minimisation, an initial trajectory { T ( 0 ) , P ( 0 ) } needs to be designed first. We consider a simplified model of (P2) by setting the order of the Bézier curves to n = 1 . The simplified problem is a path planning problem where the path is represented by piecewise line segments. The problem can be reformulated as
( P 4 ) min [ W 1 , , W N 1 ] f p a t h = i = 1 N 1 | | W i W i 1 | | 2 + | | W F W N 1 | | 2
s . t . | | W i O i | | 2 R i , i = 1 , , N 1
| | W i O i + 1 | | 2 R i + 1 , i = 1 , , N 1
z min W i z z max , i = 1 , , N 1
where W i denotes the ith handover waypoint in Equation (2) and W i z represents the z coordinate of W i .
(P4) is a second-order cone programming (SOCP) problem in which a convex quadratic function is minimised. Constraints (30b) and (30c) are second-order cone constraints. Constraint (30d) is a set of linear inequality constraints that can be transformed into second-order cone constraints. SOCPs have been widely studied in the literature and can be solved in polynomial time by interior point methods [34,35]. In this paper, the MOSEK interior-point optimiser [36] is used to solve (P4) as it is considered one of the state-of-the-art solvers for a wide range of problems, especially conic optimisation problems.
Based on the shortest feasible path [ W 0 , W 1 * , , W N 1 * , W F ] in (P4), a conservative initial trajectory can be constructed by setting the velocity and the acceleration at all the handover waypoints to [ 0 , 0 , 0 ] . For the piecewise trajectory, we have
P 1 , 0 ( 0 ) = W 0 , P i , n ( 0 ) = W i * , P N , n ( 0 ) = W F ,
Q 1 , 0 ( 0 ) = [ 0 , 0 , 0 ] , Q i , n 1 ( 0 ) ( 0 ) = [ 0 , 0 , 0 ] , Q N , n 1 ( 0 ) = [ 0 , 0 , 0 ] ,
K 1 , 0 ( 0 ) = [ 0 , 0 , 0 ] , K i , n 2 ( 0 ) = [ 0 , 0 , 0 ] , K N , n 2 ( 0 ) = [ 0 , 0 , 0 ] ,
where P i , n ( 0 ) , Q i , n 1 ( 0 ) , and K i , n 2 ( 0 ) are the location, velocity, and acceleration at the ith optimal handover waypoint W i * of the initial trajectory, respectively.
According to the expressions of Q i , j (see Equation (9)) and K i , j (see Equation (11)), we further have
P i , 0 ( 0 ) = P i , 1 ( 0 ) = P i , 2 ( 0 ) = W i 1 * , P i , n ( 0 ) = P i , n 1 ( 0 ) = P i , n 2 ( 0 ) = W i * ,
where i = 1 , , N . W 0 * and W N * refers to W 0 and W F , respectively. To ensure Equation (32) always holds, the number of control points for each Bézier curve must not be less than 6, which means n 5 . In this way, trajectory segments are connected smoothly so that Constraints (13) and (15) are satisfied.
Note that the determined partial variables (see Equation (32)) of P meet Constraints (17) and (19), which only involve the variable P . Although Constraint (21) involves both T and P , it can be met by adjusting only the value of T while P remains unchanged. Therefore, we set n = 5 to reduce the number of variables in P while maintaining the feasibility of the initial trajectory. The control points’ locations on the initial trajectory are given by
P ( 0 ) = [ W i 1 * , W i 1 * , W i 1 * , W i * , W i * , W i * ] , i = 1 , , N .
With P ( 0 ) and undetermined T , all the other constraints are satisfied apart from Constraint (21). Substituting P ( 0 ) into Constraints (21a) and (21b), we have
n | | W i * W i 1 * | | 2 T i v max , i = 1 , , N ,
n ( n 1 ) | | W i * W i 1 * | | 2 T i 2 a max , i = 1 , , N .
The minimum time duration of the ith trajectory piece satisfying Constraint (21) can be derived as
T i , min = max { n | | W i * W i 1 * | | 2 v max Equation ( 34a ) , n ( n 1 ) | | W i * W i 1 * | | 2 a max Equation ( 34b ) } .
As the location, velocity, and acceleration at all the handover waypoints are fixed, the time durations for the initial trajectory segments are given by
T ( 0 ) = [ T 1 , min , , T N , min ] .
To summarise, the initial trajectory generation method is presented in Algorithm 2. At first, the shortest line-segment path [ W 0 , W 1 * , , W N 1 * , W F ] is obtained by solving (P4) in Equation (30). Then, the control points P ( 0 ) are positioned on these key waypoints to satisfy Constraints (13), (15), (17), and (19). Next, the minimum time durations satisfying Constraint (21) are assigned to the time allocation T ( 0 ) . By this means, a feasible initial trajectory { T ( 0 ) , P ( 0 ) } is constructed, which provides a good start for later iterations.
An example of the initial trajectory is shown in Figure 6. The trajectory is the shortest line-segment path in the spatial space. As the UAV in this paper is a quadcopter, which can hover in place and fly in any direction, it follows the shortest line-segment path by stopping at all the handover waypoints (i.e., W 1 * and W 2 * in the figure) to ensure continuity in the velocity and acceleration profiles. So, the norm of the velocity and acceleration vectors decreases to zero at these key waypoints.

3.2.2. Update Operation of Control Point Locations

Although (P3) only has convex constraints, it is only a constraint satisfaction problem returning, whether the current P is feasible or not. Thus, it cannot indicate the minimisation direction and is not helpful in alternative minimisation.
According to the previous subsection, we can always construct P ( 0 ) satisfying all the constraints apart from Constraint (21). A constraint violation degree function J to Constraint (21) is defined to be minimised. (P3) can be reformulated as
( P 5 ) min P J = max i = 1 , , N j = 1 , , n { n · | | P i , j + 1 P i , j | | 2 T i · v max Formula A , n ( n 1 ) | | P i , j + 2 2 P i , j + 1 + P i , j | | 2 T i 2 · a max Formula B } s . t . ( 13 ) , ( 15 ) , ( 17 ) , ( 19 )
where J ( P ) is the maximum value of Formula A and Formula B. Formula A expresses the proportion of the 2-norm of the velocity vector to the maximum velocity. Formula B expresses the proportion of the 2-norm of the acceleration vector to the maximum acceleration. When the trajectory is feasible, the maximum value of Formula A and Formula B is less than or equal to 1. Furthermore, a smaller value of J ( P ) indicates that it is far from reaching the limit.
(P5) is also an SOCP and can be efficiently solved by the MOSEK interior-point optimiser. By minimising J ( P ) , the optimised trajectory has greater space to reduce the time durations. The velocity and acceleration profiles before and after the update operation of the control points’ locations are shown in Figure 7. In the figure, the optimised trajectory is flatter, and the maximum velocity and the maximum acceleration are smaller than those of the old trajectory.

3.2.3. Update Operation of Time Allocation

To optimise the time allocation T for the piecewise trajectory, the boundary conditions D are derived first following Equation (29). By fixing the boundary conditions, the piecewise trajectory is decoupled into several separate segments, allowing for the optimisation of the time durations separately.
As the objective of (P2) is to minimise the sum of time durations for trajectory segments, a better solution T can be obtained by minimising the time duration for each segment while maintaining feasibility. Since the initial trajectory is feasible and the update operation of the control points’ locations will only result in an improved trajectory, the trajectory prior to the update operation of time allocation must also be feasible. For the ith trajectory piece, the current time duration T i is feasible. However, when T i approaches 0, it becomes infeasible as it violates the dynamic feasibility constraints. Therefore, a binary search can be used to obtain a smaller and feasible T i , as shown in Algorithm 3.
The search process starts by defining the search range as [ T l b , T u b ] , which is initially set to [ 0 , T i ] . Then, the feasibility of the middle value T t e m p within the range is checked. If T t e m p is feasible, the search continues with the left half of the range. Otherwise, the right half is considered. Finally, when the length of the search range reduces to a pre-defined threshold, δ , the search is complete, and the upper bound of the feasible search space is returned.
Figure 8 displays the velocity and acceleration profiles before and after the update of the time allocation. Note that the trajectories prior to and following the update have the same boundary conditions [ Q 1 0 , Q 2 0 , Q 3 0 , Q 3 n 1 ] and [ K 1 0 , K 2 0 , K 3 0 , K 3 n 2 ] . However, the new trajectory has smaller time durations for trajectory segments, which brings it closer to the thresholds of dynamic feasibility constraints.

4. Numerical Results

This section presents numerical results to showcase the effectiveness of the proposed method, which consists of two computational experiments. Firstly, the proposed decoupling approach is compared with the enumeration method to demonstrate its efficiency. Secondly, the superiority of the AM-based trajectory planning algorithm is demonstrated by comparing it with other non-linear optimisation methods.
The simulations were conducted using the MATLAB_R2021b environment on a PC with an Intel Xeon Silver 4114 Processor 2.20GHz CPU and 32GB RAM.

4.1. Test Cases

To evaluate the performance across varying scales, three scenarios with different numbers of stations were designed in the computational experiments. The horizontal locations of stations were normally randomly distributed in an l × l square region. In each scenario, ten pairs of the initial and final locations of the UAV were also normally randomly distributed in an l × l square region and within the navigation range of stations, as shown in Figure 9. This generated a total of 30 test cases, which are identified by the number of stations and their index in each scenario. For instance, the notation M100_4 refers to the fourth pair of initial and final locations in the scenario M = 100 . The parameter settings for each case are given in detail in Table 2.

4.2. Performance of the Decoupling Framework

In order to showcase the efficiency of the proposed decoupling framework, a comparison was made between the enumeration approach (EM) and the proposed decoupling approach (DC), with a maximum runtime of 1800 s. Both approaches employed the proposed AM-based trajectory planning algorithm to solve the lower-level sub-problem. Since different upper-level station sequences can result in different dimensions of the lower-level UAV trajectory planning problem, the second stopping criterion, i = 1 N T i ( i t e r 1 ) i = 1 N T i ( i t e r ) > δ T , is more suitable. For this computational experiment, δ T was set to 1, and r u n t i m e max was set to + . The small number δ of 10 6 was set for the binary search.
The comparison results are presented in Table 3. In the table, three indicators are considered. The number of station sequences represents the quantity of upper-level solutions explored in the experiment. The decoupling method performs UAV trajectory planning on only one station sequence, whilst the enumeration approach explores as many upper-level solutions as possible within the maximum runtime. The mission completion time is the objective function, and Rank_DC refers to the rank of the solution obtained by the decoupling method among all the solutions obtained by both methods. The runtime refers to the time cost of the algorithms. Both of these two time-related indicators are in seconds.
In the first scenario ( M = 30 ), the runtime of the enumeration method is less than the maximum runtime, indicating that it explored all the possible station sequences and obtained the optimal solution. The decoupling method obtained a relatively good approximate solution in comparison to the optimal solution obtained by the enumeration method, which can be seen from Rank_DC. In the second scenario ( M = 100 ), a greater number of available stations and increased adjacency complexity among stations resulted in the enumeration method hitting the maximum runtime limit and therefore not exploring all potential upper-level solutions. This scenario led to the decoupling method finding better solutions than the enumeration method in some cases. In the third scenario ( M = 300 ), the larger upper solution space led to the decoupling method outperforming the enumeration method in all test cases.
A visualisation of the comparison results is shown in Figure 10, which provides an intuitive understanding of the performance of the decoupling method. The grey circles depict all the solutions obtained through EM, with the best solution highlighted in black. The orange triangles represent the approximated solution obtained via DC. From the figure, we can see that the decoupling method is capable of obtaining a relatively good station sequence at a much lower time cost. As the upper-level solution space increases, the runtime of the enumeration method increases significantly, while the decoupling method’s runtime only increases slightly. This implies that the decoupling method is more efficient and scalable for tackling larger-scale problems.

4.3. Performance of the AM-Based Trajectory Planning Method

To demonstrate the effectiveness of the trajectory planning method based on AM, Figure 11 visualises the iterative refinement of segment times in Case M30_1. The trajectories obtained in different iterations are depicted in different colours. From the figure, we can see that the initial trajectory is feasible but conservative. The trajectory obtained in the first iteration has the same time allocation as the initial trajectory but is flatter. However, in the second iteration, the trajectory is more aggressive and completes the mission in a considerably shorter time. As the algorithm proceeds, the trajectory is optimised in terms of the mission completion time and proximity to the threshold of the dynamic feasibility constraints.
Further, to evaluate its efficiency, the proposed AM-based algorithm was compared with the following algorithms:
  • Backtracking gradient descent (BGD) [22,23], in which the directional derivative is computed numerically by using finite difference and the step size is determined by backtracking line search.
  • NLopt [37,38], which is a popular and widely used open-source library for non-linear optimisation, supports both derivative-based and derivative-free optimisation, along with local and global optimisation. With no derivative information, six derivative-free local optimisation algorithms are compared, namely COBYLA, BOBYQA, NEWUOA-BOUND, PRAXIS, NELDERMEAD, and SBPLX.
These comparative algorithms were used to optimise the time allocation, while the locations of control points were optimised by MOSEK. As there were always trajectories (e.g., the initial trajectory { T ( 0 ) , P ( 0 ) } in Section 3) satisfying all the other constraints apart from the dynamic feasibility constraints, we softened and penalised the dynamic feasibility constraints in the objective function to use these algorithms to solve constrained optimisation problems, as suggested in [39]. For a fair comparison, all the algorithms started from the same initial trajectory constructed in Section 3 and stopped when the runtime exceeded the maximum runtime r u n t i m e max . In this computational experiment, r u n t i m e max was set to 3 s. The second stopping criterion of the proposed AM-based algorithm was disabled by setting δ T = . In addition, the small number δ in binary search was set to 10 6 .
Table 4 presents the comparison results of various algorithms on 30 test cases based on the mission completion time, with the best outcomes highlighted in bold. The numbers of the variables T and P are listed for each case, and it is observed that most of the variables belong to P , which can be efficiently optimised by the MOSEK interior-point optimiser when T is fixed. Therefore, it is advantageous to conduct alternative minimisation between T and P instead of optimising these variables simultaneously to solve the non-linear and non-convex UAV trajectory planning sub-problem. The comparison results show that all algorithms obtained feasible solutions, mainly due to the careful design of the feasible initial trajectory. It is noteworthy that the proposed AM-based algorithm outperformed the comparative algorithms in all test cases by finding trajectories with the shortest mission completion time within a given computational time of 3 s. Although some comparative algorithms (e.g., COBYLA) showed similar performance to the proposed AM-based algorithm in small-scale cases (e.g., Case M30_1), as the problem dimension increased, the proposed algorithm exhibited significantly better performance, with a noticeable widening of the performance gap between it and the comparative algorithms. This indicates better scalability of the proposed algorithm.
Figure 12 depicts radar graphs to illustrate the performance gap more clearly. In the figure, different algorithms are represented by different colours. This figure shows that the proposed algorithm consistently achieves the shortest mission completion time across all cases and is competitively superior to the comparative algorithms in solving this trajectory planning problem. Taking Case M300_10 as an example, the planned results obtained by all the algorithms are depicted in Figure 13. It is easy to see that the trajectory planned by the proposed AM-based algorithm is shorter and more aggressive than other algorithms while maintaining its feasibility.
The efficiency of the proposed algorithm mainly benefits from the alternative minimisation framework. With fixed time allocation, the tolerance of the dynamic feasibility constraints is maximised, offering greater space to further reduce the time durations. With fixed boundary conditions, the time durations are minimised separately to approach the threshold of the dynamic feasibility constraints while maintaining the feasibility of the trajectory. By this means, the update operation of time allocation may provide a reasonable direction of iterations and a larger step size. In addition, the update operation of time allocation is computationally cheap and thus might have more iterations within the limited runtime.
The convergence curves of the proposed AM-based algorithm and backtracking gradient descent are shown in Figure 14. In each scenario, three cases are selected for illustration. The solutions obtained via the backtracking gradient descent method and the AM-based algorithm are represented by blue squares and red circles, respectively. As shown in the figure, the backtracking gradient descent method only has a few iterations, while the proposed AM-based algorithm iterates many more times. This is because the finite difference is used in the backtracking gradient descent method, which requires solving N SOCPs for a trajectory with N segments for each gradient evaluation. In the AM-based algorithm, the time durations are minimised efficiently via binary search, and only a few feasibility check operations are required in each bisection.

5. Conclusions

This paper addressed the 3D smooth trajectory planning problem for a quadrotor that had a mission to fly from one point to another. As the distance between the UAV’s initial and final locations was larger than the navigation range of a single GCS, multiple stations were employed to guide the UAV in sequence to its destination, i.e., the NRMS technique. Given multiple available stations, there was an additional station routing problem, which was to find the optimal station sequence to build a flight corridor where the UAV could maintain a continuous connection with at least one of the stations for safety. Due to the interdependence of station routing and UAV trajectory planning, it is computationally expensive to optimise the station sequence and the UAV’s trajectory jointly. Therefore, an efficient decoupling approach was proposed. For the upper-level sub-problem, an undirected graph was first constructed, taking into account the adjacency among stations. Based on the graph, the station routing problem was transformed into an SPP and solved by the classic Dijkstra’s algorithm. For the lower-level sub-problem, alternative minimisation was conducted between the time durations for trajectory segments and the control points’ locations of curves to obtain a feasible and high-quality solution. Comparison results demonstrated the efficiency of the decoupling framework, the quality of the approximate station sequence, and the superiority of the AM-based trajectory planning algorithm.
Compared with prior studies in which trajectories were discretised into line segments, a piecewise Bézier curve was used to represent the 3D smooth trajectory that could be followed by the UAV, eliminating the need for additional path smoothing. In addition, the numerical results showed that the decoupling method could find a relatively good approximate solution in a short computational time, making it suitable for online trajectory planning when combined with receding horizon planning. Considering the fact that the use of BVLOS UAVs is heavily regulated at present due to safety concerns, the proposed approach can provide a solution to ensure a UAV is always being monitored and is under control while executing its mission, which is helpful to promote the widespread applications of BVLOS UAVs.
In future work, it is worth considering building a system consisting of multiple stations and a UAV for the implementation of the proposed approach in a real-world setting and validating its effectiveness. Moreover, other objectives (e.g., energy efficiency [40]) can be integrated into the trajectory planning model to enhance its applicability and versatility in various scenarios. A piecewise Bézier curve is still a good choice for trajectory representation, while some other methods (e.g., intelligent optimisation methods) could be considered for solving more complex models.

Author Contributions

Conceptualisation, M.Q.; methodology, B.X.; software, M.Q.; validation, M.Q.; formal analysis, M.Q.; investigation, M.Q.; resources, M.Q.; data curation, M.Q.; writing—original draft preparation, M.Q.; writing—review and editing, M.Q.; visualisation, M.Q.; supervision, L.D.; funding acquisition, L.D. All authors have read and agreed to the published version of the manuscript.

Funding

This work was supported in part by the National Outstanding Youth Talent Support Program 61822304, NSFC, the Basic Science Centre Programs of the NSFC under Grant 62088101, the Joint Funds of the National Natural Science Foundation of China under Grant U20B2073, the Beijing Advanced Innovation Centre for Intelligent Robots and Systems, the Shanghai Municipal Science and Technology Major Project (2021SHZDZX0100), and the Shanghai Municipal Commission of Science and Technology Project (19511132101).

Data Availability Statement

The data presented in this study are available in the article.

Conflicts of Interest

The authors declare no conflict of interest.

Abbreviations

The following abbreviations are used in this manuscript:
UAVUnmanned aerial vehicle
VLOSVisual line of sight
BVLOSBeyond visual line of sight
GCSGround control station
NRMSNavigation relayed by multiple stations
3DThree-dimensional
SPPShortest path problem
2DTwo-dimensional
SOCPSecond-order cone programming
AMAlternative minimisation
EMEnumeration approach
DCDecoupling approach
BGDBacktracking gradient descent

References

  1. Valavanis, K.P.; Vachtsevanos, G.J. Handbook of Unmanned Aerial Vehicles; Springer: Dordrecht, The Netherlands, 2015. [Google Scholar]
  2. Zhang, H.; Xin, B.; Dou, L.; Chen, J.; Hirota, K. A Review of Cooperative Path Planning of an Unmanned Aerial Vehicle Group. Front. Inf. Technol. Electron. Eng. 2020, 21, 1671–1694. [Google Scholar] [CrossRef]
  3. Khan, S.; Tufail, M.; Khan, M.T.; Khan, Z.A.; Iqbal, J.; Wasim, A. A Novel Framework for Multiple Ground Target Detection, Recognition and Inspection in Precision Agriculture Applications Using a UAV. Unmanned Syst. 2022, 10, 45–56. [Google Scholar] [CrossRef]
  4. Feng, K.; Li, W.; Han, J.; Pan, F. Low-Latency Aerial Images Object Detection for UAV. Unmanned Syst. 2022, 10, 57–67. [Google Scholar] [CrossRef]
  5. Zhang, R.; Dou, L.; Xin, B.; Chen, C.; Deng, F.; Chen, J. A Review on the Truck and Drone Cooperative Delivery Problem. Unmanned Syst. 2023, 12, 1–25. [Google Scholar] [CrossRef]
  6. UAV Market by Point of Sale, Systems, Platform, Function, End Use, Application, Type, Mode of Operation, Mtow, Range and Region—Global Forecast to 2027. Available online: https://www.marketsandmarkets.com/Market-Reports/unmanned-aerial-vehicles-uav-market-662.html (accessed on 30 September 2022).
  7. Operation and Certification of Small Unmanned Aircraft Systems. Available online: https://www.federalregister.gov/documents/2016/06/28/2016-15079/operation-and-certification-of-small-unmanned-aircraft-systems (accessed on 28 June 2016).
  8. Unmanned Aircraft System Operations in UK Airspace—Guidance. Available online: https://publicapps.caa.co.uk/modalapplication.aspx?appid=11&mode=detail&id=415 (accessed on 22 December 2022).
  9. Qi, M.; Dou, L.; Xin, B.; Chen, J. Optimal Path Planning for an Unmanned Aerial Vehicle under Navigation Relayed by Multiple Stations to Intercept a Moving Target. In Proceedings of the 2017 IEEE 56th Annual Conference on Decision and Control, Melbourne, Australia, 12–15 December 2017; pp. 2744–2749. [Google Scholar]
  10. Qi, M.; Dou, L.; Xin, B.; Chen, J. A Branch-and-Bound Approach for Path Planning of Vehicles under Navigation Relayed by Multiple Stations. In Proceedings of the 2017 11th Asian Control Conference, Gold Coast, Australia, 17–20 December 2017; pp. 168–173. [Google Scholar]
  11. Qi, M.; Dou, L.; Xin, B.; Chen, J. Optimal Path Planning for Vehicles under Navigation Relayed by Multiple Stations. In Proceedings of the 2017 IEEE International Conference on Systems, Man, and Cybernetics, Banff, AB, Canada, 5–8 October 2017; pp. 1465–1470. [Google Scholar]
  12. Zeng, Y.; Lyu, J.; Zhang, R. Cellular-Connected UAV: Potential, Challenges, and Promising Technologies. IEEE Wirel. Commun. 2018, 26, 120–127. [Google Scholar] [CrossRef]
  13. Song, Q.; Zeng, Y.; Xu, J.; Jin, S. A Survey of Prototype and Experiment for UAV Communications. Sci. China Inf. Sci. 2021, 64, 140301. [Google Scholar] [CrossRef]
  14. Zhang, S.; Zeng, Y.; Zhang, R. Cellular-Enabled UAV Communication: A Connectivity-Constrained Trajectory Optimization Perspective. IEEE Trans. Commun. 2019, 67, 2580–2604. [Google Scholar] [CrossRef]
  15. Fonseca, E.; Galkin, B.; Amer, R.; DaSilva, L.A.; Dusparic, I. Adaptive Height Optimization for Cellular-Connected UAVs: A Deep Reinforcement Learning Approach. IEEE Access 2023, 11, 5966–5980. [Google Scholar] [CrossRef]
  16. Zhang, S.; Zhang, R. Radio Map-Based 3D Path Planning for Cellular-Connected UAV. IEEE Trans. Wirel. Commun. 2021, 20, 1975–1989. [Google Scholar] [CrossRef]
  17. Bulut, E.; Guevenc, I. Trajectory Optimization for Cellular-Connected UAVs with Disconnectivity Constraint. In Proceedings of the 2018 IEEE International Conference on Communications Workshops, Kansas City, MO, USA, 20–24 May 2018; pp. 1–6. [Google Scholar]
  18. Zhang, S.; Zhang, R. Trajectory Optimization for Cellular-Connected UAV Under Outage Duration Constraint. J. Commun. Inf. Netw. 2019, 4, 55–71. [Google Scholar] [CrossRef]
  19. Zeng, Y.; Xu, X.; Jin, S.; Zhang, R. Simultaneous Navigation and Radio Mapping for Cellular-Connected UAV With Deep Reinforcement Learning. IEEE Trans. Wirel. Commun. 2021, 20, 4205–4220. [Google Scholar] [CrossRef]
  20. Chen, Y.J.; Huang, D.Y. Joint Trajectory Design and BS Association for Cellular-Connected UAV: An Imitation-Augmented Deep Reinforcement Learning Approach. IEEE Internet Things J. 2022, 9, 2843–2858. [Google Scholar] [CrossRef]
  21. Fontanesi, G.; Zhu, A.; Arvaneh, M.; Ahmadi, H. A Transfer Learning Approach for UAV Path Design With Connectivity Outage Constraint. IEEE Internet Things J. 2023, 10, 4998–5012. [Google Scholar] [CrossRef]
  22. Mellinger, D.; Kumar, V. Minimum Snap Trajectory Generation and Control for Quadrotors. In Proceedings of the 2011 IEEE International Conference on Robotics and Automation, Shanghai, China, 9–13 May 2011; pp. 2520–2525. [Google Scholar]
  23. Richter, C.; Bry, A.; Roy, N. Polynomial Trajectory Planning for Aggressive Quadrotor Flight in Dense Indoor Environments. In Robotics Research; Springer: Cham, Switzerland, 2016; pp. 649–666. [Google Scholar]
  24. Shomin, M.; Hollis, R. Fast, Dynamic Trajectory Planning for a Dynamically Stable Mobile Robot. In Proceedings of the 2014 IEEE/RSJ International Conference on Intelligent Robots and Systems, Chicago, IL, USA, 14–18 September 2014; pp. 3636–3641. [Google Scholar]
  25. Oleynikova, H.; Burri, M.; Taylor, Z.; Nieto, J.; Siegwart, R.; Galceran, E. Continuous-Time Trajectory Optimization for Online UAV Replanning. In Proceedings of the 2016 IEEE/RSJ International Conference on Intelligent Robots and Systems, Daejeon, Republic of Korea, 9–14 October 2016; pp. 5332–5339. [Google Scholar]
  26. Liu, S.; Watterson, M.; Mohta, K.; Sun, K.; Bhattacharya, S.; Taylor, C.J.; Kumar, V. Planning Dynamically Feasible Trajectories for Quadrotors Using Safe Flight Corridors in 3-D Complex Environments. IEEE Robot. Autom. Lett. 2017, 2, 1688–1695. [Google Scholar] [CrossRef]
  27. Gao, F.; Wu, W.; Lin, Y.; Shen, S. Online Safe Trajectory Generation for Quadrotors Using Fast Marching Method and Bernstein Basis Polynomial. In Proceedings of the 2018 IEEE International Conference on Robotics and Automation, Brisbane, Australia, 21–25 May 2018; pp. 344–351. [Google Scholar]
  28. Wang, Z.; Zhou, X.; Xu, C.; Chu, J.; Gao, F. Alternating Minimization Based Trajectory Generation for Quadrotor Aggressive Flight. IEEE Robot. Autom. Lett. 2020, 5, 4836–4843. [Google Scholar] [CrossRef]
  29. Gao, F.; Wu, W.; Gao, W.; Shen, S. Flying on Point Clouds: Online Trajectory Generation and Autonomous Navigation for Quadrotors in Cluttered Environments. J. Field Robot. 2019, 36, 710–733. [Google Scholar] [CrossRef]
  30. West, D.B. Introduction to Graph Theory; Prentice Hall: Upper Saddle River, NJ, USA, 2001. [Google Scholar]
  31. Todorović, U.; Frejo, J.R.D.; De Schutter, B. Distributed MPC for Large Freeway Networks Using Alternating Optimization. IEEE Trans. Intell. Transp. Syst. 2022, 23, 1875–1884. [Google Scholar] [CrossRef]
  32. Qi, C.; Ci, W.; Zhang, J.; You, X. Hybrid Beamforming for Millimeter Wave MIMO Integrated Sensing and Communications. IEEE Commun. Lett. 2022, 26, 1136–1140. [Google Scholar] [CrossRef]
  33. Zhuang, P.; Wu, J.; Porikli, F.; Li, C. Underwater Image Enhancement With Hyper-Laplacian Reflectance Priors. IEEE Trans. Image Process. 2022, 31, 5442–5455. [Google Scholar] [CrossRef]
  34. Alizadeh, F.; Goldfarb, D. Second-Order Cone Programming. Math. Program. 2003, 95, 3–51. [Google Scholar] [CrossRef]
  35. Sheng, M.; Zhao, C.; Liu, J.; Teng, W.; Dai, Y.; Li, J. Energy-Efficient Trajectory Planning and Resource Allocation in UAV Communication Networks under Imperfect Channel Prediction. Sci. China Inf. Sci. 2022, 65, 222301. [Google Scholar] [CrossRef]
  36. MOSEK: A Software Package for Solving Large-Scale Optimization Problems. Available online: https://www.mosek.com (accessed on 27 March 2023).
  37. The NLopt Nonlinear-Optimization Package. Available online: http://github.com/stevengj/nlopt (accessed on 4 December 2021).
  38. Ji, Y.; Wang, M.; Yu, Y.; Zhu, C. Curvature-Based r-Adaptive Isogeometric Analysis with Injectivity-Preserving Multi-Sided Domain Parameterization. J. Syst. Sci. Complex. 2023, 36, 53–76. [Google Scholar] [CrossRef]
  39. Burri, M.; Oleynikova, H.; Achtelik, M.W.; Siegwart, R. Real-Time Visual-Inertial Mapping, Re-localization and Planning Onboard MAVs in Unknown Environments. In Proceedings of the 2015 IEEE/RSJ International Conference on Intelligent Robots and Systems, Hamburg, Germany, 28 September–3 October 2015; pp. 1872–1878. [Google Scholar]
  40. Carabin, G.; Scalera, L. On the trajectory planning for energy efficiency in industrial robotic systems. Robotics 2020, 9, 89. [Google Scholar] [CrossRef]
Figure 1. Example of a UAV under NRMS from its initial location to a pre-determined destination while maintaining a connection with one of the stations.
Figure 1. Example of a UAV under NRMS from its initial location to a pre-determined destination while maintaining a connection with one of the stations.
Electronics 12 02358 g001
Figure 2. Example of a piecewise Bézier curve of degree n = 3 consisting of three segments.
Figure 2. Example of a piecewise Bézier curve of degree n = 3 consisting of three segments.
Electronics 12 02358 g002
Figure 3. Flowchart of the proposed decoupling framework, where the steps for the upper-level and lower-level sub-problems are marked in blue and orange, respectively.
Figure 3. Flowchart of the proposed decoupling framework, where the steps for the upper-level and lower-level sub-problems are marked in blue and orange, respectively.
Electronics 12 02358 g003
Figure 4. Illustration of the judgement condition for adjacency between two stations.
Figure 4. Illustration of the judgement condition for adjacency between two stations.
Electronics 12 02358 g004
Figure 5. Illustration of the upper-level station sequence planning process.
Figure 5. Illustration of the upper-level station sequence planning process.
Electronics 12 02358 g005
Figure 6. Illustration of the initial trajectory with its velocity and acceleration profiles.
Figure 6. Illustration of the initial trajectory with its velocity and acceleration profiles.
Electronics 12 02358 g006
Figure 7. Velocity and acceleration profiles before and after updating the control points’ locations.
Figure 7. Velocity and acceleration profiles before and after updating the control points’ locations.
Electronics 12 02358 g007
Figure 8. Velocity and acceleration profiles before and after updating the time allocation.
Figure 8. Velocity and acceleration profiles before and after updating the time allocation.
Electronics 12 02358 g008
Figure 9. Three scenarios with different numbers of stations and randomly generated initial and final locations for UAV trajectory planning.
Figure 9. Three scenarios with different numbers of stations and randomly generated initial and final locations for UAV trajectory planning.
Electronics 12 02358 g009
Figure 10. Performance comparison of the enumeration method (EM) and the proposed decoupling method (DC) on mission completion time (i.e., the upper part) and runtime (i.e., the lower part) in seconds.
Figure 10. Performance comparison of the enumeration method (EM) and the proposed decoupling method (DC) on mission completion time (i.e., the upper part) and runtime (i.e., the lower part) in seconds.
Electronics 12 02358 g010
Figure 11. Illustration of the iterative process of the proposed AM-based trajectory planning method. (a) Trajectories generated at different iterations; (b) Altitude, velocity, and acceleration profiles versus time at different iterations, where the dashed lines at the end of the curves indicated the mission completion time of the trajectories.
Figure 11. Illustration of the iterative process of the proposed AM-based trajectory planning method. (a) Trajectories generated at different iterations; (b) Altitude, velocity, and acceleration profiles versus time at different iterations, where the dashed lines at the end of the curves indicated the mission completion time of the trajectories.
Electronics 12 02358 g011
Figure 12. Radar graphs depicting the performance gap among eight algorithms on mission completion time in three scenarios.
Figure 12. Radar graphs depicting the performance gap among eight algorithms on mission completion time in three scenarios.
Electronics 12 02358 g012
Figure 13. Planned results of eight algorithms for lower-level UAV trajectory planning in Case M300_10. (a) Upper-level station sequence obtained by Dijkstra’s algorithm; (b) Lower-level trajectories obtained by eight algorithms; (c) Altitude, velocity, and acceleration profiles versus time obtained by eight algorithms, where the dashed lines at the end of the curves indicated the mission completion time of the trajectories.
Figure 13. Planned results of eight algorithms for lower-level UAV trajectory planning in Case M300_10. (a) Upper-level station sequence obtained by Dijkstra’s algorithm; (b) Lower-level trajectories obtained by eight algorithms; (c) Altitude, velocity, and acceleration profiles versus time obtained by eight algorithms, where the dashed lines at the end of the curves indicated the mission completion time of the trajectories.
Electronics 12 02358 g013
Figure 14. Convergence curves of backtracking gradient descent (BGD) and the proposed alternative minimisation-based algorithm (AM) on mission completion time in seconds in selected cases.
Figure 14. Convergence curves of backtracking gradient descent (BGD) and the proposed alternative minimisation-based algorithm (AM) on mission completion time in seconds in selected cases.
Electronics 12 02358 g014
Table 1. Symbols and descriptions.
Table 1. Symbols and descriptions.
SymbolDescription
StationMThe number of available stations, M Z , M > 1 ;
O m The location of the mth station, O m R 3 , m = 1 , , M ;
R m Radius of the mth station navigation range, R m R + , m = 1 , , M ;
NThe number of chosen stations, N Z , 1 < N M ;
S i The ith chosen navigation station, S i { 1 , 2 , , M } , i = 1 , , N ;
UAV u i ( t ) The ith trajectory segment with respect to the time instant, t [ t i 1 , t i ] ;
W 0 The initial location, W 0 R 3 ;
W F The final location, W F R 3 ;
W i The location of the ith handover waypoint, W i R 3 , i = 1 , , N 1 ;
z min The minimum altitude, z min R ;
z max The maximum altitude, z max R ;
v max The maximum velocity, v max R ;
a max The maximum acceleration, a max R ;
T i The time duration of the ith trajectory segment, T i R + , i = 1 , , N ;
nThe degree of Bézier curves, n Z + ;
P i , j The ( j + 1 ) th control point’s location of the ith trajectory Bézier curve,
P i , j R 3 , i = 1 , , N , j = 0 , , n ;
Q i , j The ( j + 1 ) th control point’s location of the ith velocity Bézier curve,
Q i , j R 3 , i = 1 , , N , j = 0 , , n 1 ;
K i , j The ( j + 1 ) th control point’s location of the ith acceleration Bézier curve,
K i , j R 3 , i = 1 , , N , j = 0 , , n 2 ;
Table 2. Case parameter settings in computational experiments.
Table 2. Case parameter settings in computational experiments.
ParameterValue
Scenario 1l10,000 m
M30
Scenario 2l20,000 m
M100
Scenario 3l30,000 m
M300
Stations a m , b m , m = 1 , , M U ( 0 , l ) m
h m , m = 1 , , M 0 m
R m , m = 1 , , M [500, 1500] m
UAV x 0 , y 0 , x F , y F U ( 0 , l ) m
z 0 , z F U ( z min , z max ) m
z min 100 m
z max 120 m
V max 50 m/s
a max 5 m/s 2
Table 3. Comparison results of the enumeration method (EM) and the proposed decoupling method (DC) on 30 test cases.
Table 3. Comparison results of the enumeration method (EM) and the proposed decoupling method (DC) on 30 test cases.
CaseNumber of
Station Sequences
Mission Completion Time (s)Runtime (s)
No.NameEMDCEMDCRank_DCEMDC
1M30_1421158.54168.66654.710.86
2M30_2511206.72225.31766.240.97
3M30_3691125.04140.251275.700.87
4M30_42071146.52156.082327.140.98
5M30_52281165.47175.938346.031.07
6M30_62761122.92127.316420.300.81
7M30_72881183.72189.2510477.081.10
8M30_83421148.43149.533523.760.91
9M30_95311151.50154.464918.040.98
10M30_108281201.00201.7041624.501.25
11M100_16791252.12256.9741800.161.38
12M100_27701309.11303.6911800.331.41
13M100_36271304.34304.3411801.121.42
14M100_46371391.84398.1131800.391.69
15M100_56541370.72372.7431802.021.43
16M100_63961563.69319.1611800.951.41
17M100_75131330.95332.5621801.371.55
18M100_85211487.02412.4411802.101.85
19M100_94611604.76273.5411800.531.42
20M100_104231660.74357.3911802.511.86
21M300_18161294.95267.6711802.701.31
22M300_25691472.93315.3011802.661.39
23M300_33811619.67390.3711801.001.64
24M300_44891659.90417.5611802.781.89
25M300_530911223.22503.5011804.812.02
26M300_63551807.23535.9011805.962.09
27M300_75601679.89547.3311801.312.26
28M300_83101908.12602.0511801.712.64
29M300_93851813.51688.4511801.472.54
30M300_1026211547.86695.2411801.442.68
The best results are highlighted in bold.
Table 4. Comparison results of eight algorithms for the lower-level UAV trajectory planning sub-problem on mission completion time in seconds.
Table 4. Comparison results of eight algorithms for the lower-level UAV trajectory planning sub-problem on mission completion time in seconds.
CaseNumber of VariablesCOBYLA
+ MOSEK
BOBYQA
+ MOSEK
NEWUOA-BO-
UND + MOSEK
PRAXIS
+ MOSEK
NELDERMEAD
+ MOSEK
SBPLX
+ MOSEK
BGD
+ MOSEK
AM
+ MOSEK
No.NameTPTotal
1M30_16108114180.71426.15394.94389.51360.55221.10513.20165.67
2M30_27126133409.82653.60418.18534.35610.90431.79659.03222.40
3M30_359095197.60319.05235.76179.40354.61174.76299.13135.56
4M30_46108114187.37422.99292.29211.91404.65299.13370.17151.61
5M30_58144152315.91585.62557.96347.18488.86329.86470.21172.48
6M30_66108114181.10337.84214.63224.53287.98152.31332.85123.94
7M30_78144152348.26655.56425.63485.03505.60397.71588.32185.24
8M30_87126133302.11459.13339.42309.44354.21352.46436.84144.26
9M30_96108114262.14380.60261.33172.90344.91214.56295.97149.89
10M30_108144152320.57631.50428.32457.60599.41362.97618.42196.11
11M100_18144152464.76935.87884.65773.75761.79553.95898.22253.65
12M100_210180190765.061163.71784.37999.211053.56919.301057.17299.14
13M100_310180190811.341167.38950.00998.681043.37908.021101.79300.50
14M100_413234247938.861561.071264.711403.641590.941484.901373.60394.54
15M100_511198209798.251402.801070.931018.471284.131070.691098.14368.47
16M100_612216228653.961170.53752.10958.041185.311054.76958.08314.99
17M100_713234247740.891306.021360.801150.711218.171059.641124.03327.98
18M100_8162883041059.501695.871448.791569.561785.901449.001416.82408.15
19M100_910180190649.71983.17878.66946.25878.25784.521002.30269.49
20M100_1014252266672.021405.671235.691347.411365.031293.451263.75352.58
21M300_17126133525.00876.17574.49318.79784.77592.28390.84262.49
22M300_29162171769.151081.121260.70834.04895.95908.78982.40310.98
23M300_311198209808.851463.941230.12823.631427.721255.92822.98385.72
24M300_413234247848.051688.571496.501599.491596.181449.611597.61412.89
25M300_5152702851260.661968.581481.721658.342054.451725.701465.91499.14
26M300_6173063231282.662206.581525.611773.882376.842041.901414.90531.25
27M300_7193423611331.532381.232638.532153.832505.642237.461920.59543.26
28M300_8203603801756.992607.252907.242477.682758.982451.522207.49598.70
29M300_9223964181855.982667.602760.892538.782802.322530.352172.03684.11
30M300_10234144371837.073113.593370.892701.313264.962996.731941.38693.12
The best results are highlighted in bold.
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Qi, M.; Dou, L.; Xin, B. 3D Smooth Trajectory Planning for UAVs under Navigation Relayed by Multiple Stations Using Bézier Curves. Electronics 2023, 12, 2358. https://doi.org/10.3390/electronics12112358

AMA Style

Qi M, Dou L, Xin B. 3D Smooth Trajectory Planning for UAVs under Navigation Relayed by Multiple Stations Using Bézier Curves. Electronics. 2023; 12(11):2358. https://doi.org/10.3390/electronics12112358

Chicago/Turabian Style

Qi, Mingfeng, Lihua Dou, and Bin Xin. 2023. "3D Smooth Trajectory Planning for UAVs under Navigation Relayed by Multiple Stations Using Bézier Curves" Electronics 12, no. 11: 2358. https://doi.org/10.3390/electronics12112358

APA Style

Qi, M., Dou, L., & Xin, B. (2023). 3D Smooth Trajectory Planning for UAVs under Navigation Relayed by Multiple Stations Using Bézier Curves. Electronics, 12(11), 2358. https://doi.org/10.3390/electronics12112358

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop