Robust Nonlinear Tracking Control for Unmanned Aircraft in the Presence of Wake Vortex

: The ﬂight trajectory of unmanned aerial vehicles (UAVs) can be signiﬁcantly affected by external disturbances such as turbulence, upstream wake vortices, or wind gusts. These effects present challenges for UAV ﬂight safety. Hence, addressing these challenges is of critical importance for the integration of unmanned aerial systems (UAS) into the National Airspace System (NAS), especially in terminal zones. This work presents a robust nonlinear control method that has been designed to achieve roll/yaw regulation in the presence of unmodeled external disturbances and system nonlinearities. The data from NASA-conducted airport experimental measurements as well as high-ﬁdelity Large Eddy Simulations of the wake vortex are used in the study. Side-by-side simulation comparisons between the robust nonlinear control law and both linear H ∞ and PID control laws are provided for completeness. These simulations are focused on applications involving small UAV affected by the wake vortex disturbance in the vicinity of the ground (which models the take-off or landing phase) as well as in the out-of-ground zone. The results demonstrate the capability of the proposed nonlinear controller to asymptotically reject wake vortex disturbance in the presence of the nonlinearities in the system (i.e., parametric variations, unmodeled, time-varying disturbances). Further, the nonlinear controller is designed with a computationally efﬁcient structure without the need for the complex calculations or function approximators in the control loop. Such a structure is motivated by UAV applications where onboard computational resources are limited.


Introduction
The Federal Aviation Administration (FAA) currently faces numerous operational safety challenges associated with the integration of UAS flights in the NAS. For instance, the Integrated Safety Assessment Model (ISAM) [1] requires further development and improvements in order to adequately address UAS operations for current and future risk analyses. Motivated by the desire to improve the safety of UAS operating in the NAS, the development of novel UAV flight control technologies is of critical importance. Specifically, there is a need for control system technologies that are capable of quickly recovering from unpredictable and potentially hazardous operating conditions resulting from phenomena such as airflow disturbances due to upstream wake vortex, wind gusts, or turbulence. Based on these considerations, the focus of the current work is on the development of a nonlinear control method, which demonstrates a positive and accurate UAV trajectory regulation in the presence of wake vortex disturbance, in addition to the nonlinearity in the governing UAV dynamic model.
The UAV is a multi-input/multi-output (MIMO) system, which is a highly coupled and unstable system. A number of works are devoted to the problem of the UAV following the path in a windy environment [2,3] as well as in the presence of microburst [4]. Different control strategies have been widely used for the UAVs' control, such as adaptive control [5][6][7], neural network-based [8,9] techniques, inversion techniques [10,11], optimal control [12,13], which may effectively compensate for the actuator nonlinearities. The linear quadratic regulator (LQR) [14,15] as well as robust H ∞ [16][17][18] control methods have been used for lateral and longitudinal dynamics in flight control. However, such techniques require additional computational complexity over purely robust feedback designs. Hence, the minimalism of the controller design in this work is motivated by the desire to develop control methods that are suitable for small UAVs.
This work presents a robust nonlinear flight control strategy [19][20][21] capable of the uncertain bounded wake vortex disturbance rejection in different phases of evolution, including in-ground and out-of-ground effects. The control of the systems in the presence of uncertainties and unmodeled disturbances is used nowadays for various practical applications. For example, the robust control as well as station-keeping control of the underwater vehicle in the presence of the unmodeled disturbances and model uncertainties was applied by Vu et al. [22,23] to obtain the desired trajectory tracking performance of the AUV. Haitong Xu [24] recently demonstrated the control system for coupled pathfollowing and the collision-avoided task of the ship model in the presence of static obstacles. The sliding mode control for second-order systems in application to UAVs for attitude and altitude control has been described [25,26]. The effective application of the robust control for the second order-systems [27,28] in the presence of the uncertain bounded disturbances was shown. The disturbance is modeled by applying the low-fidelity model of wake vortex pair [29,30] for the steady level flight and Monte-Carlo simulations scenario [31]. The high fidelity simulations (Large Eddy Simulation (LES)) of the wake vortex for the generator aircraft's take-off and landing phases [32][33][34] are performed to model the in-ground phase of vortex evolution. Side-by-side simulations and performance comparisons with H ∞ control method are provided for the robust nonlinear controller response at every flight phase. In addition, the comparison with the PID controller is demonstrated. The proposed control design is particularly advantageous in maintaining flight stability in the presence of a high degree of uncertainty and nonlinearity in the UAV operating conditions (e.g., flight conditions inherent in tight urban environments and terminal zones).
The remainder of this paper is organized as follows. The next three sections present the mathematical model utilized in the controller development as well as the nonlinear and linear controllers considered. Section 6 presents the wake vortex disturbance modeling methodology for low-fidelity and high-fidelity simulations. Finally, Section 7 outlines the numerical simulation and results, where the performance of the proposed controller is investigated and compared with the existing linear controllers.

