Unmanned aerial vehicles (UAVs) are making increasingly long flights today with significantly longer mission times. This requires the UAVs to have long *endurance* as well as have long range capabilities. Motivated by locomotory patterns in birds and marine animals which demonstrate a powered-coasting-powered periodic locomotory behavior, an optimal control problem is formulated to study UAV trajectory planning. The concept of differential flatness is used to reformulate the optimal control problem as a nonlinear programing problem where the flat outputs are parameterized using Fourier series. The Π test is also used to verify the existence of a periodic solution which outperforms the steady-state motion. An example of an Aerosonde UAV is used to illustrate the improvement in endurance and range costs of the periodic control solutions relative to the equilibrium flight.

## Introduction

Unmanned aerial vehicles (UAVs) are increasingly being integrated into more and more areas of application. For example, Horcher and Visser in their article [1] motivate the need for UAVs in natural resource management, monitoring forestry, and agriculture where UAVs can serve as ideal vehicles to provide low cost observations with high temporal rates. In another article, Adams and Friedland [2] survey the role of UAVs in emergency response, emergency management, and postdisaster assessment. In many of these applications and others, UAVs are required to stay in the air longer or cover more distance to complete assigned duties. This requires the UAVs to have long endurances and long ranges. Therefore, there is a need to optimize these two performance criteria. Generally, aircraft trajectories are steady-state trajectories due to the constraint of human comfort as there is a limit to the amount of acceleration and jerk that is considered comfortable for humans. Even though the human factor is not present in a UAV, the steady-state trajectory is still used. Many researchers have already shown that the steady-state optimal solutions (to maximize endurance or range) for aircraft flight can be nonoptimal [3–5] and seeking periodic solutions could benefit performance.

Motivation to look for periodic solutions can be dominantly found in nature. Several processes like the orbital motion of planets, the rhythm of the beating heart, or even the flight of some migratory birds are known to be periodic [6]. Gleiss et al. [7] even conjecture that there is an evolutionary convergence in locomotion pattern in aerial and marine animals which is characterized by a periodic interspersed active propulsion and coasting phases for optimal locomotion.

Nature has also shown that periodic locomotion is more optimal than steady-state in terms of energetics. For instance, the albatross, which covers about 100,000 miles annually in flight, exhibits periodicity while taking into account the wind shear [8]. The zebra finch is another bird that presents periodic locomotion. Periodic locomotion is not just limited to birds: mammals such as the dolphin and the elephant seal also move in similar patterns. Numerous other animals that present intermittent locomotion as opposed to steady motion can be found in article [9] by Paoletti and Mahadevan. Moreover, articles [6,10,11] have also found that periodic motion is more efficient than pure steady motion. Since higher efficiencies in energy consumption of UAVs in flight can potentially increase endurance and range, seeking periodic control trajectories for aircraft path planning finds strong motivation.

Optimal periodic solutions to improve performance have been found in various processes such as drug delivery [12], vehicle cruise control [13], methane production from micro-algae [14], and harvesting [15,16] to name a few. Various researchers have tackled the problem of optimal trajectory generation for aircrafts as well. For example, pioneering fuel optimal trajectories related work can be found in Refs. [17–21]. Several among them have also recognized the benefit of periodic control in improving fuel optimal performance. For example, Speyer [4] and Dewell and Speyer [19] showed the nonoptimality of the steady-state flight and the improvement in optimal periodic flight in a hypersonic aircraft using the energy model. In recent years, Harada and Bollino [20] and Chen and Speyer [22] have worked on the fuel optimal periodic solutions for UAVs.

