Smooth 3D Path Planning by Means of Multiobjective Optimization for Fixed-Wing UAVs

: Demand for 3D planning and guidance algorithms is increasing due, in part, to the increase in unmanned vehicle-based applications. Traditionally, two-dimensional (2D) trajectory planning algorithms address the problem by using the approach of maintaining a constant altitude. Addressing the problem of path planning in a three-dimensional (3D) space implies more complex scenarios where maintaining altitude is not a valid approach. The work presented here implements an architecture for the generation of 3D ﬂight paths for ﬁxed-wing unmanned aerial vehicles (UAVs). The aim is to determine the feasible ﬂight path by minimizing the turning effort, starting from a set of control points in 3D space, including the initial and ﬁnal point. The trajectory generated takes into account the rotation and elevation constraints of the UAV. From the deﬁned control points and the movement constraints of the UAV, a path is generated that combines the union of the control points by means of a set of rectilinear segments and spherical curves. However, this design methodology means that the problem does not have a single solution; in other words, there are inﬁnite solutions for the generation of the ﬁnal path. For this reason, a multiobjective optimization problem (MOP) is proposed with the aim of independently maximizing each of the turning radii of the path. Finally, to produce a complete results visualization of the MOP and the ﬁnal 3D trajectory, the architecture was implemented in a simulation with Matlab/Simulink/ﬂightGear.