Mathematical Model
This section describes the mathematical model utilized in controller development. The UAV dynamic model under consideration is assumed to contain parametric model uncertainty in addition to unmodeled, time-varying nonlinearities. The aircraft dynamics is modeled as a linear parameter-varying (LPV) system [35][36][37][38][39][40][41]: where ρ is a measurable parameter vector, called the scheduling parameter, x ∈ R n are the deviation states, u ∈ R m are the deviation inputs, y ∈ R p are the deviation outputs, and f (x, t) is an unknown nonlinear disturbance. Here, A(ρ) ∈ R n×n represents the system matrix, B(ρ) ∈ R n×m the input matrix, C(ρ) ∈ R p×n the output matrix, and D(ρ) ∈ R p×m the feedthrough matrix. The system is then arranged such that where A 0 , B 0 , C 0 , and D 0 are the known, nominal state space matrices. The parametric uncertainty is reflected by δ i ∈ [−1, 1], and the structural knowledge about the uncertainty is contained in the matrices A i , B i , C i , and D i [42]. The lateral dynamics of an aircraft is described by the state vector x = [v, p, r, φ, ψ] T with input vector u = [δ a , δ r ] T , where v,p,r,φ,ψ is lateral velocity, roll, and yaw rates and angles. The input parameters are aileron and rudder deflections. The maximum rate of deflection of ailerons and rudder is set to 30 deg/s, and the maximum amplitude is 25 deg. The scheduling parameter ρ can be the altitude and Mach number, and the unknown disturbance f (x, t) could be wake vortex, wind gusts, turbulent disturbances, or nonlinearities not captured in the linear model.

•
Assumption 1: We assume that the lateral dynamics is uncoupled from the longitudinal. • Assumption 2: It is assumed that the magnitude of the wake vortex disturbance is bounded.

Control Objective
The control objective is to force the UAV attitude (i.e., φ(t) and ψ(t)) to regulate to a given constant value. To quantify the control objective, the trajectory regulation error e(t) ∈ R 2 and the auxiliary regulation error r(t) = r p (t) r r (t) T ∈ R 2 are defined as: In Equation (7), α 1 , α 2 ∈ R denote positive, constant control gains. Thus, the trajectory regulation control objective can be stated mathematically as e(t) → 0, where * denotes the standard Euclidean norm of the vector argument. Note that, based on the auxiliary regulation error definitions in Equation (7), r(t) → 0 ⇒ e(t) → 0.
Note that the regulation error and auxiliary errors defined in Equations (6) and (7) are key in the contribution in the current study. The definitions of the auxiliary regulation errors enable us to recast the dynamic model in Equation (1) in a form that is amenable to roll/yaw angle angle regulation control. Indeed, it can be seen that the differentiation of r(t) produces a set of equations that render the roll and yaw angles (φ, ψ) fully controllable through the aileron and rudder deflections δ a (t), δ r (t) . Thus, the auxiliary error terms r r (t), r p (t) can be viewed as sliding surfaces, which enables us to prove our roll and yaw angle regulation results.

Robust Controller Development
In order to achieve asymptotic convergence of φ, ψ to zero with a given convergence rate in the presence of a bounded disturbance (i.e., wake vortex or wind gust), we have to drive the auxiliary regulation error r to zero in finite time. Taking into account the original dynamic model and the auxiliary regulation error r, the auxiliary control term u is designed as: whereΩ denotes a constant auxiliary matrix, and [] −1 denotes the inverse of a matrix. The feedback control gains (i.e., amplifiers) α 1 ,α 2 , k 1 ,k 2 ,β 1 ,β 2 can be tuned to adjust the closed loop regulation response to achieve the desired system performance (e.g., to achieve a faster response time). Note that the robust nonlinear control law given by Equation (8) ensures asymptotic trajectory regulation in the sense that provided the control gains k 1 and k 2 introduced in Equation (8) are selected sufficiently large, and α 1 and α 2 are selected based on the known upper bounds on the disturbance (i.e., the known maximum velocity and acceleration of the wind gust). A straightforward proof of the above statement based on Lyapunov stability analysis can be utilized and is omitted here for brevity. Details of the proof are similar to those provided in our recent results [19,20].

