# Finite Horizon Nonlinear Suboptimal Control for an Autonomous Soaring UAV

Authors: David Ángeles-Rojas; Omar-Jacobo Santos-Sánchez; Sergio Salazar; Rogelio Lozano
Journal: Mathematical Problems in Engineering
 (2022)

---

## Abstract

In this paper, a suboptimal nonlinear discrete control for nonlinear discrete affine systems is implemented in an autonomous soaring unmanned aerial vehicle (UAV) to improve energy consumption and performance. General expressions for this controller are obtained from the continuous nonlinear guidance model of a fixed-wing UAV and applied in a multiloop control structure. Simulation results in the task of trajectory tracking inside a thermal updraft and a comparative study with a proportional derivative (PD) controller validated the effectiveness of this method.


---

## Body

## 1. Introduction

In recent years, researchers have become increasingly interested in unmanned aerial vehicles (UAVs) due to the number of applications they can be used for and the advantages in using them. However, with all their advantages come some disadvantages: one of the more representatives is the flight time, which is limited to the battery capacity installed in the UAV and the power consumption of its components. To improve the endurance of fixed-wing UAVs, different techniques have been proposed in research works; one of them is to detect and identify thermals to exploit their energy and this is known as static soaring. A method to identify thermals and estimation of its parameters is presented in [1]; once the algorithm has identified its parameters, a circular flight path for soaring is generated. In [2], an algorithm to detect thermals is proposed and it can estimate its parameters by flying twice inside the thermal; the simulation results showed the effectiveness of the algorithm. A method that uses the sink polar of a UAV to detect thermals along with the measuring of total energy was proposed in [3]; for this, they use the potential energy and the kinetic energy, and the simulations on Condor soaring simulator and MatLab/Simulink showed the efficacy of the proposed method. In [4], an autonomous soaring system is proposed; meteorological data are given to the UAV to create a dynamic map of lift surfaces where the UAV could find thermals. Using this along with their system described for identifying thermals and their parameters, the UAV can exploit the updrafts.Besides exploiting energy from thermal updrafts to save energy from batteries, there is another important method to extend the flight time, and it is the implementation of optimal controllers, which can help to achieve the maximum endurance of fixed-wing UAVs. There are some works that have studied this aspect; in [5], they addressed the problem to calculate an energy optimal trajectory for a UAV, and they tried to use energy in an optimal way by minimizing the energy used for propulsion. To obtain the optimal guidance, they use the kinematic model in 2D and consider coordinated turn conditions. A dynamic programming approach is used in [6] to find an optimal trajectory for a UAV considering wind and aircraft parameters, and the cost function proposed is evaluated based on the kinematic model and forecast data. An optimized controller is obtained in [7] using deep reinforcement learning, and it is implemented to the controlling attitude of a fixed-wing UAV; its performance is compared with a proportional integral derivative (PID) controller. In [8], an optimal controller is implemented to obtain the minimum time and fuel consumption during the climb phase of an aircraft where they take into consideration the thrust and the lift coefficient as control inputs, and the results are presented in simulations. An optimal control is designed in [9] based on the linearized models of a fixed-wing UAV system and implemented to follow trajectories.An adaptive vector field control law is proposed in [10]; this is an interesting work because the performance of this guidance control law does not depend on the knowledge of the dynamics model and the wind, simulations, and experimental results demonstrate that the performance of this control law proposed is better than the standard version and can compensate the lack of knowledge of the wind vector and unmodeled course angle dynamics. In [11], adaptive formation controllers for fixed-wing UAVs and a software-in-the-loop are designed and implemented, these adaptive formation controllers are designed to lead with parametric uncertainties such as mass and inertia, and software-in-the-loop simulations show the effectiveness of these controllers handling uncertain mass and inertia in the task of path following with multiple UAVs.In this work, the suboptimal control method proposed minimizes the performance index for nonlinear discrete affine systems using the dynamic programming approach without the need to solve the Bellman equation; for linear systems, the problem has been solved but the backward solution has the problem called “course of dimensionality” where it is necessary to solve a Riccati equation, and the more large is the final horizon (final time), the more the computational power is needed. For this reason, it is common to solve the Riccati equation offline but now it is necessary to have enough storage space to save the results for every step. The control equations obtained here for a fixed-wing UAV do not have those problems and can be implemented in autopilots without much computational power or storage space available like the Pixhawk autopilot; furthermore, this suboptimal control method uses nonlinear dynamics that implies better performance.This paper is divided into 4 sections. Section1 provides a brief introduction to the problem we focused on and previous published related papers. Section 2 is dedicated to the suboptimal nonlinear discrete controller and the control expressions for the implementation using the continuous nonlinear guidance model. Simulation results and a comparative study with a PD (proportional derivative) controller are presented in Section 3. Finally, conclusions are given in Section 4.


## 2. Materials and Methods

### 2.1. Optimization Procedure of Finite Horizon