Introduction
The latest technological and scientific advances in the field of mobile robotics have enabled the area of autonomous vehicles (AV) to become a reality applied to the civil and military sectors [1,2]. Consequently, this technological branch presents a constant and vertiginous development in different fields, among them the development of new navigation and guidance techniques. The constant evolution in this field responds to new challenges in real applications [3][4][5].
It is important to note that the most common problem when determining a possible feasible path is the consideration of the AV intrinsic constraints. Therefore, the non-inclusion of kinematic and/or dynamic constraints of the AV when addressing the path planning problem may lead to non-feasible solutions that, for example, make it impossible for the AV to satisfactorily follow a path. However, including in the design all the constraints of the AV in the calculation phase of the path planning can lead to very complex optimization problems without a single solution and with very high computational costs. This paper focuses on the generation of smooth paths for fixed wing unmanned aerial vehicles (UAVs) that have high flight capabilities and extended flight-time missions, even in low power propulsion situations [19]. For these reasons, fixed-wing UAVs are suitable for use in terrain mapping applications for later action by the security forces, and in search and rescue tasks, both for the detection of people and the provision of first aid.
Due to the non-holonomic constraints of fixed-wing UAVs, the aim is to create a smooth three-dimensional curve from an initial point to a goal point through a complex 3D space with or without obstacles. To achieve this goal, it is essential to define a feasible path that minimizes flight turning effort and distance traveled. The ordered set of waypoints that will be used to generate the path to follow is defined as the control points set. In general, the problem of path planning is defined in a space region denoted as χ and split as the tripla (χ f ree , χ init , χ goal ). The movement space is defined as (W = R N : N = {2, 3}, where the obstacle region is denoted by χ obs , so that χ/χ obs is an open set denoting the collision-free space χ f ree . The initial condition χ init and the final condition χ goal are elements of χ f ree . The set of control points that define the collision-free space is calculated using specific path planning methods based on continuous and discrete environment sampling. Some examples of these techniques are: the rapidly-exploring random tree (RRT) [20][21][22][23]; probabilistic road maps (PRM) [24][25][26][27][28]; heuristic planners (genetic algorithms-GA) [29,30]; swarm intelligence [31][32][33][34]; fuzzy logic [35,36]); Voronoi diagrams [37][38][39]; artificial potential [40][41][42][43]; and recursive rewarding modified adaptive cell decomposition (RR-MACD) [44].
All the techniques mentioned build piece-wise paths in 2D or 3D to address the standard problem of path planning. These methods may provide optimal or near-optimal paths; however, they cannot guarantee smoothness and continuity, which could make it difficult to guide the UAV through the paths generated. Moreover, these techniques do not directly incorporate the operational constraints of the UAV and the environment. Therefore, this paper proposes a methodology to define feasible and smooth UAV paths, including system operational kinematic constraints.
The ability of an UAV to fly from one position to another and the consequent definition of the mission to be performed remains a challenge that requires the application of increasingly sophisticated strategies. One of the fundamental constraints to be considered in mission planning is the ability to be positionable and sensorially oriented (on the UAV or the environment) throughout the duration of the mission. This sensory location enables the construction of maps of the environment, and also enables the UAV to estimate its own current position and complete its self-location on the map. Similarly, it is important to note that good positioning and sensory orientation can be achieved by making a path plan and tracking it across or beyond the sensory detection domain. It is important to mention that the definition of smooth paths is not a new subject; various approaches have been proposed for non-holonomic UAVs, such as Dubins [45,46] in which paths are defined by connecting lines and arc-circular segments. The disadvantage is the generation of discontinuities in the connection points between segments. Another useful methodology in the literature is Clothoid curves [47][48][49], the main advantage being increases in curvature as a function of the arc-length; meanwhile the disadvantage lies in the limitation of length. These approximate curves generate continuities [50] of the type C 1 (a continuous path C 1 preserves the continuity of the tangent vector, in addition to maintaining continuous speed) and have been used in applications on mobile vehicles and on UAVs with flight limitation at a constant altitude, and examples are mentioned in [51][52][53][54][55]. An intuitive approach that ensures C 2 continuity (a continuous trajectory C 2 preserves the second order differential values at each point of the trajectory, in addition to maintaining the continuity of the acceleration vector) in 3D planning focuses on curvature and torsion zero at the junction points, as proposed in [56]. This is an interesting proposal where the 3D curve is built from a 2D curve in a (x, y) plane; the resulting length of this curve is the main parameter to build the curve in the other (x, z) plane. Finally, the Bézier, B-spline [57,58] are easy to implement polynomial curves; however, none are suitable for planning since they are sensitive to control parameters and weights [59], and do not take into account the constraints of the vehicles on which they are applied, so those curves need to be optimized. Finally, these types of parametric curves have no physical meaning and the relationship between design parameters and system variables is not defined. The above-mentioned approximation methodologies generate two different phenomena called interpolation (generation of a curve which must pass through the control points) and approximation (generation of a curve that approximates the control points, but may only go through the start and finish rather than all of them) detailed in [60,61]. The work presented here explores both phenomena in-depth to respond to the definition of feasible trajectories.
The starting point of this work is the set of 3D collision-free points generated by the various planners that guarantee that the path departs from the init point and ends at the goal point, and avoids static and dynamic obstacles. From these collision-free points, an ordered set of straight lines is built that define a first path that will later be smoothed to incorporate the feasibility and constraints of the UAV. The UAV constraints in this work are focused on its ability to turn horizontally and vertically. Therefore, for a UAV to complete a sequence of turns at a defined speed, it must determine its minimum turning radius Rp. If the turning radius is too small, the UAV will lose the trajectory; however, if the turning radius grows, the UAV can perform maneuvers with less effort.
The aim is to maintain 3D planning results, and at the same time, generate a finite set of possible 3D curves that optimize an approximate 3D curve, and simultaneously, the turns of the UAV-which is raised as a multiobjective optimization problem (MOP) [62]. This approach will result in a set of paths that meet UAV constraints expressed as dominant solutions on a Pareto n dimensional front [63]. Finally, selection criteria must be applied to determine the desired response from the point of view of curvature κ and torsion τ of the generated 3D curve. Thus, in order to verify the functionality of the proposed methodology, the results of the curves generated after the 3D curve optimization were compared with a known Bézier type approximation methodology [64].
This document is structured as follows. A brief summary of MOP concepts is given in Section 2.1. In Section 2.2, a brief description of smooth curves is given. Section 3 presents the formulation of the problem. Section 4 details the complete methodology for solving the problem. Section 5 details the experiments and results of 3D smooth path planning. Finally, conclusions and future work are presented in Section 6.