In this work, the concept of differential flatness is exploited to solve an optimal control problem to determine the desired periodic trajectories of the control inputs. Differentially flat systems are the ones where a set of outputs can be identified, which equal the number of inputs, such that all the system states and inputs can be algebraically represented in terms of the flat outputs and their derivatives. Differential flatness, first introduced by Fliess et al. [23], has been extensively used in trajectory planning of UAVs [24–26]. It has been coupled with the linear quadratic regulator-based controller to give optimal paths [27], used to control multiple UAVs [28], and even applied to obstacle avoidance trajectory generation [29–31]. It has also been used in planar vertical takeoff and landing vehicles [32], in trucks and trailers [33,34], wheeled mobile robots [35], oscillation damping [36], control of fuel-cells [37] as well as periodic drug delivery [12]. In fact, a noncomprehensive catalog of mechanical systems which are differentially flat can be found in the article by Murray et al. [38]. This paper adds to the pool of knowledge by employing the concept of differential flatness to determine endurance and range optimal UAV periodic trajectories.

The paper is structured in the following way. Section 2 presents the dynamic model of the UAV used and the assumptions on which the formulation is based. Section 3 gives an overview of the concept of differential flatness. Section 4 summarizes the Π test. Section 5 shows the formulation of the optimization problem for both the steady-state solution and the periodic control solution. Section 6 presents the results of the Π test and the optimal control problem illustrated on the Aerosonde UAV. Finally, Sec. 7 provides concluding remarks.

## Unmanned Aerial Vehicles Model

A point-mass aircraft model is used to represent UAV dynamics. The model captures the dynamic effects encountered in civil aviation [39] and is widely accepted in the community. The model, however, makes the following assumptions:

The bank angle (

*ϕ*) and the angle of attack (*α*) can be changed instantaneously (i.e., the UAV actuator dynamics are relatively much faster).The lift and drag of the UAV changes instantaneously with respect to the angle of attack (

*α*).The mass of the UAV remains constant while seeking solutions to optimization problems, as the change in fuel weight in the period of interest is small relative to the gross weight of the UAV.

*x*,

*y*, and

*z*are the position of the aircraft center of gravity in the earth reference frame,

*γ*is the flight-path angle,

*χ*is the heading angle as illustrated in Fig. 1,

*V*is the aircraft speed,

*g*is the acceleration due to gravity, and

*m*is the mass of the aircraft.

*T*is the thrust,

*ϕ*is the bank angle,

*n*is the load factor, and

*D*is the drag. The load factor (

*n*) is given by

where *L* is the lift, *ρ* is the density of air, *S* is the planform area of the UAV wing, *α* is the angle of attack, and *C _{L}*

_{0}and

*C*are the aerodynamic coefficients of lift. The expression for lift can be found in Ref. [41].

_{Lα}*D*) is also taken from the literature [41] and is given by

*e*is the Oswald efficiency factor, AR is the wing aspect ratio, and

*C*

_{D}_{0}is the constant parasitic drag. Wing stall is not taken into consideration while modeling the aerodynamic forces on the aircraft. Consequently, a limit is placed on the angle of attack to ensure accuracy of the model. It is assumed that the position coordinates of the UAV are available for measurement. Hence, the inputs and outputs of the model can be conveniently grouped as

## Differential Flatness

*flat outputs*) and its derivatives. For a system that meets these conditions, the state and control variables can be expressed as algebraic functions of the flat outputs, while the flat outputs can be parametrized using basis functions in time. This allows us to seek solutions to constrained optimal control problems by solving a nonlinear optimization problem where the objective is to determine the coefficients of these basis functions under certain constraints. These constraints can be derived after mapping the state and the control constraints to the coefficient space. To illustrate this, consider the dynamic system

**is the output vector of interest. It should be pointed out that one of the conditions for a system to be differentially flat is that the number of the outputs should be equal to the number of inputs (therefore, we have $y\u2208Rm$). Under this condition, the system given by Eq. (10) is differentially flat if a flat-output vector $z\u2208Rm$ can be identified of the form**

*y*and (.)^{(}^{a}^{)} refers to the *a*th derivative of (.). The components of the flat outputs should also be differentially independent [42]. It should also be noted that the flat output ** z** need not be the output of interest

**. It is possible that the flat outputs do not have any physical meaning or significance. The flat output now allows one to plan optimal trajectories with boundaries and dynamic constraints in flat output space, which can then later be mapped back to the state space and the input space.**

*y**x*,

*y*, and