Consider the nonlinear discrete affine system(1)x¯k+1=f0x¯k+f1x¯kuk,where f0x¯k∈ℝn, f1x¯k∈ℝn×m, and uk∈ℝm. Now, the performance index is defined as in [12](2)J=12x¯NHx¯N+12∑k=0N−1x¯TkQx¯k+uTkRuk,where H,Q≥0 and R>0 are symmetric matrices of appropriate dimensions, the horizon is tf=TsN with Ts defining the sample period, N is the final step, and k∈ℕ, where k=0,1,…,N. The suboptimal control sequence is obtained as in [13] using the dynamic programming approach in discrete time proposed in [12]; in the last step, the performance index is(3)JN,N∗=12x¯NHx¯N.This term does not depend on the control lawuN, then we can say JN,N∗ is the optimal value of J in the discrete time N. For the next step, the performance index is(4)JN−1,N∗=minuN−1JN∗+12x¯TN−1Qx¯N−1+12uTN−1RuN−1.An explicit form ofJN−1,N∗ can be obtained by using (1)(5)JN−1,N∗=minuN−112f0x¯N−1+f1x¯N−1uN−1T×Hf0x¯N−1+f1x¯N−1uN−1+12x¯TN−1Qx¯N−1+12uTN−1RuN−1.To find the optimal value ofJN−1,N, it is calculated in the traditional way with the partial derivative of (5) with respect to uN−1 equal to zero.(6)u∗N−1=−f1Tx¯N−1Hf1x¯N−1+R−1f1Tx¯N−1Hf0x¯N−1.Equation (6) is the control law that minimizes JN−1,N in the discrete time N−1. With the condition set before R>0, we can be sure the inverse(7)f1Tx¯N−1Hf1x¯N−1+R−1,exists. With the purpose of verifying that (6) is minimum, JN−1,N is derived partially with respect to uN−1 for a second time; then, we have(8)∂2JN−1,N∂u2N−1=f1Tx¯N−1Hf1x¯N−1+R>0.The control given by equation (6) is the optimal control in local sense. For the next step N−2, the performance index is(9)JN−2,N∗=minuN−212x¯NHx¯N+12x¯TN−1Qx¯N−1+12uTN−1RuN−1+12x¯TN−2Qx¯N−2+12uTN−2RuN−2.Using (1) in (9), we have(10)JN−2,N∗=minuN−212f0x¯N−1+f1x¯N−1uN−1THf0x¯N−1+f1x¯N−1uN−1+12f0x¯N−2+f1x¯N−2uN−2TQf0x¯N−2+f1x¯N−2uN−2+12f0x¯N−3+f1x¯N−3uN−3TQf0x¯N−3+f1x¯N−3uN−3+12uTN−1RuN−1+12uTN−2RuN−2.Again we proceed in the traditional way to find the minimal value ofJN−2,N. In this case, the second partial derivative is derived with respect to uN−2; then, we have that(11)uN−2=−f1Tx¯N−2Qf1x¯N−2+R−1f1Tx¯N−2Qf0x¯N−2.Repeating this task for a few more steps, we realized that the general expressions are as follows.Fork=1,(12)uN−k=−f1Tx¯N−kHf1x¯N−k+R−1f1Tx¯N−kHf0x¯N−k,and, for k=2,…,N,(13)uN−k=−f1Tx¯N−kQf1x¯N−k+R−1f1Tx¯N−kQf0x¯N−k.This suboptimal control sequence guarantees that an approximation to the minimal value of the performance index (2) is reached.


### 2.2. Mathematical Model of a Fixed-Wing UAV

The previously exposed control strategy is implemented in a three-layer control structure composed of an inner loop, a middle loop, and an outer loop.


#### 2.2.1. Continuous Nonlinear Model

For the inner loop, we used a part of the dynamic model presented in [14, 15], where due to the symmetry of the UAV Ixz is very much smaller than Ixx, Iyy, and Izz and can be neglected:(14)θ..=Izz−IxxIyyψvϕ˙+1Iyyτθ,ϕ¨=Iyy−IzzIxxψvθ.+1Ixxτϕ,where Ixx, Iyy, and Izz are the inertia parameters, ψv is the velocity in angle yaw, θ is the pitch angle, ϕ is the roll angle, and τθ and τϕ are control inputs for pitch and roll, respectively. The purpose of this inner loop is to follow the commands given by the middle loop and outer loop.The middle loop is dedicated to control the course angleχ. As it is described in [16], the course angle χ is the direction of travel relative to the axis x, the yaw or heading angle ψ is the direction that the UAV is pointed, and the side-slip angle β is the difference between the course and the heading angles. Figure 1 depicts this.Figure 1
Definition of the headingψ, side-slip β, and course χ angles.In this paper, we consider coordinated turn conditions, absence of wind, and side-slip. Then, we have thatψ=χ; with this, the kinematic model is given by(15)ψ˙=gVatanϕ,ϕ˙=1τsϕϕc−ϕ,where g is the magnitude of gravity at the sea level, τsϕ is the constant time of the system, ϕc is the control angle commanded, and Va is the airspeed which under the assumption of absence of wind this value is the UAV ground speed Va=Vg.Finally, for the outer loop, we considered the next kinematic model for guidance based on the one proposed in [16](16)x˙=Vacosψcosθ,y˙=Vasinψcosθ,z˙=Vasinθ,Va.=Tm−Dm−gsinθ,ψ˙=1τsψψc−ψ,θ.=1τsθθc−θ,where x, y, and z denote the inertial position of the aircraft in a 3D inertial frame, m is the mass of the UAV, D represents the drag force, ψ is the heading angle (or course angle under the assumptions of this paper), θ is the pitch angle, and τsθ and τsψ are constants time of the system. In this kinematic model, we have 3 control inputs: T is the thrust motor and θc and ψc are commanded angles. In this paper, we use the course angle for trajectory tracking, but with the suppositions made before we have that ψ=χ.In order to simplify the process to obtain the control expressions for the aircraft, the model is divided into 3 subsystems as follows:SubsystemVa:(17)Va.=Tm−Dm−gsinθ.This subsystem is dedicated to control the velocity of the aircraftVa and its control input is the thrust force T.Subsystemz:(18)z˙=Vasinθ,θ.=1τsθθc−θ.This subsystem is for control the altitude position, and the control input isθc.Subsystemxy:(19)x˙=Vacosψcosθ,y˙=Vasinψcosθ,ψ˙=1τsψψc−ψ.This subsystem is devoted to controlling the position in the planex−y and has as a control input ψc.In the trajectory tracking problem, there are basically two kinds of paths to follow, straight-line segments and circular orbits; with a combination of these, more complex paths can be tracked. Then, it is necessary to represent the model given in (19) in a form where we can know the cross-track error in an easy way.For straight-line segments, we used the matrix rotation given in (20) to transform the inertial frame to the path frame:(20)Rline=cosψtsinψt−sinψtcosψt,where ψt is the course angle of the path. Using (20), we calculated a new representation of (19), where the important dynamics are(21)y˙ψ=Vacosθsinψ−ψt,ψ˙=1τsψψc−ψ.For the circular orbit tracking problem, due to its nature, it is better to represent (19) in a polar coordinate system; for this, we use the matrix rotation given by(22)Rorbit=cosφsinφ−sinφcosφ,where φ is the phase angle of the aircraft. Using (22), a new representation of (19) is calculated, where the important dynamics are(23)d˙=Vacosθcosψ−φ,ψ˙=1τsψψc−ψ,where d is the distance from the UAV to the center of the coordinated system.