Observer Design
The control design presented in the current result is based on the assumption that only attitude measurements (i.e., φ(t) and ψ(t)) are available for feedback. Considering this, an observer (or estimator) will be employed, which estimates the complete system state using only the available attitude measurements. This section summarizes the design of the observer.
In Equation (1), the explicitly defined output y can be denoted as the sufficiently differentiable vector function h(x). To facilitate the subsequent observer design and analysis, a vector H(x) ∈ R n of output derivatives is defined as: where L f i h(x) denotes the i-th Lie derivative of the output function h(x) along the direction of the vector field Ax. An observer that estimates the full state x(t) of the system in Equation (1) using only measurements of y(t) can be designed as [21] In Equation (11), (11), M(x) ∈ R n×n denotes a diagonal matrix with positive elements defined as: where m i (x) ∈ R,x ∈ R n i = 1, . . . , n. Similar to the design in [43], it can be shown that the observer design in Equation (11) achieves finite-time estimation of the state x(t). Specifically, it can be shown that, through judicious selection of the diagonal matrix M(x),x(t) = x(t) for any t ≥ t 1 . By including the additional term in the observer design, it follows that, for t ≥ t 1 , the system converges to the sliding manifold σ = V(t) − H(x) , and the observer Equation (11) exactly estimates the state x(t) of the dynamic system in Equation (1).

H ∞ Linear Controller
The interconnection used to synthesize the linear controllers is depicted in Figure 1.
Here, e represents the control objective outputs, i.e., the measurable attitude outputs φ(t) and ψ(t). The controller K ∞ takes the estimated state vector as feedback measurements. The objective is to minimize the H ∞ -norm between the disturbances and wind gust d, w and the lateral attitude and control action e, u while compensating for model uncertainty. The performance specifications in Figure 1 are in the form of closed-loop transfer functions designed to be small through feedback. With this, the mathematical objective is to drive the H ∞ -norm of the closed-loop transfer function, γ, to be less than 1. Hence, the design task is to find the appropriate weighting functions to meet the performance objectives. For this control design, the closed-loop transfer functions of the aircraft are where Constant weights, W m , are used to model the disturbances to each control input. This weight is selected as W m = 0.4 for all control inputs. Furthermore, a multiplicative input weight is included to avoid destabilizing unmodeled frequency modes outside the control bandwidth. The uncertainty model weight is described by W u = 1.22(s + 0.48)/(s + 0.224).
On the other hand, the disturbance rejection of critical variables is enforced by a constant performance weight W p . A constant gain attenuation of W p = 1 is assigned to each control objective in order to ensure damping. Similarly, wind gust disturbances are attenuated with a constant weight W d = 1, and noise is accounted for with the constant weight W n = 0.05 for each measurement. Once all the weights are defined, the next step is to find an optimal K ∞ such that the system gain ||F CL (G p , K ∞ )|| ∞ < 1. The Robust Control Toolbox from MATLAB is used to design the H ∞ controller with a system gain of γ = 0.99. Hence, the designed controller achieves the control objectives specified in the control interconnection.

PID Controller
To provide a further comparison of the closed-loop system under wake vortex disturbances, the roll and yaw (i.e., φ(t) and ψ(t)) regulation control objective was tested using a PID controller. The Control Toolbox from MATLAB is used to design a PID controller. To adapt the Matlab PID controller simulation to our MIMO system, PID control laws are applied to each of the roll and yaw states separately. The parameters for proportional (P), integral (I) and derivative (D) gains, as well as filter coefficient (N) used in the transfer function for the PID controller Equation (16), are shown in Table 1. The Matlab "Transfer Function" tuning method was used to tune the controller.