*z*serving as the flat outputs, permitting the states and inputs to be written as a function of the flat-outputs (

**) and their derivatives as follows:**

*Y*The flatness of the point-mass model can now be exploited to determine periodic optimal control trajectories by solving a nonlinear programing problem in the flat output space (by determining the coefficients of the basis functions used to parameterize the flat outputs).

## The Π Test

*τ*and an output function

where $G\u2032$ is the transpose of *G*, and $H\xafxx,\u2009H\xafux,\u2009H\xafxu,\u2009and\u2009H\xafuu$ are second-order partial derivatives of the Hamiltonian evaluated at the steady-state solutions. The Π test determines whether small sinusoidal perturbations of the input ** u** around the optimal steady solutions improves performance. Performance improvements are possible if the self-adjoint matrix Π(

*ω*) is found to have negative eigenvalue for any frequency range in the domain:

*ω*> 0. The Π test will be used to verify the existence of a periodic solution which outperforms the steady-state motion for the UAV problem.

## Problem Formulation

In this work, two cost functions are considered which reflect how far a UAV can fly (range) and how long it can stay up (endurance) in the air. This section presents a derivation of these costs, a list of the imposed constraints, and a comprehensive formulation of the optimization problem that needs to be solved to optimize endurance and range.

### Cost Function Formulation.

According to Ref. [48], endurance can be defined as the total time that an aircraft stays in the air on a tank of fuel. To improve performance of an aircraft, it is often desirable to maximize the endurance of the aircraft. To maximize the endurance, the average rate of fuel consumption per unit time should be minimized.

A major factor influencing fuel consumption is the thrust specific fuel consumption (TSFC) [48]. The TSFC is a characteristic of the power-plant. It is a measure of the fuel efficiency of an engine design with respect to thrust output and is mathematically represented as the fuel consumed per unit of thrust generated per unit time.

*σ*. The mass of fuel consumed over an incremental time

*dt*is then given by

*T*is the thrust. Integrating Eq. (23) with respect to time gives the mass of fuel used in a given time period

*T*

_{0}–

*T*. Therefore, the mass of fuel consumed (

_{f}*M*) is

_{f}*x-axis*direction is chosen to be that direction for illustration. Hence, the speed of the aircraft

*V*is substituted by the velocity in the

*x*-

*direction*making the final cost function (henceforth referred to as the range cost function) be

### Optimal Steady-State Flight.

For the optimal steady-state solutions, the cost functions given in Eqs. (25) and (28) are written in their steady-state form. This leads to the following optimization problems for endurance and range cost functions:

*Endurance*

*Range*

It should be noted that in all the optimization problems posed in this section, ($x\u02d9=0$) is not considered to be a constraint. This is because the steady-state flight objective is to fly without losing altitude in the *x*-direction. The other constraints are imposed to ensure flight in a constant configuration.

### Optimal Periodic Flight and Trajectory Parametrization Based on Flatness.

For the optimal periodic solutions, the cost functions given in Eqs. (25) and (28) are rewritten here to pose two optimization problems (*P _{E}*) and (

*P*), respectively, as

_{R}*Endurance (P*

_{E})*Range (P*

_{R})The constraint on the angle of attack (*α*) is imposed to maintain the validity of the UAV model which requires a small angle of attack.

where $\omega =(2\pi /Tf\u2212T0)$ The highest order derivative of the flat output is parametrized and lower order derivatives are obtained by successive integration. The periodicity constraint on the states “*y*” and “*z*” will allow only the constants of integration on the second integral to be nonzero, while all constants of integration for the *x* state can be nonzero. This is done to avoid polynomial terms in time in the states *y* and *z*. The coefficients are then determined by discretizing the time domain, imposing all the physical constraints at each distinct instant in time, imposing the periodicity constraints and finally solving the nonlinear optimization problem. The nonlinear optimization problem eventually results in solving for 6 *N* + 5 variables (6 *N* Fourier coefficients, four constants of integration, and the period *T _{f}* with

*T*

_{0}= 0).

## Results