#### 2.2.2. Discrete Nonlinear Model

The nonlinear discrete affine representation of the continuous nonlinear models was obtained by using the Euler approximation. For the inner loop, the nonlinear discrete affine state-space representation of the dynamic model (14) is(24)θ1k+1θ2k+1=Tsθ2k+θ1kTsIzz−IxxIyyψvkϕ2k+θ2k+0TsIyyτθk,ϕ1k+1ϕ2k+1=Tsϕ2k+ϕ1kTsIyy−IzzIxxψvkθ2k+ϕ2k+0TsIxxτϕk,respectively, where Ts is the sample time. The nonlinear discrete affine representation of (15) is given by(25)ψk+1ϕk+1=TsgVaktanϕk+ψk1−Tsτsϕϕk+0Tsτsϕϕck.The subsystemVa given in (17) has the discrete representation as follows:(26)Vak+1=−TsmD−Tsgsinθk+Vak+TsmTk.For subsystem (18), its nonlinear discrete affine representation is(27)zk+1θk+1=TsVaksinθk+zk1−Tsτsθθk+0Tsτsθθck.The subsystemxy has two versions: one for straight-line trajectory tracking and the other one for circular orbit tracking; the first one given in (21) has the following discrete representation:(28)yψk+1ψk+1=TsVakcosθksinψk−ψt+yψk1−Tsτsψψk+0Tsτsψψck,and the second one given in (23) is represented in discrete time as(29)dk+1ψk+1=TsVakcosθkcosψk−φk+dk1−Tsτsψψk+0Tsτsψψck.The controllability of all discrete subsystems was revised.


### 2.3. Suboptimal Control for a Fixed-Wing UAV

With the general control law (13) presented in Section 2, the suboptimal nonlinear discrete control sequence for a UAV can be obtained. In the inner loop for control the pitch and roll angles consider the matrices Qθ and Qϕ as follows:(30)Qθ=Qθ11Qθ12Qθ21Qθ22,Qϕ=Qϕ11Qϕ12Qϕ21Qϕ22,with the scalars Rθ and Rϕ, the control laws are(31)τθN−k=−TsIyy2Qθ22+Rθ−1Ts2IyyQθ21+TsIyyQθ22θ2k+TsIyyQθ21θ1k+Ts2Iyy2Izz−IxxQθ22ψvkϕ2k,τϕN−k=−TsIxx2Qϕ22+Rϕ−1Ts2IxxQϕ21+TsIxxQϕ22ϕ2k+TsIxxQϕ21ϕ1k+Ts2Ixx2Iyy−IzzQϕ22ψvkθ2k,for pitch and roll, respectively. To control the system given in (25), consider the matrix Qψ defined as follows:(32)Qψ=Qψ11Qψ12Qψ21Qψ22,and the scalar Rψ; then, the control law for this system is(33)ϕcN−k=−Ts2τsϕ2Qψ22+Rψ−1Ts2τsϕgVakQψ21tanϕk+TsτsϕQψ21ψk+Tsτsϕ−Ts2τsϕ2Qψ22ϕk.The outer loop is divided into three subsystems; for subsystem (26), scalars QT and RT are used. Then, the control law is given by(34)TN−k=−Tsm2QT+RT−1−Tsm2QTD−Ts2mQTgsinθk+TsmQTVak.In order to the control altitude in subsystem (27), consider the matrix Qz defined as(35)Qz=Qz11Qz12Qz21Qz22.and the scalar Rz; then, the control law for this subsystem is(36)θcN−k=−Ts2τsθ2Qz22+Rz−1Ts2τsθQz21Vaksinθk+TsτsθQz21zk+Tsτsθ−Ts2τsθ2Qz22θk.To control the subsystem for straight-line tracking (28), consider the matrix Qy as(37)Qy=Qy11Qy12Qy21Qy22,and the scalar Ry; with this, the control law is given by(38)ψcN−k=−Ts2τsψ2Qy22+Ry−1Ts2τsψQy21Vakcosθksinψk−ψt+TsτsψQy21yψk+Tsτsψ−Ts2τsψ2Qy22ψk.Finally, to control the subsystem for circular orbit tracking (29), the matrix Qd is defined as(39)Qd=Qd11Qd12Qd21Qd22,and the scalar Rd; we have that the control law for this subsystem is(40)ψcN−k=−Ts2τsψ2Qd22+Rd−1Ts2τsψQd21Vakcosθkcosψk−φk+TsτsψQd21dk+Tsτsψ−Ts2τsψ2Qd22ψk.


## 2.1. Optimization Procedure of Finite Horizon