Wake Vortex Modeling
In order to determine the loads induced on the aircraft flying through the wake vortex, one needs to model the flow field induced by the vortex pair as well as the forces and moments acting on the airplane during the interaction. By the interaction, the authors mean the wake-vortex/aircraft contact when sweeping the follower aircraft (FA) through the wake in the direction perpendicular to the generator aircraft (GA) velocity vector, as shown in Figure 2. The general block diagram is shown in Figure 3. Knowing the wake characteristics of the generator aircraft with elliptically loaded wings, one can obtain the vortex flow field created by the vortex pair. The vortex pair in this study is modeled by Burnham-Hallock model [44] where tangential velocities induced by the wake are calculated with a constant core radius, The evolution of the wake vortex was simulated by the Wake Vortex Safety System (WVSS) [45] in the out-of-ground effect zone (OGE) as well as in the in-ground effect zone (IGE), where the interaction of the wake vortex and surface boundary layer occurs. The forces and moments acting on the follower aircraft (FA) in both zones are calculated using the velocity field obtained from the simulations using the strip theory and form the disturbance f (x, t) (Equation (1)). In the rest of the paper, the steady level flight regime denotes the wake vortex/aircraft interaction in the OGE zone, and take-off/landing regime describes the interaction in the IGE zone. Furthermore, Monte-Carlo simulations for the steady-level flight regime are performed for the OGE zone. Based on the data from NASA-conducted airport experimental measurements, the real cone of uncertainty of the wake vortex pair evolution (position and strength) was modeled using the WVSS. Thus, the controller performance was investigated within the experiment-based cone of uncertainty. The controller ability to withstand the wake vortex induced disturbances within the terminal zone was demonstrated. The full procedure of the Monte-Carlo approach and validation with the experimental data is discussed in our previous work [31]. Hence, the steady-level flight phase is modeled in the out-of-ground zone.
In the vicinity of the ground, the behavior of the vortex is based on the data from the 3D LES simulations performed in OpenFOAM software and described in detail in our studies focusing on the wake vortex propagation near different types of ground surfaces [32,33]. The velocity fields are extracted from the 3D simulation and are used for the calculation of forces and moments acting on the airplane. When interacting with each other and with the surface boundary layer, the wake vortices generate a turbulent velocity field, which can significantly affect the aircraft. The interaction of the aircraft with the primary and secondary vortices as well as turbulent velocity field is demonstrated.

Estimation of the Maximum Allowable Roll Moment on the Wing
The maximum loading acting on the wing was obtained based on the classical methods of mechanics of materials [46]. Given the cross-section of the wing shown in Figure 4, the section properties were calculated using balsa wood material ( Figure 5). The ultimate stress of the wing was estimated using where M is the internal moment acting on a given section, y is the vertical distance to the neutral axis, and I xx is the area moment of inertia. According to the measurements taken from the wing drawing, the approximate ultimate stress for the wing was calculated based on the data from Table 2 using the linear interpolation for the particular size of the thickest stringer and was equal to 9558 psi. For each simulation, the stress was estimated for every stringer.

Linear, Parameter-Varying Model
For further comparison of linear vs. nonlinear control designs, here, we employ a linear, parameter-varying (LPV) model based on a small fixed-wing UAV. The UAV selected for numerical simulations is the Ultrastick 120 for which the aerodynamics and flight dynamics have been widely studied [47,48]. Figure 6 shows a picture of the UAV airframe. The LPV model is created by linearizing the nonlinear equations of motion about a set of equilibrium points. Here, the LPV model of the Ultrastick 120 holds a constant pitch/yaw angle and varying airspeed.