The effectiveness of the proposed method is shown using simulation. The parameters used in the simulation are the Aerosonde UAV parameters from Beard and McLain [41]. These parameters are given in Table 1.

Parameters | Value |
---|---|

Air density (ρ) | 1.2682 kg/m^{3} |

Planform area (S) | 0.55 m^{2} |

Coefficient of parasitic drag (C_{d0}) | 0.03 |

Coefficient of lift at zero α (C_{l0}) | 0.28 |

Coefficient of α induced lift (C_{lα}) | 3.45 |

Mass (m) | 13.5 kg |

Gravitational acceleration (g) | 9.81 m/s^{2} |

Oswald efficiency factor (e) | 0.9 |

AR | 15.2445 |

Parameters | Value |
---|---|

Air density (ρ) | 1.2682 kg/m^{3} |

Planform area (S) | 0.55 m^{2} |

Coefficient of parasitic drag (C_{d0}) | 0.03 |

Coefficient of lift at zero α (C_{l0}) | 0.28 |

Coefficient of α induced lift (C_{lα}) | 3.45 |

Mass (m) | 13.5 kg |

Gravitational acceleration (g) | 9.81 m/s^{2} |

Oswald efficiency factor (e) | 0.9 |

AR | 15.2445 |

### Optimal Steady-State Flight.

For the steady profile, the cost, the values of the states, and input for the optimal solution obtained from the optimizer are given in Table 2.

### Π Test.

To construct the Π test, the system equations given by Eqs. (1)–(6) are linearized about the steady-state solutions. The expressions for the partial derivatives of the Hamiltonian and the linearized equations are then obtained. The equations for $x\u02d9$ and $y\u02d9$ are not considered as considering them to be steady implies that the aircraft is motionless. Since $V\u02d9$ is considered steady, it implies that $x\u02d9$ and $y\u02d9$ are constant but not both zero.

where ** f** is the system dynamics vector (i.e., a vector comprising the right-hand side of Eqs. (3)–(6)), $J\u0302$ is the integrand of any of the cost functions given in Eqs. (25) and (28), and $\lambda $ is the costate of the steady-state solution. To test local properness using the Π test, the self-adjoint matrix Π(

*ω*) is derived using Eq. (22).

According to Wolfe and Speyer [49], if the self-adjoint matrix Π(*ω*) is partially negative for some *ω* > 0, then the optimal periodic control problem is locally proper. This can be checked by observing the minimum eigenvalue of the self-adjoint matrix. If a *ω* > 0 can be found such that the minimum eigenvalue is negative, it can be asserted that the system is locally proper.

After deriving the Π(*ω*) matrix and observing the corresponding eigen values of the matrix, it can be seen that the minimum eigenvalue falls below zero for a finite range of *ω* (presented in Fig. 2(a)). Hence, it is inferred that a periodic solution will produce a better performance than a static optimum.

Once again, the minimum eigenvalue of the corresponding Π(*ω*) matrix is observed (presented in Fig. 2(b)) to fall below 0 for values of *ω* > 0: motivating the need to look for a periodic solution with a better cost value.

### Optimal Periodic Flight.

This subsection presents the results of the optimization problems *P _{E}* and

*P*. The effect of increasing the complexity of the parametrization on the optimal control design is studied and is presented in Sec. 6.3.1. Analysis on the final solution profiles obtained is investigated subsequently in Secs. 6.3.2 and 6.3.3.

_{R}#### Convergence Analysis.

*P _{E}* and

*P*are solved with varying number of Fourier terms in the parametrization to examine the effect of

_{R}*N*on the cost. From the definitions in Eqs. (33)–(35), the number of terms in the Fourier function is given by 2

*N*. The value of

*N*is varied between 1 and 20 and the effect on both cost functions is examined.

The plot for the evolution of the endurance cost as a function of *N* is shown in Fig. 3. All values of *N* yielded a better cost than the steady-state. Although it might appear that for *N *=* *1, the steady-state cost and the periodic cost are the same: they are in fact marginally different. *J _{N}*