Consider the nonlinear discrete affine system(1)x¯k+1=f0x¯k+f1x¯kuk,where f0x¯k∈ℝn, f1x¯k∈ℝn×m, and uk∈ℝm. Now, the performance index is defined as in [12](2)J=12x¯NHx¯N+12∑k=0N−1x¯TkQx¯k+uTkRuk,where H,Q≥0 and R>0 are symmetric matrices of appropriate dimensions, the horizon is tf=TsN with Ts defining the sample period, N is the final step, and k∈ℕ, where k=0,1,…,N. The suboptimal control sequence is obtained as in [13] using the dynamic programming approach in discrete time proposed in [12]; in the last step, the performance index is(3)JN,N∗=12x¯NHx¯N.This term does not depend on the control lawuN, then we can say JN,N∗ is the optimal value of J in the discrete time N. For the next step, the performance index is(4)JN−1,N∗=minuN−1JN∗+12x¯TN−1Qx¯N−1+12uTN−1RuN−1.An explicit form ofJN−1,N∗ can be obtained by using (1)(5)JN−1,N∗=minuN−112f0x¯N−1+f1x¯N−1uN−1T×Hf0x¯N−1+f1x¯N−1uN−1+12x¯TN−1Qx¯N−1+12uTN−1RuN−1.To find the optimal value ofJN−1,N, it is calculated in the traditional way with the partial derivative of (5) with respect to uN−1 equal to zero.(6)u∗N−1=−f1Tx¯N−1Hf1x¯N−1+R−1f1Tx¯N−1Hf0x¯N−1.Equation (6) is the control law that minimizes JN−1,N in the discrete time N−1. With the condition set before R>0, we can be sure the inverse(7)f1Tx¯N−1Hf1x¯N−1+R−1,exists. With the purpose of verifying that (6) is minimum, JN−1,N is derived partially with respect to uN−1 for a second time; then, we have(8)∂2JN−1,N∂u2N−1=f1Tx¯N−1Hf1x¯N−1+R>0.The control given by equation (6) is the optimal control in local sense. For the next step N−2, the performance index is(9)JN−2,N∗=minuN−212x¯NHx¯N+12x¯TN−1Qx¯N−1+12uTN−1RuN−1+12x¯TN−2Qx¯N−2+12uTN−2RuN−2.Using (1) in (9), we have(10)JN−2,N∗=minuN−212f0x¯N−1+f1x¯N−1uN−1THf0x¯N−1+f1x¯N−1uN−1+12f0x¯N−2+f1x¯N−2uN−2TQf0x¯N−2+f1x¯N−2uN−2+12f0x¯N−3+f1x¯N−3uN−3TQf0x¯N−3+f1x¯N−3uN−3+12uTN−1RuN−1+12uTN−2RuN−2.Again we proceed in the traditional way to find the minimal value ofJN−2,N. In this case, the second partial derivative is derived with respect to uN−2; then, we have that(11)uN−2=−f1Tx¯N−2Qf1x¯N−2+R−1f1Tx¯N−2Qf0x¯N−2.Repeating this task for a few more steps, we realized that the general expressions are as follows.Fork=1,(12)uN−k=−f1Tx¯N−kHf1x¯N−k+R−1f1Tx¯N−kHf0x¯N−k,and, for k=2,…,N,(13)uN−k=−f1Tx¯N−kQf1x¯N−k+R−1f1Tx¯N−kQf0x¯N−k.This suboptimal control sequence guarantees that an approximation to the minimal value of the performance index (2) is reached.


## 2.2. Mathematical Model of a Fixed-Wing UAV

The previously exposed control strategy is implemented in a three-layer control structure composed of an inner loop, a middle loop, and an outer loop.


### 2.2.1. Continuous Nonlinear Model

For the inner loop, we used a part of the dynamic model presented in [14, 15], where due to the symmetry of the UAV Ixz is very much smaller than Ixx, Iyy, and Izz and can be neglected:(14)θ..=Izz−IxxIyyψvϕ˙+1Iyyτθ,ϕ¨=Iyy−IzzIxxψvθ.+1Ixxτϕ,where Ixx, Iyy, and Izz are the inertia parameters, ψv is the velocity in angle yaw, θ is the pitch angle, ϕ is the roll angle, and τθ and τϕ are control inputs for pitch and roll, respectively. The purpose of this inner loop is to follow the commands given by the middle loop and outer loop.The middle loop is dedicated to control the course angleχ. As it is described in [16], the course angle χ is the direction of travel relative to the axis x, the yaw or heading angle ψ is the direction that the UAV is pointed, and the side-slip angle β is the difference between the course and the heading angles. Figure 1 depicts this.Figure 1
Definition of the headingψ, side-slip β, and course χ angles.In this paper, we consider coordinated turn conditions, absence of wind, and side-slip. Then, we have thatψ=χ; with this, the kinematic model is given by(15)ψ˙=gVatanϕ,ϕ˙=1τsϕϕc−ϕ,where g is the magnitude of gravity at the sea level, τsϕ is the constant time of the system, ϕc is the control angle commanded, and Va is the airspeed which under the assumption of absence of wind this value is the UAV ground speed Va=Vg.Finally, for the outer loop, we considered the next kinematic model for guidance based on the one proposed in [16](16)x˙=Vacosψcosθ,y˙=Vasinψcosθ,z˙=Vasinθ,Va.=Tm−Dm−gsinθ,ψ˙=1τsψψc−ψ,θ.=1τsθθc−θ,where x, y, and z denote the inertial position of the aircraft in a 3D inertial frame, m is the mass of the UAV, D represents the drag force, ψ is the heading angle (or course angle under the assumptions of this paper), θ is the pitch angle, and τsθ and τsψ are constants time of the system. In this kinematic model, we have 3 control inputs: T is the thrust motor and θc and ψc are commanded angles. In this paper, we use the course angle for trajectory tracking, but with the suppositions made before we have that ψ=χ.In order to simplify the process to obtain the control expressions for the aircraft, the model is divided into 3 subsystems as follows:SubsystemVa:(17)Va.=Tm−Dm−gsinθ.This subsystem is dedicated to control the velocity of the aircraftVa and its control input is the thrust force T.Subsystemz:(18)z˙=Vasinθ,θ.=1τsθθc−θ.This subsystem is for control the altitude position, and the control input isθc.Subsystemxy:(19)x˙=Vacosψcosθ,y˙=Vasinψcosθ,ψ˙=1τsψψc−ψ.This subsystem is devoted to controlling the position in the planex−y and has as a control input ψc.In the trajectory tracking problem, there are basically two kinds of paths to follow, straight-line segments and circular orbits; with a combination of these, more complex paths can be tracked. Then, it is necessary to represent the model given in (19) in a form where we can know the cross-track error in an easy way.For straight-line segments, we used the matrix rotation given in (20) to transform the inertial frame to the path frame:(20)Rline=cosψtsinψt−sinψtcosψt,where ψt is the course angle of the path. Using (20), we calculated a new representation of (19), where the important dynamics are(21)y˙ψ=Vacosθsinψ−ψt,ψ˙=1τsψψc−ψ.For the circular orbit tracking problem, due to its nature, it is better to represent (19) in a polar coordinate system; for this, we use the matrix rotation given by(22)Rorbit=cosφsinφ−sinφcosφ,where φ is the phase angle of the aircraft. Using (22), a new representation of (19) is calculated, where the important dynamics are(23)d˙=Vacosθcosψ−φ,ψ˙=1τsψψc−ψ,where d is the distance from the UAV to the center of the coordinated system.