Highly Turbulent Case and Take-Off/Landing Cases
The velocity fields are extracted from the LES simulation of the wake vortex in the ground vicinity for two cases: highly turbulent and take-off/landing case. The corresponding vorticity fields are shown in Figure 8. The red dashed line displays the path along which the follower aircraft was swept through the wake vortex behind the follower aircraft ( Figure 2). Figures 9-14 show the states' deviations due to wake-vortex-aircraft interaction for a number of points behind the generator aircraft. The secondary vortices are being destroyed due to the interaction with the primary vortices. The destruction is accompanied by the generation of the turbulence with the opposite vorticity, which triggers the breaking up of the vortex pair and the formation of the smaller vortices. A similar process takes place when the vortex pair interacts with the turbulent atmosphere. The larger the turbulence intensity, the more pronounced the decay is. This case demonstrates the behavior and control of the aircraft interacting with highly turbulent upstream wake vorticity. Figures 9 and 10 demonstrate the response of the aircraft to the disturbance created by the highly turbulent velocity fields at H = 54.5 m from the ground surface and expressed in terms of roll rate disturbance and horizontal velocity disturbance. Roll rate p and yaw rate reactions are similar for both controllers in terms of amplitude and convergence rate. The horizontal velocity v oscillations tend to be bigger for the nonlinear controller. The rudder and ailerons' deflections ( Figure 11) are also bigger in the case of the nonlinear controller. In general, the highly turbulent case is characterized by sharp peaks in each state, which correspond to intense oscillations in the velocity field and the resulting forces and moments acting on the aircraft.
The take-off/landing case is shown in Figure 8b and is characterized by fully formed secondary vortices revolving over the primary ones. This case is different from a highly turbulent one by vortex strength. The trajectory of the follower airplane crosses both of them at H = 37 m from the ground. In this case, the nonlinear controller outperforms the linear one: the maximum amplitude of the states' deviation is lower, and the convergence time is shorter (Figures 12 and 13a). However, this advantage is reached at the expense of energy spent on surface deflection ( Figure 14). The wake vortex is characterized by: The convergence rate for roll angle φ, yaw angle ψ, and roll rate (Figures 12 and 13a) is higher in the case of the nonlinear controller. The maximum deviation for all the states is higher for the H ∞ controller, which makes the nonlinear controller more efficient. The aileron deviations shown in Figure 14 are almost the same for both controllers; however, the rudder deflection is bigger in the case of the nonlinear controller, which makes it less energy efficient in this case.

Interaction with Wake Vortex: Steady Level Flight Phase
Several sets of Monte Carlo simulations were implemented to model the steady level flight of the aircraft far from the ground in the OGE zone, as indicated in Figure 3, in the presence of the wake vortex. As long as the real wake vortex is affected by the stochastic nature of the turbulent flow generated by the airplane as well as by the interaction with the atmosphere, the exact position of the vortex pair cannot be predicted by the deterministic model. For this reason, the probabilistic approach [31] was used to generate the cone of uncertainty (for wake vortex position and strength), which predicts not the determined position and strength of the wake but the range within which the wake vortex probability to occur is the highest. The first set of simulations was performed for a nominal trim point. The second set corresponds to simulations with parametric uncertainty. Furthermore, structural analysis to estimate the maximum load on the wing was performed. The maximum allowable roll angle deviation less than 90 degrees was set, and the limit of the roll moment was also taken into account.