_{= 1}= 0.08632 which is a 0.29% improvement over the steady-state cost. Figure 3 also shows that the cost value saturates as the number of Fourier terms increases. The magnitude of differential improvement in the cost keeps dropping with increasing values of

*N*. This can be attributed to the fact that, at low values of

*N*, there are not enough degrees-of-freedom in the parametrization while solving

*P*to capture the global optimal solution. However, this fact is remedied as more and more Fourier terms are included. As a result, we see the cost converging.

_{E}Figures 4(a)–4(c) present the control profiles (in normalized time) for a few selected *N* values. Considering all trajectories are periodic and are indifferent to shifts in time, the peaks of the thrust for all the *N* plots have been aligned for ease of comparison. It should be noted that all the other plots have been shifted appropriately corresponding to their shifts in thrust. The blue dotted lines represent the corresponding steady-state values and have been included for reference. It is interesting to note that in Fig. 4(a), as the number of Fourier terms is increased, the thrust profiles approach a scenario where the thrust is active for a single brief interval of time and is then inactive for the rest of the maneuver. In Fig. 4(b), a similar pattern is observed for the angle of attack (*α*). *α* remains at its maximum allowable limit throughout the maneuver when the thrust is inactive while it sharply drops in magnitude before restoring itself at its previous maximum during the thrusting phase.

Figures 5(a)–5(c) present the potential energy, kinetic energy, and the total energy profiles for the maneuver, respectively. It should be noted that the mean of the altitude is the zero potential energy reference, which is why the potential and total energy charts display negative values. In all of these plots, it can be seen that as the number of Fourier terms increases, the profile tends toward a particular shape with less oscillations. Moreover, on careful observation, one can learn that the total energy keeps dropping at all times except when the thrust is active where it rapidly increases. Intuitively, this makes sense: the UAV continuously loses energy while it overcomes drag and this lost energy is then pumped back into the system when the thrust is active.

Similar figures and patterns can also be observed when optimizing for range. Figure 6 shows the convergence of the range cost with increasing values of *N*.

Figures 7(a)–7(c) show the plots of the control profiles obtained after solving *P _{R}*. Once again, the plots have been shifted and time normalized for ease of visual comparison. Figures 8(a)–8(c) present the energy profiles for the maneuvers. Similar to endurance, it can be seen in these plots that as the number of Fourier terms is increased the profiles tend to converge to a particular shape: the thrust approaches a distinct active and an inactive region, and the energy correspondingly is pumped up and gradually drained in those regions, respectively. The angle of attack (

*α*) (Fig. 7(b)), however, has a more involved trajectory compared to the

*P*case. A distinctive double trough is observed during the thrusting phase followed by another gentler valley in the inactive phase.

_{E}#### Endurance.

To analyze the optimal solution for endurance, the result with 40 Fourier terms (*N *=* *20) was chosen. Figures 9(a)–9(c) present plots of the control input profiles, while Figs. 10(a)–10(c) present a subset of the state trajectories corresponding to the periodic solution. Figure 11 shows the trajectory of the optimal periodic flight: the starting point of the maneuver is demarcated by the circle on the left, while the end point is marked by the circle on the right.

The thrust profile illustrated in Fig. 9(a) demonstrates a distinct active propulsion phase amidst the coasting of the UAV. This active phase is shaded in gray in Figs. 9 and 10. The spatial trajectory is also plotted differently to distinguish the active and the inactive phase in Fig. 11. The dotted line represents the coasting phase, while the solid line represents the active phase. During the active phase, the UAV accumulates potential energy as shown by the increase in the altitude of the UAV in Fig. 11. During the coasting portion of the periodic control, Fig. 10(a) shows that the velocity of the UAV is essentially constant. The only increase in the velocity is during the active propulsion phase which also corresponds to an increase in kinetic energy. It can also be seen from Fig. 9(b) that the increase in velocity is accompanied by a reduction in the angle of attack. Examining the equation of both lift and drag, it is evident that the reduction in drag corresponding to lowering the angle of attack is more significant compared to the loss of lift due to the square term of *α* in the drag equation.