### 2.2.2. Discrete Nonlinear Model

The nonlinear discrete affine representation of the continuous nonlinear models was obtained by using the Euler approximation. For the inner loop, the nonlinear discrete affine state-space representation of the dynamic model (14) is(24)θ1k+1θ2k+1=Tsθ2k+θ1kTsIzz−IxxIyyψvkϕ2k+θ2k+0TsIyyτθk,ϕ1k+1ϕ2k+1=Tsϕ2k+ϕ1kTsIyy−IzzIxxψvkθ2k+ϕ2k+0TsIxxτϕk,respectively, where Ts is the sample time. The nonlinear discrete affine representation of (15) is given by(25)ψk+1ϕk+1=TsgVaktanϕk+ψk1−Tsτsϕϕk+0Tsτsϕϕck.The subsystemVa given in (17) has the discrete representation as follows:(26)Vak+1=−TsmD−Tsgsinθk+Vak+TsmTk.For subsystem (18), its nonlinear discrete affine representation is(27)zk+1θk+1=TsVaksinθk+zk1−Tsτsθθk+0Tsτsθθck.The subsystemxy has two versions: one for straight-line trajectory tracking and the other one for circular orbit tracking; the first one given in (21) has the following discrete representation:(28)yψk+1ψk+1=TsVakcosθksinψk−ψt+yψk1−Tsτsψψk+0Tsτsψψck,and the second one given in (23) is represented in discrete time as(29)dk+1ψk+1=TsVakcosθkcosψk−φk+dk1−Tsτsψψk+0Tsτsψψck.The controllability of all discrete subsystems was revised.


## 2.2.1. Continuous Nonlinear Model

For the inner loop, we used a part of the dynamic model presented in [14, 15], where due to the symmetry of the UAV Ixz is very much smaller than Ixx, Iyy, and Izz and can be neglected:(14)θ..=Izz−IxxIyyψvϕ˙+1Iyyτθ,ϕ¨=Iyy−IzzIxxψvθ.+1Ixxτϕ,where Ixx, Iyy, and Izz are the inertia parameters, ψv is the velocity in angle yaw, θ is the pitch angle, ϕ is the roll angle, and τθ and τϕ are control inputs for pitch and roll, respectively. The purpose of this inner loop is to follow the commands given by the middle loop and outer loop.The middle loop is dedicated to control the course angleχ. As it is described in [16], the course angle χ is the direction of travel relative to the axis x, the yaw or heading angle ψ is the direction that the UAV is pointed, and the side-slip angle β is the difference between the course and the heading angles. Figure 1 depicts this.Figure 1
Definition of the headingψ, side-slip β, and course χ angles.In this paper, we consider coordinated turn conditions, absence of wind, and side-slip. Then, we have thatψ=χ; with this, the kinematic model is given by(15)ψ˙=gVatanϕ,ϕ˙=1τsϕϕc−ϕ,where g is the magnitude of gravity at the sea level, τsϕ is the constant time of the system, ϕc is the control angle commanded, and Va is the airspeed which under the assumption of absence of wind this value is the UAV ground speed Va=Vg.Finally, for the outer loop, we considered the next kinematic model for guidance based on the one proposed in [16](16)x˙=Vacosψcosθ,y˙=Vasinψcosθ,z˙=Vasinθ,Va.=Tm−Dm−gsinθ,ψ˙=1τsψψc−ψ,θ.=1τsθθc−θ,where x, y, and z denote the inertial position of the aircraft in a 3D inertial frame, m is the mass of the UAV, D represents the drag force, ψ is the heading angle (or course angle under the assumptions of this paper), θ is the pitch angle, and τsθ and τsψ are constants time of the system. In this kinematic model, we have 3 control inputs: T is the thrust motor and θc and ψc are commanded angles. In this paper, we use the course angle for trajectory tracking, but with the suppositions made before we have that ψ=χ.In order to simplify the process to obtain the control expressions for the aircraft, the model is divided into 3 subsystems as follows:SubsystemVa:(17)Va.=Tm−Dm−gsinθ.This subsystem is dedicated to control the velocity of the aircraftVa and its control input is the thrust force T.Subsystemz:(18)z˙=Vasinθ,θ.=1τsθθc−θ.This subsystem is for control the altitude position, and the control input isθc.Subsystemxy:(19)x˙=Vacosψcosθ,y˙=Vasinψcosθ,ψ˙=1τsψψc−ψ.This subsystem is devoted to controlling the position in the planex−y and has as a control input ψc.In the trajectory tracking problem, there are basically two kinds of paths to follow, straight-line segments and circular orbits; with a combination of these, more complex paths can be tracked. Then, it is necessary to represent the model given in (19) in a form where we can know the cross-track error in an easy way.For straight-line segments, we used the matrix rotation given in (20) to transform the inertial frame to the path frame:(20)Rline=cosψtsinψt−sinψtcosψt,where ψt is the course angle of the path. Using (20), we calculated a new representation of (19), where the important dynamics are(21)y˙ψ=Vacosθsinψ−ψt,ψ˙=1τsψψc−ψ.For the circular orbit tracking problem, due to its nature, it is better to represent (19) in a polar coordinate system; for this, we use the matrix rotation given by(22)Rorbit=cosφsinφ−sinφcosφ,where φ is the phase angle of the aircraft. Using (22), a new representation of (19) is calculated, where the important dynamics are(23)d˙=Vacosθcosψ−φ,ψ˙=1τsψψc−ψ,where d is the distance from the UAV to the center of the coordinated system.