MultiObjetive Optimization
The optimization problem (OP) attempts to determine a solution that represents the optimal value (minimum or maximum) of a function, such as f : X → R, where X is a feasible decision vector, being min( f (x)) : x ∈ X. However, for problems where simultaneous optimization of more than one objective is necessary, i.e., multiobjective optimization (MOP), the function is shaped f : x → R k , where k ≥ 2 is the number of objectives. Therefore, the value vector of the target function could be defined as f : However, there is not usually a single X that generates an optimum that simultaneously satisfies each of the k objectives, due to the conflict between the objectives. The aim is to find a situation in which all objectives are satisfactorily within acceptable parameters. The MOP solution leads to points where any improvement in one target results in the degradation of any other target (one or more). Thus, these points are represented as a Pareto front [63], where all the points of the front are equally optimal.
Therefore, as expressed in [62] the MOP can be established as subject to: where θ ∈ R n is the decision vector, D is the decision space; J(θ) ∈ R m is the target vector; g(θ) and h(θ) are constraint vectors; and finally, θ il is the upper boundary and θ iu is the lower boundary of the decision space. Consequently, there is no single optimal model; in fact, there is a set of optimal solutions with different trade-offs between objectives, where none is better than the others. Using the definition of dominance, the Pareto set Θ P is the set of each non-dominated solution.
In this way, the Pareto domination is defined in case a solution θ 1 dominates another solution θ 2 ; that is, ( where J i (θ), i ∈ B := [1 · · · m] are the objectives to be optimized. Therefore, the optimal set of Pareto Θ P is given by where Θ p and J(Θ p ) are MOP solutions. However, in most cases they are unreachable because Θ P normally includes infinite solutions. Therefore, a finite set of Θ * P from Θ P and another finite set of J(Θ * p ) from J(Θ p ) represent satisfactory solutions. Starting from J(Θ * p ), the decision maker selects a solution according to the established preferences. For example, a certain point in the Pareto front that is close to the ideal point (called utopia point) J ideal .
Hence, an appropriate methodology for characterizing MOP is known as the elitist multiobjective evolutionary algorithm ( −MOGA) [62], which makes a distributed approach to Pareto's front. The aim of −MOGA is to find an intelligent distributed convergence towards a set of -Pareto; i.e., determine θ * P along the Pareto front J(Θ P ). The target space is split into a fixed number of boxes. Therefore, for each dimension i ∈ B, cells n_box i wide are created i , where Each box can be occupied by a single solution; therefore, this grid produces an intelligent distribution and preserves the diversity of J(Θ * P ). In addition, it is important to note that only the occupied boxes are verified, avoiding the need to use other clustering techniques to obtain adequate distributions. On the other hand, for a solution θ ∈ D, box i (θ) is defined as Then, box(θ) = {box 1 (θ), · · · , box m (θ)}. A solution θ 1 with value J(θ 1 ) dominates the solution θ 2 with value J(θ 2 ), denoted by θ 1 ≺ θ 2 , only if box(θ 1 ) ≺ box(θ 2 ) ∨ box(θ 1 ) = box(θ 2 ) and θ 1 ≺ θ 2 .
Hence, it is possible to control the maximum number of solutions to characterize the Pareto front. Finally, due to the definition of the box, the anchor points J i (θ i * ) are assigned a value of box i (θ i * ) = 0, whereby J i (θ i * ) = J min i . Therefore, no solution θ can -dominate them because, by applying the definition of the box, their box i (θ) ≥ 1.
The above process delivers two defined sets of responses: (a) the Pareto front Θ * P deploys a finite set of minimum values as the optimal path response within the search space; (b) the corresponding optimal points J(Θ * j ).

3D Curves for UAVs
A non-holonomic UAV [65] can perform flights in 3D Euclidean space. Nevertheless, to complete each movement sequence (horizontal and vertical), a set of UAV flight constraints must overcome. An UAV attempts to perform 3D movements at a defined velocity, meaning that the UAV is moving continuously, attempting to maintain that velocity. However, an UAV has a maximum capacity of turn and elevation at a defined velocity. Therefore, the aim is to build a 3D smooth curve inside the UAV flight turning boundaries, in such a way to reach a complete 3D smooth curve tracking.
A smooth curve can be defined as the representation of a continuous function, which can be expressed as C C C : I → X, where I is the curve interval composed of real numbers, while X represents the topological space. If the topological space is three-dimensional X = R 3 , then C C C : [a, b] → R 3 , is a differentiable injectable and continuous function, whose arc-length s is independent of the parametrization C C C. If a UAV is considered as a particle that travels along the envelope of the C C C curve at a defined speed v, this particle suffers changes in its local coordinate system due to the set of rotations of the C C C curve. Therefore, C C C must not allow rotation changes outside of the intrinsic rotation capabilities of the UAV.
It is important to mention that the formulation known as Frenet-Serret [66,67] describes the kinematic properties of particles that move along three-dimensional Euclidean space R 3 continuous and differentiable C C C(s), and parametetrized by its arc-length s (the arc-length is an invariant Euclidean characteristic of the curve). However, if we assume a curve given by a series of points along the Euclidean space as r(t), where the parameter t does not need the arc-length, then the tangent (T), normal (N) and bi-normal (B) derived vectors can be described, as all are mutually perpendicular (orthogonal base). Therefore, according to the theory of differential geometry of curves [68], the following equals can be defined as: where r (t) = dr(t) dt , r (t) = d 2 r(t) dr 2 and r (t) = d 3 r(t) dr 3 are the position vector derivatives r(t). These three vectors configure a navigation reference system of the UAV. Similarly, it is important to mention that the tangent vector T(t) is parallel to velocity, while the normal vector N(t) is represented by the change of direction per time unit of velocity.
The curvature terms κ (change of direction of the vector tangent T(t) to the curve r(t)) and torsion τ (change of direction of the vector bi-normal B(t)) are defined as: κ indicates a direct correlation with the horizontal rotation capability of the UAV, while τ indicates the elevation capability of the UAV. Therefore, the triedro Frenet-Serret can be defined in matricial notation as a skew-symmetric matrix: where the point over the variable indicates the derivative with respect to the parameter of arc-length s.
In summary, the space curve according to the formulation Frenet-Serret, defines a smooth curve that will be built from a spatial point and its tangent vector; this curve will be generated in the Euclidean 3D space depending on the pre-defined values of κ and τ, and finalize at a spatial 3D point after completing the arc-length s. However, our goal was to start from a defined point p init , and build a curve that touches a target point p goal ; therefore, the values of κ and τ have to fit so that when complete, the arc-length s will touch p goal . That gives us an infinite number of possible values of κ and τ with which we could meet that goal. Even if the maximum and minimum values of κ and τ are bounded, the complexity of the problem is high. In [69] it is mentioned that the complexity in a 2D environment becomes NP-hard, and the need is demonstrated for path planning algorithms that generate short paths with bounded curvatures in complicated environments. The aim is to find for possible approximate values of κ and τ that do not exceed the turn capabilities of the UAV. In other words, starting from Figure 1, and assuming a configuration of the UAV as a triplet (ρ, κ, τ), where ρ = [P 1 , · · · , P 5 ] is a dimensional vector that specifies n collision-free points, κ and τ are a set of curvatures and torsion along the path. The smooth curve starts from P 1 and reaches out to P 5 , and approaches the remaining ρ without affecting them, with values of κ and τ within the boundaries established by the maneuverability capabilities of the UAV. In summary, the selection of ρ goal points which determine radii of the tangent curves to the ρ points, will be obtained by solving a multiobjective optimization problem (MOP). This MOP is stated in such a way that the values of all radii will be maximized simultaneously. Obviously, the optimizer handles these values, taking into account that they are in conflict (as radius of one of them is increased, consequently, the adjacent radius are reduced. Therefore, the MOP solver will try to find a trade-off solution that guarantees the best set of points ρ goal for all tangent curves between control points ρ. > Figure 1. Perspective of the 3D flight problem for a fixed wing UAV, where C G represents the position vector of the center of gravity of the UAV. The global coordinate system CS g is placed at the origin, the orientation of the local body coordinate system CS b expressed by Euler roll, pitch and yaw angles respectively, which have been defined by three unitary orthogonal vectors aligned with the three axles of the vehicle and centered at C G . Finally, the angular velocities along the local axis of the UAV X, Y and Z are represented by p, q and r, respectively. The set of collision-free points P i is represented by black dots; the blue line describes the discrete path built from 3D path planning; the red dotted line is the new smooth optimized path that the UAV could follow.

Problem Definition
Let us assume a workspace W = R 3 , where it is possible to define a set of static or dynamic obstacles, such as ground or aerial boxes of different dimensions and locations (see Figure 1). The in-flight UAV receives data from its control station regarding environmental conditions and performs the necessary calculations to determine the best smooth 3D trajectory. Relevant data include the set of ordered 3D flight waypoints that are collision-free ρ = [P 1 , · · · , P 5 ] in the environment. The intrinsic maneuverability capabilities are determined by a Rp turning radius (which determines the vertical and horizontal turning limitations) defined by their flight speed. The aim is to start from ρ init and reach ρ goal in such a way that the UAV approaches the direct trajectory marked by the ordered sequence ρ. Therefore, ρ = P i (x i , y i , z i ) where (i = 1, · · · , n) and n is the total set of collision-free spaces and can be expressed as a discrete interpolation sequence ρ = f (t i ) → R, where f (t i ) is a set of nodes in 3D space. Therefore, it is possible to establish a set of sub-intervals (n − 1) between i = 1 and i = n partitioned in [a, b], defined as: A linear union between pairs of points then results in L L L : [a, b] → (x, y, z). And can be expressed as a set of straight lines that mark a direct flight path L L L(t) split into (n − 1) piece-wise.
Therefore, L L L(t) is a linear interpolation function for the discrete sequence ρ = f (t i ). In the same way, between the ρ points, there is a subset of (n − 1) straight lines that join the init and the goal of the trajectory along the collision-free flight space.
However, a non-holonomic UAV cannot perform every type of maneuver defined by L L L(t). In general, it is desirable to perform maneuvers with a high turning radius. Therefore, the approach presented in this work builds a smooth trajectory from ρ, that attempts to avoid inappropriate maneuvers using low values of κ and τ included within the boundaries of the UAV flight turn, while simultaneously closing in on the trajectory L L L(t).
Let us assume, from Figure 1, that the blue line denoted as L L L(t) is the direct path between the collision-free points of the environment, and the red dotted line is a smooth 3D spatial curve defined as C C C(t). The construction of this 3D smooth curve C C C(t) is done by joining a set of segments that can be of two types: spherical curves S (defined from a sphere of radius Rp) or straight lines L. Thus, each S segment is defined by three continuous points of ρ; this segment S has two points of tangency, one for each pair of adjacent straight lines L L L(t) formed by the current set of three points ρ. Hence, each S segment may have infinite solutions, with each radius Rp resulting in different tangent points on the lines L L L(t). Therefore, for each S segment, an infinite set of spheres can be defined, which will be linked through the relevant L segments or another S segment. Obviously, this approach to the problem suggests the existence of infinite combinations for the S and L segments. The way to address this issue has been through the approach of an MOP.

Methodology
This section describes the proposed methodology for the generation of 3D smooth trajectories. The proposed method is split into two parts, first detailing how the S segments were obtained, and then describing the union with the L segments.

Definition of Spherical Segment
Let us assume that from the result of a path planning, ρ = [P 1 , · · · , P n ] is the set of collision-free points of the environment described in Figure 2 (red points). This set of points is defined as As indicated above, Figure 2b shows an osculating sphere (oS) [68] defined with a minimum turning radius value Rp, located between the set of the first 3 P i and tangential to the straight lines L L L(t) formed between the same set of P i . Therefore, taking into account the number of collision-free points ρ, the set of spheres is equal to G i : i = {1, · · · , n − 2}, as shown in Figure 2c (orthogonal view). Figure 2b shows the first G i : i = 1 located among the three first P i s. Therefore, it is possible to define a plane π i : i = 1 between the same points P i , which will have an angle in relation to the location of the current set of points P i , as can be seen in the Figure 3a,b. The importance of the definition of this plane is given by the fact that within it is contained the center of G i with radius Rp. In this way, there is a self-contained curve S i (as a series of points along the Euclidean space) on the surface of the sphere and tangent to L L L(t) with t 2 and t 3 in plane π i , as shown in Figure 3; hence, the S i curve segment (black line) is defined as: where, x 0 , y 0 and z 0 together represent the center of G i . The curve S i performs a horizontal and vertical path due to the angular ranges of ψ and ϕ, which implies variations in the values of κ and τ (these have a direct connection to Rp and the arc-length of S i ). Consequently, if the value of Rp grows, S i also grows, while κ and τ decrease.

Remark 1.
If the plane π i is parallel to the horizontal plane (x, y) of the environment, then τ = 0 because the movements of the UAV will be horizontal. In the same way, if π i is parallel to the vertical plane (x, z) of the environment, then κ = 0.
However, before applying Equation (17), it is necessary to determine the location of the points (x 0 , y 0 , z 0 ) so that G i is tangent at a point on its surface with L L L(t), as shown in Figure 2b on the points (t 2 and t 3 ). Nevertheless, it should be noted that there is an angle between each pair of L L L(t), and this leads to G i approaching or moving away from the lines and their tangent points. Therefore, the geometric analysis applied to arrive at an optimal solution is detailed below.
First, a vector direction in space can be defined as − → v = p − q : p ∧ q ∈ R 3 . Therefore, starting from the known data ρ = [P 1 , · · · , P n ], taking Figure 2 as an example, where it is assumed that the collision-free initial points are (P i : {i = 1, · · · , 3}), a first set of two vectors is defined as: Just like a perpendicular vector from − → u i to − → v i , denoted as − → η , the normal vector is defined as: Consequently, the parametric equation of the π i plane containing three points is defined as: In the same way, the Euclidean distance defined between two points p and q is given by Therefore, two distances can be defined as du j : p = P (i+1) , q = P (i) and dv j : p = P (i+1) , q = P (i+2) . Finally, the angle between − → v i and − → u i is defined by: Therefore, with Equation (22) and Rp known, the tangential points at the lines L L L(t) can be located at a distance defined as: Thereby, two spatial points defined as pUi i and pUg i located in the direction of the vector − → u i provide that Pi = P i+1 , Pg = P i and d(p, q) = du i ; thus In the same way, two points pVi i and pVg i can be defined in the vector direction − → v i , so long as Pi = P i+1 , Pg = P i+2 and d(p, q) = dv i , according to Equation (24). Hence, the perpendicular bisector of pUi i and pVi i on the plane π i determines the center of the sphere (x 0 , y 0 , z 0 ). Figure 4a shows the application of the Equations (18)- (24), which can be repeated throughout the successive collision-free points ρ (this first Algorithm 1 is summarized in pseudocode). Between the centers of the spheres (x 0 , y 0 , z 0 ) and the points of intersection with L L L(t) are the displacement angles ϕ and ψ of the segment S i , as can be seen in the Figure 3, where pUi i = t 2 and pVi i = t 3 . (22), the angle formed between the points of intersection on the sphere G i , seen from its center towards the vertical or horizontal component, does not exceed 90 • in any case; that is, (0 • < ϕ < 90 • ) and (0 • < ψ < 90 • ).

Remark 2. Regardless of the angle condition produced by the pair of straight lines L L L(t) denoted in the Equation
(pVi i , pVg i ) =intersectPoint.(P i+1 , P i+2 , σ i ); 10: if i==1 then 12: t k = P 1 ; 13: t k+1 = pUi i ; 14: else 15: t k = pUg j ; 16: t k+1 = pUi i ; 17: end if 18: if i==n-2 then 19: t k+2 = pVi i ; 20: t k+3 = P n ; 21: end if 22: k = k + 2; 23: end for Algorithm 1 summarizes the geometric procedure followed. In line 4, a loop is started that performs through all the collision-free points marked by rho. From line 5 to line 9, the necessary computations are made to determine the center of G i (line 10). The definitions of the intervals containing the S and L segments are in lines 12, 13, 15, 16, 19 and 20. The described process shows the geometric analysis for the location of the set of spheres G i defined with constant radius Rp, as can be seen in Figure 2c. In addition, there is a set of four segments S and another set of five segments L, being segments S-those comprised by the intervals [t 2 , t 3 ], [t 4 , t 5 ], [t 6 , t 7 ] and [t 8 , t 9 ], and the segment L included by the intervals [t 1 , t 2 ], [t 3 , t 4 ], [t 5 , t 6 ], [t 7 , t 8 ] and [t 9 , t 10 ]. The goal now is to increase the radius Rp in each segment, so that the values of κ and τ along the curve are minimized, and the solution is to increase the radius Rp i in each G i .
The solution adopted in this work is to move the intersection point of each sphere G i in the direction of the adjacent segment L L L(t). Consequently, G i : i = 1 approximates symmetrically to the intervals t 1 and t 4 , G i : i = 2 makes the corresponding approximation to the intervals t 3 and t 6 , etc. Therefore, in Figure 4b, the segments L L L(t) can be seen adjacently to G i : i = 1, denoted as [t 1 ≡ P 1 , t 2 ≡ pUi j , t 3 ≡ pVi j , t 4 ≡ pVg j ]. Thus, between the adjacent intervals [t 1 , t 2 ] a vector is defined − → u i = t 2 − t 1 , and associated with this vector is a spatial point p i defined by the parametric equation: where θ i defines the space point p i along − → u i and within the intervals [t 1 , t 2 ]. Therefore, the value of the distance σ i from P i+1 to p i is defined according to the equation on (21), being p = P i+1 and q = p i+2 . A symbolic space point is defined q i between the intervals [t 3 , at the same distance σ i . Then σ i also has the angle φ i , and according to the Equation (22) it is possible to define a new Rp i according to Equation (23), which would have a higher radius value. Finally, the perpendicular bisector of p i and q i on the plane π i determines the center of G i (x 0 , y 0 , z 0 ) (see Figure 4b). Therefore, as the center of G i is defined, the Equation (17) defines the segment S i -and over this segment we find lower values of κ and τ according to the Equations (12) and (13).

Multiobjective Problem Definition (MOP)
Given Equation (25), it is important to note that any value of θ i between 0 and 1, defines a space point between the interval [t 1 , t 2 ]. In the same way, it is important to mention that within the boundaries of θ i , there is an infinite number of space points with an infinite number of radius Rp i and its corresponding infinite number of G i , with which its corresponding segments S i , can be built. Therefore, to obtain an optimal solution, a multiobjective problem (MOP) is solved using evolutionary algorithms [70] based on the concept of -dominance [71]. To do this, it is necessary to define the decision variables, the initial conditions of the process, the constraints of the MOP and the index vector to be optimized to represent the Pareto front. If it is assumed that the number of spheres G i is equal to m, and the number of objectives for each G i is equal to two, then J ideal (θ) = [J 1 (θ), J 2 (θ), · · · , J 2 * m (θ)] is the objectives vector, where J i denotes the i th objective.
where J A i and J B i depend on the decision variables vector θ. Assume D as a decision space within a subset R D , where θ is the decision variable vector composed of a set of θ i for all i ∈ 1 < i < m, where θ i is [0, 1] D . Consequently, the MOP problem can be stated as: where from Equation (12) from Equation (23) In summary, the aim is to find an optimal 3D smooth curve that minimizes κ and τ in each of the possible S i . It is important to mention that the adjacent spheres G i can grow into each other, until a maximum of q i ∈ − → v i ≡ p i+1 ∈ − → u i , which implies a decrease in the total set of segments, as described Figure 2d, where the green dotted line shows the set of S i segments belonging to G i .
An example of reconstruction according to the response Θ * P can be seen in Figure 2d, where the S reconstruction is made in four segments, defined by the boundaries [t 2 , t 3 ], [t 4 , t 5 ], [t 5 , t 6 ] and [t 7 , t 8 ]; the S segments belonging to C C C(t) are defined according to Equation (17).
In contrast, and with reference to Figure 2d, the L segments are defined by the rest of the boundaries, those boundaries being [t 1 , t 2 ], [t 3 , t 4 ], [t 6 , t 7 ] and [t 8 , t 9 ].

Definition of Straight-Line Segment
A segment path of L in a straight-line can be described by two points in the Euclidean space. Figure 2d shows an example of an L segment defined by the [t 1 , t 2 ] points, where the direction of the line is given by the flight path of the UAV. Therefore, − → u ( Figure 5) is a unit vector that points in the direction of the desired orientation, and with d defined as the distance between t 1 and t 2 according to Equation (21). Therefore, the L segments will be described, in general, as: Finally, the interpolation of S and L build a final 3D smooth curve on the plane (x, y, z).

Experiments and Results
This section presents the results of the computer simulation using Matlab/Simulink software and flightGear flight simulator for graphical visualization.
In this section we analyze five scenarios in 3D space, similar to the methodology proposed in [44]. Recursive rewarding modified adaptive cell decomposition (RR-MACD) splits the 3D environment like a discrete mesh of collision-free voxels. In particular, it places the UAV within an initial voxel and determines which of the adjacent voxels is the most suitable to make a displacement. To determine the best choice of displacement, the set of adjacent voxels have a set of associated constraints (conditions such as distance, vertical, and horizontal movement angles, battery consumption, etc.) to be satisfied before performing a discrete displacement. The voxel that minimizes the total effort and satisfies the constraints will be the next collision-free point. The RR-MACD methodology gives two sets of results based on the defined constraints. The results presented in [44] are shown in summary form in Table 1, where the first column shows the scenario number. The second column shows RR-MACD with four constraints and the RR-MACD with 10 constraints in the third column shows the conditions to solve the path planning problem. The 3D control points reflected in Table 1, ρ x (F) ≈ ρ, are the starting points for analyzing the method described in this paper for generating 3D smooth curves. Finally, it is important to note that the algorithms have been run on an Intel(R) Core(TM) i7-4790 3.60 GHz CPU (Manufacturer: Gigabyte Technology Co., Ltd., Model: B85M-D3H) with 8Gb RAM and S.O. Ubuntu Linux 16.04 LTS. The algorithms were programmed in MATLAB version 9.4.0.813654 (R2018a). Table 1. 3D path planning results. The number of free-collision voxels within a defined environment is indicated as S f ree , while the number of discrete 3D nodes built by [44] is denoted by ρ x (F) (3D final path). It is important to mention that the characteristics of the UAV assumed in the experiments have been taken from [65], whose study has been carried out on a fixed wing UAV model kadett 2400, represented by six states (x, y, z, φ, θ, ψ), where the first three states define the position vector of the UAV's global coordinate system, located at the origin of its center of gravity. The last three are the Euler angles of roll, pitch and yaw respectively, which define the orientation of the UAV.

RR-MACD 10 Constraints
Finally, the aim of the simulations is for the UAV to maintain its continuous flight at a defined constant speed of 18 m/s, within an established minimum radius of curvature Rp = 33 m to achieve smooth behavior and without maneuvers that could endanger the integrity of the aircraft.
It is important to remember that because the number of ρ x (F) = [P 1 , · · · , P n ] collision-free points is greater than five, a proper visualization method is essential to the decision making process for the final solution. Thus, the graphical representation method called level diagram [72] has been used, which consists of representing each objective and each design parameter in separate diagrams, all of which are synchronized with their y axis. Synchronization is made with the normalized distance of each point from the Pareto front to the ideal point. Therefore, with a brief training this representation offers a good visual understanding of the compromise reached on the Pareto front.

Bézier
To compare the results, an approximation of curves has been applied using Bézier [68], which enables the generation of trajectories for non-holonomic systems through a set of discrete control points, where the curve in a multidimensional environment can be described as: Bernstein-Bézier generates a finite vector of points belonging to the curve and guarantees to go through the first and last control points (translated as ρ = [P 1 , · · · , P n ]) and remaining inside the convex envelope.
It is important to mention that the interpolation of segments S and L builds a set of 3D smooth curves, all of which are possible solutions. It is, therefore, necessary to address a decision stage (decision maker (DM)) that selects one of them; i.e., a point on the Pareto front. In this work, the selection criteria based on the shortest distance to the ideal point have been used. Figures 6 and 7 show the selected point in red from J(Θ * p ) and Θ * p , which have been selected using the ∞-norm standard. Figure 8a shows the build of the 3D smooth curve C C C(t), while Figure 8b shows the best optimization in terms of curvature κ and total torsion τ, according to Equations (12) and (13). It is important to note that the mathematical mean of the values of the geometric variables κ and τ tends to be low. However, in some particular cases an increase is detected due to the change in direction of the flight; i.e., when the UAV is terminating an S segment in one direction and another S segment begins in the opposite direction. The curve generated by Bézier B(t) is shown as a yellow line, a more direct route can be seen between the init and goal point. However, this curve is very close to the bottom obstacle. To solve this, different authors propose modifying the control points ρ, or adding new points within the points established initially, as remarked in Section 1. Figure 9 shows the set of four additional examples from Table 1, in order to represent the functionality of the algorithm. Similarly, it is important to note that the number of ρ was different in each experiment, as were the altitudes, which guaranteed movement and 3D planning. It is also important to mention that the first environment shown in Figure 9a has smaller flight dimensional characteristics, so the turning radius in this example was set at Rp = 3 m with an average flight speed of 1.7 m/s.  To describe the different groups of trajectories generated by L L L(t), C C C(t) or B(t) displayed in Figure 9, Table 2 shows the flight results from the init to goal point in terms of distances; according to one of the results of each environment set by Table 1. It can be seen that the set of greater distances corresponding to the path in the form of the straight line marked by L L L(t), C C C(t) reduces the distance by L L L(t). As B(t) makes an approximation (as a mathematical expression) between the set of ρ of each environment, its route is the shortest. The column "EAA Error (meters)" shows the approximate absolute error EAA = 1 n ∑ n i=1 |A − B|, where A = L L L(t) and B = C C C(t) ∧ B = B(t). Therefore, the results of the column EAA Error (meters) show a greater approximation by C C C(t) and which results in a better avoidance of obstacles. Table 2. Flight distance of the curves. The column "Flight distance (meters)" shows the distance in meters at the initial and final free collision points marked by ρ. The column "EAA Error (meters)" shows the average error in meters along the trajectories.

L(t) C C C(t) B(t) L(t) vs C C C(t) L(t) vs B(t) #1
182   Similarly, Table 3 shows a set of results referring to the five environments analyzed. The averages of κ and τ generated along each smooth curve shows that B(t) exceeds C C C(t). However, in the first environment there is a collision caused by the B(t) curve. Table 3. Average results of κ and τ along the curves C C C(t) and B(t). The column "Collision" indicates collision (x) or no collision (o) of a curve against an obstacle.

Conclusions and Future Works
This paper describes an approach to the generation of a continuous 3D smooth path that enables consideration of the operational constraints of fixed-wing UAVs.
Firstly, the document describes the formulation of the problem by defining two types of segments within the trajectory: the S segments as a set of sphere segments that ensure a continuous and minimum curvature profile, and then the definition of L segments that generally connect S.
Next, we proposed the resolution of an MOP problem to obtain the numerical values of the parameters of the trajectory, given that the problem has infinite feasible solutions. When solving an MOP problem, the DM stage is essential to finally select the desired point from the Pareto set of optimal solutions.
It is important to remember that with methods such as classic Bézier or B-splines curves, you can define the number of samples along the path. However, the distance measured between one point and the next is not the same or even close (the difference can be large). These types of curves are useful in relatively simple environments (few obstacles); however, as the number of obstacles grows, the control points increase due to trajectory planning. Consequently, the construction of the curve can cause collisions. This work enables a constant approach distance between pairs of contiguous points.
The kinematic constraints of UAVs have been considered in this work, in the same way that dynamic constraints could be calculated mathematically. However, an important consideration that can enhance the generation of new trajectories in a new job is to increase the optimization criterion by improving variables such as energy consumption or incomplete data in dynamic environments.
Connected with the results shown in this paper, several future works arise. For example, integration of dynamic obstacles (UAVs swarms or other aircraft systems) into the flying area. From the optimization point of view, the proposal can be improved by taking into account dynamic constraints (i.e, inertia, wind disturbances, torque forces, etc.) into the MOP problem. Similarly, new cost indexes as flight time or/and spent energy could be added to the optimization problem. Finally, implementation of the proposed algorithm under real conditions (UAV in an outdoor environment) and its application to different uses (such as satellite trajectory generation [73]) will be investigated.