Figure 12 shows the plot of kinetic energy versus the potential energy. The dashed line represents the thrust inactive region, while the solid line represents the thrust active region. A color gradient also indicates the value of the thrust. Lines of constant total energy are plotted in the background for reference. From the plot, it can be seen that for the part of the trajectory where the thrust is not active the kinetic energy remains constant while potential energy changes. The strategy is, thus, to use the thrust to add kinetic energy to the system as it can be seen that the kinetic energy changes only when the thrust is active.

The optimal cost resulting from this trajectory is 0.05583 kg/s with a period of 309.2170 s: corresponding to a percentage improvement of 35.51% over the steady-state solution. This shows that with periodic control, a better performance can be obtained for the endurance of an unmanned aerial vehicle.

#### Range.

The optimal solution using 40 Fourier terms (*N *=* *20) was chosen. Figures 13(a)–13(c) show the control profile for the periodic solution, while Figs. 14(a)–14(c) show the plot of the state evolution of the periodic solution. The shaded region indicates the time period where thrust was active. Figure 15 shows the trajectory of the optimal periodic flight, and once again, the start and end points are demarcated by the circles on the right and left, respectively. The dotted part of the trajectory indicates the coasting phase, while the solid part represents the active phase. From the thrust profile illustrated in Fig. 13(a), a distinct active propulsion phase where the thrust is active can be seen, while the UAV is coasting at every other time. Examining Figs. 13(b) and 14(a), it can be seen that the increase in velocity is accompanied by a reduction in the angle of attack. This is a consequence of the lift been linear in *α*, while drag is quadratic in *α*. The reduction in drag corresponding to lowering the angle of attack is more significant compared to the loss of lift. In a portion of the nonactive phase, the UAV gained velocity by rapidly losing altitude hence converting potential energy into kinetic energy. During the active phase, the UAV gained kinetic energy and this kinetic energy is converted back to potential energy by increasing the altitude as shown in Fig. 15. Therefore, for range, there is a cyclic conversion between kinetic energy and potential energy where a rapid drop in altitude is used to gain velocity; hence, kinetic energy and the velocity gain are used to increase altitude thereby converting the kinetic energy to potential energy. The energy losses due to drag are compensated for by thrust in the active phase.

Figure 16 shows the plot of kinetic energy versus the potential energy. The dashed line shows the region where the thrust is inactive, while the color gradient indicates the value of the thrust. Lines of constant total energy are also indicated in the plot. From the plot, it can be seen that during the phase where the thrust is inactive part of the potential energy is converted into kinetic energy with some energy lost due to the effect of drag. Hence in this, both height reduction and thrust are used as mechanisms to increase velocity.

The optimal cost corresponding to the optimal range trajectory is 0.001935 kg/m with a period of 46.0305 s. This shows a percentage improvement of 51.95% over the steady-state solution. This shows that with the periodic control, a better performance can be obtained for the range of an unmanned aerial vehicle.

## Conclusion

In this work, three-dimensional optimal periodic solutions to enhance the endurance and range of the UAV were determined. The optimal trajectory generation was accomplished by exploiting the differential flatness of the system. The effect of the number of terms in the Fourier function was also considered and it was discovered that as the number of terms in the Fourier function increased, the cost converged.

The Π test was also carried out to examine whether periodic inputs can improve the performance of the UAV. The result showed that for both the considered cost functions, improvements are possible by giving periodic inputs which was shown in the simulation section to be true.

From the cost functions examined, it can be seen that optimal periodic endurance solution had an improvement of 35.51% with *N *=* *20, while range had an improvement of 51.95% with *N *=* *20. This improvement can be attributed to the steady-state solution involving constant usage of thrust as opposed to the periods of peaks and rests in the periodic solution. It should be noted that other locally optimal periodic solutions may exist and may be reached depending on the initial guess for coefficient of the Fourier functions as the optimization is very sensitive to initial conditions.

## Funding Data

Division of Civil, Mechanical and Manufacturing Innovation, National Science Foundation (NSF) (Grant No. CMMI-1537210; Funder ID: 10.13039/100000147).