## 2.2.2. Discrete Nonlinear Model

The nonlinear discrete affine representation of the continuous nonlinear models was obtained by using the Euler approximation. For the inner loop, the nonlinear discrete affine state-space representation of the dynamic model (14) is(24)θ1k+1θ2k+1=Tsθ2k+θ1kTsIzz−IxxIyyψvkϕ2k+θ2k+0TsIyyτθk,ϕ1k+1ϕ2k+1=Tsϕ2k+ϕ1kTsIyy−IzzIxxψvkθ2k+ϕ2k+0TsIxxτϕk,respectively, where Ts is the sample time. The nonlinear discrete affine representation of (15) is given by(25)ψk+1ϕk+1=TsgVaktanϕk+ψk1−Tsτsϕϕk+0Tsτsϕϕck.The subsystemVa given in (17) has the discrete representation as follows:(26)Vak+1=−TsmD−Tsgsinθk+Vak+TsmTk.For subsystem (18), its nonlinear discrete affine representation is(27)zk+1θk+1=TsVaksinθk+zk1−Tsτsθθk+0Tsτsθθck.The subsystemxy has two versions: one for straight-line trajectory tracking and the other one for circular orbit tracking; the first one given in (21) has the following discrete representation:(28)yψk+1ψk+1=TsVakcosθksinψk−ψt+yψk1−Tsτsψψk+0Tsτsψψck,and the second one given in (23) is represented in discrete time as(29)dk+1ψk+1=TsVakcosθkcosψk−φk+dk1−Tsτsψψk+0Tsτsψψck.The controllability of all discrete subsystems was revised.


## 2.3. Suboptimal Control for a Fixed-Wing UAV

With the general control law (13) presented in Section 2, the suboptimal nonlinear discrete control sequence for a UAV can be obtained. In the inner loop for control the pitch and roll angles consider the matrices Qθ and Qϕ as follows:(30)Qθ=Qθ11Qθ12Qθ21Qθ22,Qϕ=Qϕ11Qϕ12Qϕ21Qϕ22,with the scalars Rθ and Rϕ, the control laws are(31)τθN−k=−TsIyy2Qθ22+Rθ−1Ts2IyyQθ21+TsIyyQθ22θ2k+TsIyyQθ21θ1k+Ts2Iyy2Izz−IxxQθ22ψvkϕ2k,τϕN−k=−TsIxx2Qϕ22+Rϕ−1Ts2IxxQϕ21+TsIxxQϕ22ϕ2k+TsIxxQϕ21ϕ1k+Ts2Ixx2Iyy−IzzQϕ22ψvkθ2k,for pitch and roll, respectively. To control the system given in (25), consider the matrix Qψ defined as follows:(32)Qψ=Qψ11Qψ12Qψ21Qψ22,and the scalar Rψ; then, the control law for this system is(33)ϕcN−k=−Ts2τsϕ2Qψ22+Rψ−1Ts2τsϕgVakQψ21tanϕk+TsτsϕQψ21ψk+Tsτsϕ−Ts2τsϕ2Qψ22ϕk.The outer loop is divided into three subsystems; for subsystem (26), scalars QT and RT are used. Then, the control law is given by(34)TN−k=−Tsm2QT+RT−1−Tsm2QTD−Ts2mQTgsinθk+TsmQTVak.In order to the control altitude in subsystem (27), consider the matrix Qz defined as(35)Qz=Qz11Qz12Qz21Qz22.and the scalar Rz; then, the control law for this subsystem is(36)θcN−k=−Ts2τsθ2Qz22+Rz−1Ts2τsθQz21Vaksinθk+TsτsθQz21zk+Tsτsθ−Ts2τsθ2Qz22θk.To control the subsystem for straight-line tracking (28), consider the matrix Qy as(37)Qy=Qy11Qy12Qy21Qy22,and the scalar Ry; with this, the control law is given by(38)ψcN−k=−Ts2τsψ2Qy22+Ry−1Ts2τsψQy21Vakcosθksinψk−ψt+TsτsψQy21yψk+Tsτsψ−Ts2τsψ2Qy22ψk.Finally, to control the subsystem for circular orbit tracking (29), the matrix Qd is defined as(39)Qd=Qd11Qd12Qd21Qd22,and the scalar Rd; we have that the control law for this subsystem is(40)ψcN−k=−Ts2τsψ2Qd22+Rd−1Ts2τsψQd21Vakcosθkcosψk−φk+TsτsψQd21dk+Tsτsψ−Ts2τsψ2Qd22ψk.


## 3. Results

