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 [35] 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. [1721]. 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 [2426]. 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 [2931]. 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.

The differential equations governing the dynamics of the UAV (can be found in Ref. [40]) and are given by
x˙=Vcosγcosχ
(1)
y˙=Vcosγsinχ
(2)
z˙=Vsinγ
(3)
γ˙=gV(ncosϕcosγ)
(4)
χ˙=gVnsinϕcosγ
(5)
V˙=TDmgsinγ
(6)
where 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
n=LmgandL=12ρS(CL0+CLαα)V2
(7)

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 CL0 and C are the aerodynamic coefficients of lift. The expression for lift can be found in Ref. [41].

Fig. 1
Description of the aircraft states, inputs, coordinates, and forces
Fig. 1
Description of the aircraft states, inputs, coordinates, and forces
Close modal
The expression for drag (D) is also taken from the literature [41] and is given by
D=12ρS(CD0+(CL0+CLαα)2πeAR)V2
(8)
where e is the Oswald efficiency factor, AR is the wing aspect ratio, and CD0 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
U=[ϕ,α,T]andY=[x,y,z]
(9)
respectively.

Differential Flatness

A system is said to be differentially flat if its states and inputs can be expressed as functions of its outputs (also known as 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
x˙=f(x,u)withy=h(x)
(10)
where xRn is the state vector, uRm is the input vector, and y 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 yRm). Under this condition, the system given by Eq. (10) is differentially flat if a flat-output vector zRm can be identified of the form
z=η(x,u,u˙,u¨,u(p))
(11)
where
{x=x(z,z˙,z¨,.z(l))u=u(z,z˙,z¨,.z(l))
(12)

and (.)(a) refers to the ath 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 y. 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.

The UAV model (Eqs. (1)(6)) is differentially flat [31], with x, y, and z serving as the flat outputs, permitting the states and inputs to be written as a function of the flat-outputs (Y) and their derivatives as follows:
V=x˙2+y˙2+z˙2
(13)
γ=arcsin(z˙V)
(14)
χ=arcsin(y˙Vcosγ)
(15)
ϕ=arctan(χ˙Vcosγgcosγ+Vγ˙)
(16)
α=(2mρSV2cosϕCLα)(gcosγ+Vγ˙)CL0CLα
(17)
T=D+mV˙+mgsinγ
(18)

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

Before we look for better periodic solutions to an optimal control problem, there is a need to investigate if an optimal periodic solution exists which is better than the optimal steady-state solution. This investigation is done by conducting a test popularly known as the Π test. When a periodically varying input provides a better cost than the steady-state solution, the optimal periodic control problem is said to be proper, and if small periodic perturbations of the input around the steady-state solution leads to improvements in performance, then it is said to be locally proper. The Π test is a frequency domain criterion that determines local properness (first introduced in Ref. [43] and revisited in Ref. [44]). The test is valuable in situations where the first-order conditions for optimality prove deficient in giving information. More work on the Π test can also be found in the articles [4547]. To construct the Π test, the first step is to determine the steady-state solution of the optimal control problem of interest. Then, the system equations are linearized about the steady-state solution and the expressions for the partial derivatives of the Hamiltonian are obtained. For example, given a dynamic system (Eq. (10)), a cost function to be minimized
minJ=1τ0τg(y)dtτ>0
(19)
over the inputs, an optimal period τ and an output function
y=ν(x,u)
(20)
the Hamiltonian is defined as
H(x,u,y,λ,μ)=g(y)+λf(x,u)+μ(ν(x,u)y)
where λRn and μRm are Lagrange multipliers. The transfer function of the system linearized at the optimal steady-state solution is given by the equation
G(s)=(sIfx¯)1fu¯
(21)
where fx¯ and fu¯ are the partial derivatives of the system equation with respect to the states and the inputs, respectively. The bar denotes that the partial derivative is evaluated at the steady-state solution. For the Π test, a self-adjoint matrix Π(ω) is obtained using the expression
Π(ω)=G(jω)H¯xxG(jω)+H¯uxG(jω)+G(jω)H¯xu+H¯uu
(22)

where G is the transpose of G, and H¯xx,H¯ux,H¯xu,andH¯uu 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.

Let the TSFC be represented by σ. The mass of fuel consumed over an incremental time dt is then given by
dm=σTdt
(23)
where T is the thrust. Integrating Eq. (23) with respect to time gives the mass of fuel used in a given time period T0Tf. Therefore, the mass of fuel consumed (Mf) is
Mf=T0TfσTdt
(24)
This means that the average fuel consumption per unit time (i.e., the final cost function henceforth referred to as the endurance cost function) is given by
J=Mftime=1TfT0T0TfσTdt
(25)
The second objective of interest is the range of the aircraft. Range can be defined as the total distance moved by an aircraft on a tank of fuel. To maximize the range, the time average rate of fuel consumption per unit distance should be minimized. To derive the range equation, we note that
V=dRdtdt=dRV
where R is the range. Therefore, Eq. (23) can be rewritten as
dm=σTVdR
(26)
Then, the time average fuel consumption per unit range is given as
1TfT0dmdRdt=1TfT0T0TfσTVdt
(27)
Equation (27) only captures the average fuel consumption per unit distance along the direction the aircraft travels. However, a more appropriate cost function would be to determine the average fuel consumption per unit displacement along a particular direction of interest (could be the direction in which the aircraft was supposed to travel). In this work, the 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
J=1TfT0T0Tf[σTVcos(γ)cos(χ)]dt
(28)

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
minα,ϕ,T,z,γ,χ,V[σT]
(29)
Range
minα,ϕ,T,z,γ,χ,V[σTVcos(γ)cos(χ)]
(30)
subject to
z˙=0,γ˙=0,χ˙=0,V˙=0

It should be noted that in all the optimization problems posed in this section, (x˙=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 (PE) and (PR), respectively, as

Endurance (PE)
minα,ϕ,T,(TfT0)1TfT0T0Tf[σT]dt
(31)
Range (PR)
minα,ϕ,T,(TfT0)1TfT0T0Tf[σTVcos(γ)cos(χ)]dt
(32)
subject to
Periodic constraints
y(T0)=y(Tf)
z(T0)=z(Tf)
V(T0)=V(Tf)
γ(T0)=γ(Tf)
χ(T0)=χ(Tf)
Physical constraints on altitude, velocity, thrust, angle of attack, and displacement are also imposed
2000z0
V0
140T0
π18απ18
x(t+δt)x(t)

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

This optimization problem is solved by making use of the differential flatness of the UAV model and parametrizing the flat outputs by a linear combination of basis functions in time casting it into a nonlinear optimization problem. Since, we seek a periodic solution to the problem, the natural choice of basis functions are elements of the Fourier series. Typically, the highest necessary derivative of the flat outputs is expressed as an unbiased Fourier expansion (see Ref. [12]). Therefore, in this case, the flat outputs are written as
x¨=i=1N[a2i1sin(iωt)+a2icos(iωt)]
(33)
y¨=i=1N[b2i1sin(iωt)+b2icos(iωt)]
(34)
z¨=i=1N[c2i1sin(iωt)+c2icos(iωt)]
(35)

where ω=(2π/TfT0) 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 Tf with T0 = 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.

Table 1

Values of parameters used in simulation

ParametersValue
Air density (ρ)1.2682 kg/m3
Planform area (S)0.55 m2
Coefficient of parasitic drag (Cd0)0.03
Coefficient of lift at zero α (Cl0)0.28
Coefficient of α induced lift (C)3.45
Mass (m)13.5 kg
Gravitational acceleration (g)9.81 m/s2
Oswald efficiency factor (e)0.9
AR15.2445
ParametersValue
Air density (ρ)1.2682 kg/m3
Planform area (S)0.55 m2
Coefficient of parasitic drag (Cd0)0.03
Coefficient of lift at zero α (Cl0)0.28
Coefficient of α induced lift (C)3.45
Mass (m)13.5 kg
Gravitational acceleration (g)9.81 m/s2
Oswald efficiency factor (e)0.9
AR15.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.

Table 2

Values of optimal steady-state solution

StatesInputs
FunctionCostzγχVTαϕ
Endurance0.0865700020.74857.21440.17450
Range0.00402600024.04998.06880.10910
StatesInputs
FunctionCostzγχVTαϕ
Endurance0.0865700020.74857.21440.17450
Range0.00402600024.04998.06880.10910

Π 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˙ and y˙ are not considered as considering them to be steady implies that the aircraft is motionless. Since V˙ is considered steady, it implies that x˙ and y˙ are constant but not both zero.

The Hamiltonian for the system is given by
H=Ĵ+λf

where f is the system dynamics vector (i.e., a vector comprising the right-hand side of Eqs. (3)(6)), Ĵ is the integrand of any of the cost functions given in Eqs. (25) and (28), and λ 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.

The results from the Π test for endurance and range are illustrated subsequently. For the endurance cost function, the Hamiltonian is given by
H=σT+λf
(36)

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.

Fig. 2
(a) Minimum eigenvalue of Π(jω) for endurance and (b) minimum eigenvalue of Π(jω) for range
Fig. 2
(a) Minimum eigenvalue of Π(jω) for endurance and (b) minimum eigenvalue of Π(jω) for range
Close modal
For the range cost function, the Hamiltonian is given by
H=σTV+λf
(37)

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 PE and PR. 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.

Convergence Analysis.

PE and PR are solved with varying number of Fourier terms in the parametrization to examine the effect of 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. JN= 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 PE 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.

Fig. 3
Effect of the N-count on the endurance cost
Fig. 3
Effect of the N-count on the endurance cost
Close modal

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.

Fig. 4
Control solution profiles for PE (varying N): (a) thrust, (b) angle of attack (α), and (c) bank angle (ϕ)
Fig. 4
Control solution profiles for PE (varying N): (a) thrust, (b) angle of attack (α), and (c) bank angle (ϕ)
Close modal

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.

Fig. 5
Energy profiles for PE (varying N): (a) potential energy, (b) kinetic energy, and (c) total energy
Fig. 5
Energy profiles for PE (varying N): (a) potential energy, (b) kinetic energy, and (c) total energy
Close modal

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.

Fig. 6
Effect of the N-count on the range cost
Fig. 6
Effect of the N-count on the range cost
Close modal

Figures 7(a)7(c) show the plots of the control profiles obtained after solving PR. 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 PE case. A distinctive double trough is observed during the thrusting phase followed by another gentler valley in the inactive phase.

Fig. 7
Control solution profiles for PR (varying N): (a) thrust, (b) angle of attack (α), and (c) bank angle (ϕ)
Fig. 7
Control solution profiles for PR (varying N): (a) thrust, (b) angle of attack (α), and (c) bank angle (ϕ)
Close modal
Fig. 8
Energy profiles for PR (varying N): (a) potential energy, (b) kinetic energy, and (c) total energy
Fig. 8
Energy profiles for PR (varying N): (a) potential energy, (b) kinetic energy, and (c) total energy
Close modal

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.

Fig. 9
Control profiles after solving PE for N = 20: (a) thrust, (b) angle of attack (α), and (c) bank angle (ϕ)
Fig. 9
Control profiles after solving PE for N = 20: (a) thrust, (b) angle of attack (α), and (c) bank angle (ϕ)
Close modal
Fig. 10
State trajectories after solving PE for N = 20: (a) velocity, (b) heading angle (χ), and (c) flight path angle (γ)
Fig. 10
State trajectories after solving PE for N = 20: (a) velocity, (b) heading angle (χ), and (c) flight path angle (γ)
Close modal
Fig. 11
Optimal trajectory for endurance with 40 Fourier terms
Fig. 11
Optimal trajectory for endurance with 40 Fourier terms
Close modal

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.

Fig. 12
Plot of endurance kinetic energy versus potential energy
Fig. 12
Plot of endurance kinetic energy versus potential energy
Close modal

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.

Fig. 13
Control profiles after solving PR for N = 20: (a) thrust, (b) angle of attack (α), and (c) bank angle (ϕ)
Fig. 13
Control profiles after solving PR for N = 20: (a) thrust, (b) angle of attack (α), and (c) bank angle (ϕ)
Close modal
Fig. 14
Plot showing the range: (a) velocity (b) heading angle, and (c) flight path angle for 40 Fourier terms
Fig. 14
Plot showing the range: (a) velocity (b) heading angle, and (c) flight path angle for 40 Fourier terms
Close modal
Fig. 15
Optimal trajectory for range with 40 Fourier terms
Fig. 15
Optimal trajectory for range with 40 Fourier terms
Close modal

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.

Fig. 16
Plot of range kinetic energy versus potential energy
Fig. 16
Plot of range kinetic energy versus potential energy
Close modal

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).

References

1.
Horcher
,
A.
, and
Visser
,
R. J.
,
2004
, “
Unmanned Aerial Vehicles: Applications for Natural Resource Management and Monitoring
,”
Council on Forest Engineering Proceedings
, Hot Springs, AR, Apr. 27–30.
2.
Adams
,
S. M.
, and
Friedland
,
C. J.
,
2011
, “
A Survey of Unmanned Aerial Vehicle (UAV) Usage for Imagery Collection in Disaster Research and Management
,”
Ninth International Workshop on Remote Sensing for Disaster Response
, Stanford, CA, Sept. 15–16, p.
8
.https://www.researchgate.net/publication/266465037_A_Survey_of_Unmanned_Aerial_Vehicle_UAV_Usage_for_Imagery_Collection_in_Disaster_Research_and_Management
3.
Speyer
,
J. L.
,
1976
, “
Nonoptimality of the Steady-State Cruise for Aircraft
,”
AIAA J.
,
14
(
11
), pp.
1604
1610
.
4.
Speyer
,
J. L.
,
1996
, “
Periodic Optimal Flight
,”
J. Guid. Control Dyn.
,
19
(
4
), pp.
745
755
.
5.
Menon
,
P. K.
,
Sweriduk
,
G. D.
, and
Bowers
,
A. H.
,
2007
, “
Study of Near-Optimal Endurance-Maximizing Periodic Cruise Trajectories
,”
J. Aircr.
,
44
(
2
), pp.
393
398
.
6.
Rayner
,
J. M.
,
Viscardi
,
P. W.
,
Ward
,
S.
, and
Speakman
,
J. R.
,
2001
, “
Aerodynamics and Energetics of Intermittent Flight in Birds
,”
Am. Zool.
,
41
(
2
), pp.
188
204
.
7.
Gleiss
,
A. C.
,
Jorgensen
,
S. J.
,
Liebsch
,
N.
,
Sala
,
J. E.
,
Norman
,
B.
,
Hays
,
G. C.
,
Quintana
,
F.
,
Grundy
,
E.
,
Campagna
,
C.
,
Trites
,
A. W.
, and
Block
,
B. A.
,
2011
, “
Convergent Evolution in Locomotory Patterns of Flying and Swimming Animals
,”
Nat. Commun.
,
2
, p.
352
.
8.
Traugott
,
J.
,
Nesterova
,
A.
, and
Sachs
,
G.
,
2013
, “
The Flight of the Albatross
,”
IEEE Spectrum
,
50
(
7
), pp.
46
54
.
9.
Paoletti
,
P.
, and
Mahadevan
,
L.
,
2014
, “
Intermittent Locomotion as an Optimal Control Strategy
,”
Proc. R. Soc. London, A
,
470
(
2164
), p.
20130535
.
10.
Wu
,
T. Y.
,
2011
, “
Fish Swimming and Bird/Insect Flight
,”
Annu. Rev. Fluid Mech.
,
43
(
1
), pp.
25
58
.
11.
Williams
,
T. M.
,
2001
, “
Intermittent Swimming by Mammals: A Strategy for Increasing Energetic Efficiency During Diving
,”
Am. Zool.
,
41
(
2
), pp.
166
176
.
12.
Varigonda
,
S.
,
Georgiou
,
T. T.
,
Siegel
,
R. A.
, and
Daoutidis
,
P.
,
2008
, “
Optimal Periodic Control of a Drug Delivery System
,”
Comput. Chem. Eng.
,
32
(
10
), pp.
2256
2262
.
13.
Li
,
S. E.
,
Xu
,
S.
,
Li
,
G.
, and
Cheng
,
B.
,
2014
, “
Periodicity Based Cruising Control of Passenger Cars for Optimized Fuel Consumption
,”
IEEE
Intelligent Vehicles Symposium
, Dearborn, MI, June 8–11, pp.
1097
1102
.
14.
Bayen
,
T.
,
Mairet
,
F.
,
Martinon
,
P.
, and
Sebbah
,
M.
,
2015
, “
Analysis of a Periodic Optimal Control Problem Connected to Microalgae Anaerobic Digestion
,”
Optim. Control Appl. Methods
,
36
(
6
), pp.
750
773
.
15.
He
,
F.
,
Leung
,
A.
, and
Stojanovic
,
S.
,
1994
, “
Periodic Optimal Control for Competing Parabolic Volterra-Lotka-Type Systems
,”
J. Comput. Appl. Math.
,
52
(
1–3
), pp.
199
217
.
16.
Xiao
,
Y.
,
Cheng
,
D.
, and
Qin
,
H.
,
2006
, “
Optimal Impulsive Control in Periodic Ecosystem
,”
Syst. Control Lett.
,
55
(
7
), pp.
558
565
.
17.
Burrows
,
J. W.
,
1983
, “
Fuel-Optimal Aircraft Trajectories With Fixed Arrival Times
,”
J. Guid. Control Dyn.
,
6
(
1
), pp.
14
19
.
18.
Chakravarty
,
A.
,
1985
, “
Four-Dimensional Fuel-Optimal Guidance in the Presence of Winds
,”
J. Guid. Control Dyn.
,
8
(
1
), pp.
16
22
.
19.
Dewell
,
L.
, and
Speyer
,
J.
,
1993
, “
An Investigation of the Fuel-Optimal Periodic Trajectories of a Hypersonic Vehicle
,”
Guidance, Navigation and Control Conference
, Monterey, CA, Aug. 9–11, p.
3753
.
20.
Harada
,
M.
, and
Bollino
,
K.
,
2009
, “
Fuel Optimization of Figure-8 Flight for Unmanned Aerial Vehicles
,”
AIAA
Paper No. 2009–6011.
21.
Sachs
,
G.
, and
Lesch
,
K.
,
1990
, “
Fuel Savings by Optimal Aircraft Cruise With Singular and Chattering Control
,”
Anal. Optim. Syst.
,
144
, pp.
590
599
.
22.
Chen
,
R. H.
, and
Speyer
,
J. L.
,
2007
, “
Improved Endurance of Optimal Periodic Flight
,”
J. Guid. Control Dyn.
,
30
(
4
), pp.
1123
1133
.
23.
Fliess
,
M.
,
Lévine
,
J.
,
Martin
,
P.
, and
Rouchon
,
P.
,
1995
, “
Flatness and Defect of Non-Linear Systems: Introductory Theory and Examples
,”
Int. J. Control
,
61
(
6
), pp.
1327
1361
.
24.
Chamseddine
,
A.
,
Zhang
,
Y.
,
Rabbath
,
C. A.
,
Join
,
C.
, and
Theilliol
,
D.
,
2012
, “
Flatness-Based Trajectory Planning/Replanning for a Quadrotor Unmanned Aerial Vehicle
,”
IEEE Trans. Aerosp. Electron. Syst.
,
48
(
4
), pp.
2832
2848
.
25.
Cowling
,
I. D.
,
Yakimenko
,
O. A.
,
Whidborne
,
J. F.
, and
Cooke
,
A. K.
,
2007
, “
A Prototype of an Autonomous Controller for a Quadrotor UAV
,”
European Control Conference
(
ECC
), Kos, Greece, July 2–5, pp.
4001
4008
.
26.
Kingston
,
D.
,
Beard
,
R.
,
McLain
,
T.
,
Larsen
,
M.
, and
Ren
,
W.
,
2003
, “
Autonomous Vehicle Technologies for Small Fixed-Wing UAVs
,”
AIAA
Paper No. 2003–6559.
27.
Whidborne
,
I. D. C. J. F.
, and
Cooke
,
A. K.
,
2006
, “
Optimal Trajectory Planning and Lqr Control for a Quadrotor UAV
,”
International Conference Control
, Glasgow, Scotland, Aug. 30–Sept. 1, Paper No. 125.https://www.researchgate.net/publication/238677911_OPTIMAL_TRAJECTORY_PLANNING_AND_LQR_CONTROL_FOR_A_QUADROTOR_UAV
28.
Hao
,
Y.
,
Davari
,
A.
, and
Manesh
,
A.
,
2005
, “
Differential Flatness-Based Trajectory Planning for Multiple Unmanned Aerial Vehicles Using Mixed-Integer Linear Programming
,” IEEE
American Control Conference
(
ACC
), Portland, OR, June 8–10, pp.
104
109
.
29.
Singh
,
L.
, and
Fuller
,
J.
,
2001
, “
Trajectory Generation for a UAV in Urban Terrain, Using Nonlinear MPC
,” IEEE
American Control Conference
(
ACC
), Vol.
3
, Arlington, VA, June 25–27, pp.
2301
2308
.
30.
Zhang
,
C.
,
Wang
,
N.
, and
Chen
,
J.
,
2010
, “
Trajectory Generation for Aircraft Based on Differential Flatness and Spline Theory
,”
International Conference on Information, Networking and Automation
(
ICINA
), Kunming, China, Oct. 18–19, pp.
110
114
.
31.
Alturbeh
,
H.
, and
Whidborne
,
J. F.
,
2014
, “
Real-Time Obstacle Collision Avoidance for Fixed Wing Aircraft Using b-Splines
,”
UKACC International Conference in Control
(
CONTROL
), Loughborough, UK, July 9–11, pp.
115
120
.
32.
Sira-Ramirez
,
H.
, and
Fliess
,
M.
,
1999
, “
Regulation of Non-Minimum Phase Outputs in a PVTOL Aircraft
,”
37th IEEE Conference on Decision and Control
(
CDC
), Tampa, FL, Dec. 18, pp.
4222
4227
.
33.
Deligiannis
,
V.
,
Davrazos
,
G.
,
Manesis
,
S.
, and
Arampatzis
,
T.
,
2006
, “
Flatness Conservation in the n-Trailer System Equipped With a Sliding Kingpin Mechanism
,”
J. Intell. Rob. Syst.
,
46
(
2
), pp.
151
162
.
34.
Ryu
,
J.-C.
,
Agrawal
,
S. K.
, and
Franch
,
J.
,
2008
, “
Motion Planning and Control of a Tractor With a Steerable Trailer Using Differential Flatness
,”
ASME J. Comput. Nonlinear Dyn.
,
3
(
3
), p.
031003
.
35.
Ryu
,
J.-C.
, and
Agrawal
,
S. K.
,
2011
, “
Differential Flatness-Based Robust Control of Mobile Robots in the Presence of Slip
,”
Int. J. Rob. Res.
,
30
(
4
), pp.
463
475
.
36.
Zimmert
,
N.
, and
Sawodny
,
O.
,
2010
, “
Active Damping Control for Bending Oscillations of a Forklift Mast Using Flatness Based Techniques
,”
American Control Conference
(
ACC
), Baltimore, MD, June 30–July 2, pp.
1538
1543
.
37.
Thounthong
,
P.
,
Pierfederici
,
S.
,
Martin
,
J.-P.
,
Hinaje
,
M.
, and
Davat
,
B.
,
2010
, “
Modeling and Control of Fuel Cell/Supercapacitor Hybrid Source Based on Differential Flatness Control
,”
IEEE Trans. Veh. Technol.
,
59
(
6
), pp.
2700
2710
.
38.
Murray
,
R. M.
,
Rathinam
,
M.
, and
Sluis
,
W.
,
1995
, “
Differential Flatness of Mechanical Control Systems: A Catalog of Prototype Systems
,”
ASME International Congress and Exposition
, San Francisco, CA, Nov. 12–17.
39.
Menon
,
P.
,
Sweriduk
,
G.
, and
Sridhar
,
B.
,
1999
, “
Optimal Strategies for Free-Flight Air Traffic Conflict Resolution
,”
J. Guid. Control Dyn.
,
22
(
2
), pp.
202
211
.
40.
Bicchi
,
A.
, and
Pallottino
,
L.
,
2000
, “
On Optimal Cooperative Conflict Resolution for Air Traffic Management Systems
,”
IEEE Trans. Intell. Transp. Syst.
,
1
(
4
), pp.
221
231
.
41.
Beard
,
R. W.
, and
McLain
,
T. W.
,
2012
,
Small Unmanned Aircraft: Theory and Practice
,
Princeton University Press
,
Oxfordshire, UK
.
42.
Rigatos
,
G. G.
,
2015
,
Differential Flatness Theory and Flatness-Based Control
,
Springer International Publishing
,
Cham, Switzerland
, pp.
47
101
.
43.
Bittanti
,
S.
,
Fronza
,
G.
, and
Guardabassi
,
G.
,
1973
, “
Periodic Control: A Frequency Domain Approach
,”
IEEE Trans. Autom. Control
,
18
(
1
), pp.
33
38
.
44.
Bernstein
,
D.
, and
Gilbert
,
E.
,
1980
, “
Optimal Periodic Control: The π Test Revisited
,”
IEEE Trans. Autom. Control
,
25
(
4
), pp.
673
684
.
45.
Bernstein
,
D.
,
1985
, “
Control Constraints, Abnormality, and Improved Performance by Periodic Control
,”
IEEE Trans. Autom. Control
,
30
(
4
), pp.
367
376
.
46.
Colonius
,
F.
,
1985
,
Optimal Periodic Control
,
Fachbereiche Mathematik/Informatik, Elektrotechnik/Physik. Forschungsschwerpunkt Dynamische Systeme, Universität Bremen
, Bremen, Germany.
47.
Speyer
,
J.
, and
Evans
,
R.
,
1984
, “
A Second Variational Theory for Optimal Periodic Processes
,”
IEEE Trans. Autom. Control
,
29
(
2
), pp.
138
148
.
48.
Anderson
,
J. D.
,
2012
,
Introduction to Flight
, 7th ed.,
McGraw-Hill
, New York.
49.
Wolfe
,
J. D.
, and
Speyer
,
J. L.
,
2003
, “
The Periodic Optimality of LQ Controllers Satisfying Strong Stabilization
,”
Automatica
,
39
(
10
), pp.
1735
1746
.