Nominal Trim Point Simulations
The Monte Carlo simulations for the wake vortex pair were performed based on the real data case from NASA airport experimental measurements [49] wake vortex data set. The approach is based on perturbing the initial wake vortex conditions (b 0 , V 0 , z 0 ) in the deterministic model and generating profiles of the ambient parameters (such as the EDR and potential temperature profiles) using the probability density functions (PDFs). The PDFs are obtained by applying the maximum likelihood estimation method to density histograms corresponding to a particular data set. The vertical profiles are truncated at the heights of the vortex generation and the mean values are calculated from the zero height to that value. A total of 400 perturbations for vertical profiles were generated using PDFs and used as input parameters along with other parametric perturbations. More details of the method can be found in Ref. [34]. The characteristics of the wake generator aircraft The distance to the generator's aircraft path line as well as 400 wake vortex realizations were used as an input to the linear and nonlinear control system. The vertical position of the follower aircraft was fixed at z = 175 m. The separation distance between the vortices in the vortex pair, the generator aircraft's weight, and the vertical distance to the generator aircraft's flight path are the input parameters for the controllers' Monte Carlo simulations. Figure 15 shows the nondimensionalized circulation and vertical coordinate of the vortex pair obtained in the first set of simulations, which forms the cone of uncertainty in the real wake vortex evolution. The response of the follower aircraft at the distances of 5.46 nm, 4.16 nm, and 2.34 nm, which correspond tot = 5,t = 3.84, andt = 2.14, is considered. Figure 15 demonstrates the cone of the uncertainty of the wake vortex in terms of vertical coordinate and its strength. Blue dashed lines show the cut that corresponds to the moment of time considered. The plots of the states for each Monte Carlo simulation (t = 5) are shown in Figure 16. The average Root Mean Square Error (RMSE) and Average Maximum Deviations (AMD) of deflections surfaces metrics are chosen to demonstrate the performance of the controller. The values of the average RMSE and AMD for each state are summarized in Tables 3-5. The RMSE for all the states shows that the nonlinear controller is more effective than the linear one in terms of states' deviations. However, the nonlinear controller uses the rudder more actively and thus expends more energy for control. Furthermore, from the comparison of the RMSE values, one can see that the nonlinear controller works better than the linear one.      The advantage of the nonlinear robust controller is the ability to suppress the disturbance at different flight conditions and aircraft parameter's deviations. The LPV simulation is performed for several trim points (equilibrium points), which corresponds to the change of parameters of the aircraft and the flight conditions. At the same time, the linear and nonlinear controllers are designed for the nominal trim point. The results are presented in Table 6. The states' average RMSE and AMD are summarized fort = 2.14, equilibrium point #7, which corresponds to 27 m/s flight with 3.5 deg angle of attack (Figure 7). It is clear that the nonlinear controller outperforms the linear one. However, this advantage is due to the use of a higher amount of energy spent on the rudder deflection. This subsection demonstrates the performance comparison between the nonlinear robust controller and the proportional-integral-derivative controller (PID) in the presence of parametric variations. We model these conditions by choosing different equilibrium (trim) points (Figure 7), as was conducted in the previous section. Both controllers are tuned based on the nominal trim point conditions. In order to compare the performance of two controllers, the Ultrastick 120 is set to follow the B-733 aircraft at the distance where the maximum circulation of the wake is 60% of Γ 0 . The parameters of generator aircraft are: • Circulation Γ 0 = 230 m 2 /s; • Span b = 28.9 m; • Velocity V = 65.7 m/s. Figures 17 and 18 show the performance of the PID controller vs. the robust nonlinear controller for the nominal trim point. In this case, both controllers show the ability to suppress the disturbance. However, at trim point #1, the PID controller is not able to compensate for the disturbance, as shown in Figures 19 and 20. Indeed, the nonlinear robust controller's performance is at the same level as it is at the nominal trim point. This demonstrates the effectiveness of the robust nonlinear control law in the presence of the disturbance.

Conclusions
A robust nonlinear control method is developed and is demonstrated to achieve roll/yaw regulation in the presence of unmodeled external disturbances (i.e., wake vortices and wind gusts) and system nonlinearities. The results of numerical simulations were provided to demonstrate the performance of the proposed nonlinear control law. For comparison, the same control objective was simulated using a linear H ∞ control law as well as a PID controller. It was shown that the nonlinear control method compensates for the disturbances similarly or more effectively than the linear controller. The PID controller is not capable of rejecting the disturbances in the presence of system nonlinearities. The results showed that the nonlinear control design outperformed the linear control method for the simulated trajectory regulation objective under the tested levels of uncertainty.
The performance of the nonlinear controller was successfully tested in the presence of the wake vortex disturbance in the vicinity of the ground as well as far from the ground, which corresponds to take-off/landing and steady level flight conditions. The interaction of wake vortex-induced turbulence field and aircraft was modeled using high-fidelity Large Eddy Simulations (LES) simulations, which reflect real flight conditions in the terminal zones. The Monte Carlo simulations for the out-of-ground interaction were performed based on the data from the NASA experimental measurements. Using this data, the cone of uncertainty for the wake-vortex evolution was modeled and employed in the simulation. The performance of the robust nonlinear controller and H ∞ linear controller were compared for both in-ground and out-of-ground zones. The nonlinear controller outperformed the linear one in terms of average Root Mean Square Error (RMSE) and Average Maximum Deviation (AMD).
The results demonstrate the capability of the proposed nonlinear controller to asymptotically reject wake vortex disturbance in the presence of a high degree of uncertainty and nonlinearity in the UAV operating conditions. Moreover, the nonlinear controller is designed with a computationally efficient structure, which is motivated by UAV applications where onboard computational resources are limited.
As a future step, the authors are planning to test the performance of proposed controller on the Ultrastick-120 model in a wind tunnel for the set of wake vortex strengths.

Conflicts of Interest:
The authors declare no conflict of interest.