In this section, simulations, results, and a comparative study with a PD controller are presented. Simulations were developed using the software Matlab/Simulink version 9.9 (R2020b) and X-Plane 11.


### 3.1. Simulation Results

Two simulations were performed: the first one using the suboptimal nonlinear discrete control laws developed in this work and the second one using PD controllers with their parameters overshot and settling time heuristically tuned. Both controllers were tested under the same conditions; the task to perform was to track an orbit trajectory inside two thermal updrafts, and the simulation lasted 1000 seconds with initial positionx=0m, y=0m, and z=80m. At the beginning, the reference on z was 200 m, but after the second 120, there was not any reference on z in order to be able to gain altitude with the thermal updrafts. In the first 80 seconds, the plane takes off and then it starts to do the task, and the sampling time was 20 ms.Thermal updrafts were set with a diameter of 1500 m and their centers on the positionx=5000m, y=1150m for the first one and x=78000m, y=−4100m for the second one, using a plugin for the simulator. Wind disturbances were set with the intensity of 2 knots in directions south, east, and west.In order to implement the controller described in Section3, the following penalization matrices were considered.For orientation,(41)Hθ=Qθ=Hϕ=Qϕ=50000000200002000001000,Rθ=Rϕ=0.3,Hψ=Qψ=2500015150.01,Rψ=1.2.For subsystemVa,(42)HT=QT=2700,RT=0.1.For the subsystemz,(43)Hz=Qz=20101010Rz=5.For subsystemxy,(44)Hy=Qy=950030300.01Ry=0.5.The position of the plane in the simulations using the suboptimal nonlinear discrete control is shown in Figure2, and the position with the PD control in Figure 3; in these figures, we can see that the altitude at the end of the simulation was higher using the suboptimal control due to the angle of roll during the trajectory tracking, which means we can exploit in a better way the thermal updrafts using the suboptimal control. The final altitude of the plane with the suboptimal control was 1071 m and with the PD was 719 m.Figure 2
Position of the plane with the suboptimal nonlinear discrete control.Figure 3
Position of the plane with the PD control.The angular orientation is shown in Figures4 and 5 for the suboptimal control and the PD, respectively; the average angle of roll during the orbit tracking inside a thermal updraft for the suboptimal control was 15 degrees, and for the PD was 18 degrees; from this difference on the roll angle, we have a different final altitude, and the area of the plane affected for the thermal updraft is greater using the suboptimal controller. This means more lift from the thermal updraft and then more altitude. It is clear that using the suboptimal controller we can exploit in a better way thermal updrafts.Figure 4
Orientation of the plane with the suboptimal nonlinear discrete control.Figure 5
Orientation of the plane with the PD control.The trajectory of the plane for both controllers is shown in Figures6 and 7. In these figures, we can see the position of the thermals too; due to the scale of the figures, it is difficult to see the differences in the results. The first 80 seconds are for the take-off phase and after that, the plane starts to follow the trajectory; that is the reason why the angle of course does not change since the beginning.Figure 6
Trajectory of the plane with the suboptimal nonlinear discrete control.Figure 7
Trajectory of the plane with the PD control.Control signals are represented in values from −1 to 1 in the case of elevator and ailerons which means the grade of deflection of the control surface and from 0 to 1 in the case of thrust which means the power of the motor, 0 means the motor is off and 1 means motor in its maximum power. Control signals for the suboptimal control are shown in Figure8 and for the PD control in Figure 9.Figure 8
Control signals with the suboptimal nonlinear discrete control.Figure 9
Control signals with the PD control.The desired airspeed was 44 m/s; in Figures10 and 11, we can see the performance of the suboptimal control and the PD control, respectively.Figure 10
Airspeed with the suboptimal nonlinear discrete control.Figure 11
Airspeed with the PD control.From the last two figures, we can see that the error in the desired airspeed is less using the suboptimal control law than using the PD control, and it should be noted that there is a little stationary error in the PD controller. Finally, the tracking error in the trajectory is shown in Figures12 and 13; again due to the scale, it is difficult to see the details. The tracking error was less with the suboptimal control law than the PD control.Figure 12
Trajectory tracking error with the suboptimal nonlinear discrete control.Figure 13
Trajectory tracking error with the PD control.


### 3.2. Comparative Study

In order to know which controller presented the best performance and energy consumption, a comparative study was made. To compute the energy consumption the index,Ju=∑uk was used, where uk is the control signal in step k, and the index used for the performance was Je=∑ek, where ek is the error in step k.Table1 shows the energy consumption of elevator, ailerons, and motor during this tracking trajectory task; as we can see, the energy consumption of the elevator is similar for both controllers, but in the case of aileron there was less energy consumption using the suboptimal control, and in the case of the motor it was the suboptimal control the one using less energy too; it is important to save energy in this element because the motor is the one that consumes much more significant energy during a flight than the rest of electronic components, according to [17] in the order of kilowatts versus watts.Table 1
Comparative of saving energy.
ControlSuboptimalJuPDJuElevator4.8608 × 1034.9482 × 103Ailerons2.2172 × 1034.5410 × 103Motor2.7591 × 1043.1181 × 104Results in performance are shown in Table2; both controllers, suboptimal and PD, had a similar performance tracking the trajectory, but in the case of the desired reference in airspeed the one with better performance was the suboptimal control law. This better performance is because the suboptimal control law uses the nonlinear model that represents more precise dynamics than the linear model which was used for the PD.Table 2
Comparative of performance.
ControlSuboptimalJePDJeAirspeed2.8078 × 1043.8336 × 104Trajectory17.7597 × 10618.2099 × 106


## 3.1. Simulation Results

Two simulations were performed: the first one using the suboptimal nonlinear discrete control laws developed in this work and the second one using PD controllers with their parameters overshot and settling time heuristically tuned. Both controllers were tested under the same conditions; the task to perform was to track an orbit trajectory inside two thermal updrafts, and the simulation lasted 1000 seconds with initial positionx=0m, y=0m, and z=80m. At the beginning, the reference on z was 200 m, but after the second 120, there was not any reference on z in order to be able to gain altitude with the thermal updrafts. In the first 80 seconds, the plane takes off and then it starts to do the task, and the sampling time was 20 ms.Thermal updrafts were set with a diameter of 1500 m and their centers on the positionx=5000m, y=1150m for the first one and x=78000m, y=−4100m for the second one, using a plugin for the simulator. Wind disturbances were set with the intensity of 2 knots in directions south, east, and west.In order to implement the controller described in Section3, the following penalization matrices were considered.For orientation,(41)Hθ=Qθ=Hϕ=Qϕ=50000000200002000001000,Rθ=Rϕ=0.3,Hψ=Qψ=2500015150.01,Rψ=1.2.For subsystemVa,(42)HT=QT=2700,RT=0.1.For the subsystemz,(43)Hz=Qz=20101010Rz=5.For subsystemxy,(44)Hy=Qy=950030300.01Ry=0.5.The position of the plane in the simulations using the suboptimal nonlinear discrete control is shown in Figure2, and the position with the PD control in Figure 3; in these figures, we can see that the altitude at the end of the simulation was higher using the suboptimal control due to the angle of roll during the trajectory tracking, which means we can exploit in a better way the thermal updrafts using the suboptimal control. The final altitude of the plane with the suboptimal control was 1071 m and with the PD was 719 m.Figure 2
Position of the plane with the suboptimal nonlinear discrete control.Figure 3
Position of the plane with the PD control.The angular orientation is shown in Figures4 and 5 for the suboptimal control and the PD, respectively; the average angle of roll during the orbit tracking inside a thermal updraft for the suboptimal control was 15 degrees, and for the PD was 18 degrees; from this difference on the roll angle, we have a different final altitude, and the area of the plane affected for the thermal updraft is greater using the suboptimal controller. This means more lift from the thermal updraft and then more altitude. It is clear that using the suboptimal controller we can exploit in a better way thermal updrafts.Figure 4
Orientation of the plane with the suboptimal nonlinear discrete control.Figure 5
Orientation of the plane with the PD control.The trajectory of the plane for both controllers is shown in Figures6 and 7. In these figures, we can see the position of the thermals too; due to the scale of the figures, it is difficult to see the differences in the results. The first 80 seconds are for the take-off phase and after that, the plane starts to follow the trajectory; that is the reason why the angle of course does not change since the beginning.Figure 6
Trajectory of the plane with the suboptimal nonlinear discrete control.Figure 7
Trajectory of the plane with the PD control.Control signals are represented in values from −1 to 1 in the case of elevator and ailerons which means the grade of deflection of the control surface and from 0 to 1 in the case of thrust which means the power of the motor, 0 means the motor is off and 1 means motor in its maximum power. Control signals for the suboptimal control are shown in Figure8 and for the PD control in Figure 9.Figure 8
Control signals with the suboptimal nonlinear discrete control.Figure 9
Control signals with the PD control.The desired airspeed was 44 m/s; in Figures10 and 11, we can see the performance of the suboptimal control and the PD control, respectively.Figure 10
Airspeed with the suboptimal nonlinear discrete control.Figure 11
Airspeed with the PD control.From the last two figures, we can see that the error in the desired airspeed is less using the suboptimal control law than using the PD control, and it should be noted that there is a little stationary error in the PD controller. Finally, the tracking error in the trajectory is shown in Figures12 and 13; again due to the scale, it is difficult to see the details. The tracking error was less with the suboptimal control law than the PD control.Figure 12
Trajectory tracking error with the suboptimal nonlinear discrete control.Figure 13
Trajectory tracking error with the PD control.


## 3.2. Comparative Study

In order to know which controller presented the best performance and energy consumption, a comparative study was made. To compute the energy consumption the index,Ju=∑uk was used, where uk is the control signal in step k, and the index used for the performance was Je=∑ek, where ek is the error in step k.Table1 shows the energy consumption of elevator, ailerons, and motor during this tracking trajectory task; as we can see, the energy consumption of the elevator is similar for both controllers, but in the case of aileron there was less energy consumption using the suboptimal control, and in the case of the motor it was the suboptimal control the one using less energy too; it is important to save energy in this element because the motor is the one that consumes much more significant energy during a flight than the rest of electronic components, according to [17] in the order of kilowatts versus watts.Table 1
Comparative of saving energy.
ControlSuboptimalJuPDJuElevator4.8608 × 1034.9482 × 103Ailerons2.2172 × 1034.5410 × 103Motor2.7591 × 1043.1181 × 104Results in performance are shown in Table2; both controllers, suboptimal and PD, had a similar performance tracking the trajectory, but in the case of the desired reference in airspeed the one with better performance was the suboptimal control law. This better performance is because the suboptimal control law uses the nonlinear model that represents more precise dynamics than the linear model which was used for the PD.Table 2
Comparative of performance.
ControlSuboptimalJePDJeAirspeed2.8078 × 1043.8336 × 104Trajectory17.7597 × 10618.2099 × 106


## 4. Conclusions

This paper has considered the problem of flight time endurance in fixed-wing UAVs, and it was proposed to use a suboptimal nonlinear discrete control to save energy and improve performance, in combination with navigation in thermal updrafts. Control equations were obtained for dynamic and kinematic nonlinear models of a fixed-wing UAV and tested in simulations; these simulation results showed that this control technique consumes less energy than a PD control and has similar performance. Also, using this suboptimal nonlinear discrete control to navigate in thermal updrafts helps to exploit in a better way the updraft gaining more altitude than the PD controller.



---

*Source: 2214217-2022-03-12.xml*