Robust Control of UAV with Disturbances and Uncertainty Estimation

Abstract

In this work, a nonlinear estimator-based robust controller is designed for the position and yaw control of a quadrotor with uncertainty estimation. This controller ensures the tracking of desired references in the presence of parameters variation and external disturbances, making use of high-order sliding mode (HOSM) estimators to estimate these perturbations that can be canceled by the control, thus improving the dynamic behavior of the controlled system. Its performance is evaluated making use of a Simcenter Amesim quadrotor based on physical models generated from experimental data in a co-simulation framework with Matlab–Simulink used to implement the designed controller with FPGA implementation. A challenging and generic maneuver with time-varying wind disturbances and uncertainty model parameters is considered.

1. Introduction

In recent years UAVs have attracted attention due to the vast range of applications that can be addressed, and they have been widely utilized by researchers, law enforcement and security, delivery companies, farmers, etc. Furthermore, UAVs are evolving from being only important for the sensing action [1] to being fundamental also for the actuator action. In the upcoming years, there will therefore be a remarkable development of drones for functional applications such as payload delivery services to transport deformable linear objects such as hoses and goods to carry organs for organ transplants, which helped reduce car journeys to just a few minutes of flight. Furthermore, a recent study showed that the delivery of UAVs could even help reduce greenhouse gas emissions caused by the freight industry [2]. While the application interest of UAVs is growing, on the other hand, there is a need to make the system more and more efficient. The main limitation that occurs in all applications, however, is the limited amount of energy that batteries or other energy sources provide. In fact, in the drones on the market, the average flight time is between 15 and 45 min, which is a time that highly limits the real application. Improvement and optimization of the source energy, lightening of the drone, optimization of the path and improvement of the performance in the presence of disturbances that modify its behavior and significantly decrease performance from an energy point of view are some of the actions investigated to solve the problem. The execution of trajectories that allow for lower consumption is important. In fact, one of the factors that often influences the consumption is linked to disturbances, such as wind, which leads to a departure from the ideal trajectory, and it is then necessary to compensate through additional time and higher battery consumption. In the same way, the model on which the control is based is an approximation, the parameters are time-varying, and it is therefore necessary to intervene also on the parametric uncertainties. It is a challenge that must be faced from different points of view and that requires a multidisciplinary effort in order to achieve the goal. Motivated by the European project COMP4DRONES (see [3]), it was necessary to study methodologies for controlling the flight of a drone in precision agriculture, to monitor the health status of plants, and possibly spray fertilizer to the plants with time limits, trying to minimize the energy and to follow desired trajectory in presence of disturbances and parametric uncertainties. This represents a complex problem with different objectives. In [4], a robust flight controller based on linear active disturbance rejection control is proposed for the stability control of an aerial robot quadrotor under wind gusts. The nonlinear dynamical model of the quadrotor, considering the wind disturbance, is firstly established through Newton–Euler method. In this control scheme, the linear extended state observer serves as a compensator which can effectively reject the wind gusts. In this case, the disturbance due to the wind is therefore compensated, but a perfect knowledge of the model and its parameters is assumed. In [5], a robust dynamic sliding mode control algorithm using a nonlinear disturbance observer for system dynamics is investigated. The proposed method is applied to provide a rapid adaptation and strictly robust performance for the attitude and altitude control of UAVs. The procedure of the proposed method consists of two stages. First, a nonlinear disturbance observer is applied to estimate the exogenous perturbation. Second, a robust dynamic sliding mode controller integrated with the estimated values of disturbances is presented by a combination of a PID sliding surface and super twisting technique to compensate for the effect of these perturbations on the system. Again, external noise is considered but not model-related uncertainties. Furthermore, the control is tested on an ideal plant. In [6], a novel extension of the classic Nonlinear Dynamic Inversion (NLDI) control architecture for wind disturbance rejection is designed by using adaptive artificial neural networks. In [7], an adaptive controller for a quadrotor UAV for carrying unknown payloads while tracking any trajectory is proposed. The proposed adaptive controller is robust to modeling uncertainties and does not require any a priori knowledge of the bounds of the uncertainties. The controller is also robust to time-varying delays without any constraint on the derivative of the time delay. The simulations are conducted on the simple basic model and do not take into account possible disturbances, such as wind, which influence performance in a non-negligible way. In [8], the deep reinforcement learning method, and robust deep deterministic policy gradient, are proposed for developing a controller that allows robust flying of an unmanned aerial vehicle in dynamic uncertain environments. While a very interesting work, it is tested in a simplified scenario, and the wind and the parametric uncertainties of the model are not directly considered. Furthermore, there is no formal proof of the stability of the controlled system.

The contribution of the paper is dual: on the one hand, it designs a controller that is robust both to the uncertainties of the model and to external disturbances such as wind or aerodynamic factors using HOSM estimators; on the other side, this controller is tested in a co-simulation environment between Matlab–Simulink and Simcenter Amesim on FPGA. In fact, to our best knowledge, results concerning the quadrotor tracking problem in the case of variations of all the parameters emerging in the dynamics and in the presence of environmental disturbances and/or unmodeled dynamics influencing dynamics of the UAV are not available in the research literature. Moreover, after preliminary tests on Matlab–Simulink implementation, that have not been included in this paper, the proposed controller is tested on the quadrotor model provided by the Simcenter Amesim developed by Siemens, based on physical models, in a co-simulation environment with FPGA implementation. As a matter of fact, Simcenter Amesim is proven to accurately represent the dynamic behavior of the drone, it was extensively validated with experimental tests, and it is used by important aerospace enterprises worldwide. In fact, the Simcenter Amesim drone is characterized by components with physical models, unmodeled dynamics, time-varying parameters, disturbances, and parametric uncertainties that make the problem complicated and validate the research conducted meaningfully. The work is organized as follows. In Section 2, the quadrotor mathematical model is presented. The estimators of disturbances and a nonlinear controller to ensure robustness and asymptotic convergence are designed in Section 3. In Section 4, the performance of the controller and estimators are presented with respect to Matlab and Simcenter Amesim co-simulation environment implementation with FPGA implementation. Some considerations and future research developments conclude the work in Section 5.

2. Mathematical Model of A Quadrotor

The quadrotor considered in this work consists of a rigid frame equipped with four rotors (see [9,10] for a detailed description). The rotors generate the propeller force $F_i=b·{\omega }_{p,i}^2$, which is proportional to the propeller angular velocity ${\omega }_{p,i}$, $i=1,2,3,4$. Propellers 1 and 3 rotate counterclockwise; propellers 2 and 4 rotate clockwise. Quadrotor orientation with Euler angles is used. Let us indicate with $RC(O,e_1,e_2,e_3)$, $R\mathrm{\Gamma }(\mathrm{\Omega },{\epsilon }_1,{\epsilon }_2,{\epsilon }_3)$) the frames fixed with the Earth and the quadrotor, respectively, and $\mathrm{\Omega }$ is coincident with the center of mass of the quadrotor (see Figure 1).

The absolute position of the quadrotor in $RC$ is described by $p={(x,y,z)}^T$, whereas its attitude is described by the Euler angles $\alpha ={(\varphi ,\theta ,\psi )}^T$, with $\varphi \in [-\pi /2,\pi /2)$, $\theta \in (-\pi /2,\pi /2)$, $\psi \in [-\pi ,\pi )$ the roll, pitch, yaw angles, respectively. The sequence 3–2–1 is here considered [11]. Moreover, $v={(v_x,v_y,v_z)}^T$, $\omega ={({\omega }_1,{\omega }_2,{\omega }_3)}^T$ are the linear and angular velocities of the center of mass of the quadrotor, expressed in $RC$ and in $R\mathrm{\Gamma }$, respectively. The translation dynamics, expressed in $RC$, and rotation dynamics, expressed in $R\mathrm{\Gamma }$, of the quadrotor are

$$\begin{array}{cc}\hfill \dot{p}& =v\hfill \\ \hfill \dot{v}& =\frac{1}{m}\left(f_p+f_g+f_d\right)\hfill \\ \hfill \dot{\alpha }& =M\left(\alpha \right)\omega \hfill \\ \hfill \dot{\omega }& =J^{-1}\left(-\tilde{\omega }J\omega +{\tau }_p+{\tau }_{gy}+{\tau }_d\right)\hfill \end{array}$$

(1)

where m is the mass of the quadrotor, $f_p=\mathcal{R}\left(\alpha \right){\phi }_p$ is the force in $RC$ excerted by the propellers (${\phi }_p$ is the same force, but expressed in $R\mathrm{\Gamma }$), J (positive definite symmetric matrix in ${\mathbb{R}}^{3\times 3}$, expressed in $R\mathrm{\Gamma }$) is the inertia matrix of the quadrotor, and

$$\tilde{\omega }=\left(\begin{array}{ccc}0& -{\omega }_3& {\omega }_2\\ {\omega }_3& 0& -{\omega }_1\\ -{\omega }_2& {\omega }_1& 0\end{array}\right)$$

(2)

is the so–called dyadic representation of $\omega$. Furthermore,

$${\phi }_p=\left(\begin{array}{c}0\\ 0\\ u_p\end{array}\right),u_p=\sum _{i=1}^4{F}_i,f_g=\left(\begin{array}{c}0\\ 0\\ -mg\end{array}\right)$$

(3)

$$\begin{array}{c}\hfill {\tau }_p=\left(\begin{array}{c}\hfill {\tau }_1\\ \hfill {\tau }_2\\ \hfill {\tau }_3\end{array}\right)=\left(\begin{array}{c}\ell b({\omega }_{p,2}^2-{\omega }_{p,4}^2)\\ \ell b({\omega }_{p,3}^2-{\omega }_{p,1}^2)\\ c({\omega }_{p,1}^2-{\omega }_{p,2}^2+{\omega }_{p,3}^2-{\omega }_{p,4}^2)\end{array}\right)\end{array}$$

(4)

are the input forces and moments produced by the propellers (expressed in $R\mathrm{\Gamma }$), where ℓ is the distance between the center of mass $CG$ to the rotor shaft, and b ([b]=$\frac{N{s}^2}{ra{d}^2}$) and c ([c]=$\frac{N{s}^{2}m}{ra{d}^2}$) are, respectively, the so–called thrust and drag factor. Obviously,

$$\begin{array}{c}\left(\begin{array}{c}{F}_1\\ F_2\\ F_3\\ F_4\end{array}\right)={\left(\begin{array}{cccc}1& 1& 1& 1\\ 0& \ell & 0& -\ell \\ -\ell & 0& \ell & 0\\ c& -c& c& -c\end{array}\right)}^{-1}\left(\begin{array}{c}{u}_p\\ {\tau }_1\\ {\tau }_2\\ {\tau }_3\end{array}\right),\left(\begin{array}{c}{\omega }_{p,1}^2\\ {\omega }_{p,2}^2\\ {\omega }_{p,3}^2\\ {\omega }_{p,4}^2\end{array}\right)=\left(\begin{array}{cccc}\frac{1}{4b}& 0& \frac{-1}{2b\ell }& \frac{1}{4c}\\ \frac{1}{4b}& \frac{1}{2b\ell }& 0& \frac{-1}{4c}\\ \frac{1}{4b}& 0& \frac{1}{2b\ell }& \frac{1}{4c}\\ \frac{1}{4b}& \frac{-1}{2b\ell }& 0& \frac{-1}{4c}\end{array}\right)\left(\begin{array}{c}{u}_p\\ {\tau }_1\\ {\tau }_2\\ {\tau }_3\end{array}\right),\end{array}$$

(5)

$${\omega }_{p,i}=\sqrt{\frac{F_i}{b}},0\le F_i\le F_{i,max},0\le {\omega }_{p,i}\le {\omega }_{p,i,max},i=1,2,3,4$$

where $F_{i,max}$ and ${\omega }_{p,i,max}$ are, respectively, the maximum forces and the angular velocities related to each propellers due to physical limitations. Moreover, $f_g$ in (3) is the force due to the gravity, which is expressed in $RC$. The vectors expressed in $R\mathrm{\Gamma }$ are transformed into vectors in $RC$ by the rotation matrix

$$\begin{array}{c}{\mathcal{R}}_{ib}=\mathcal{R}\left(\alpha \right)=\left(\begin{array}{ccccc}{c}_{\theta }{c}_{\psi }& & s_{\varphi }{s}_{\theta }{c}_{\psi }-c_{\varphi }{s}_{\psi }& & c_{\varphi }{s}_{\theta }{c}_{\psi }+s_{\varphi }{s}_{\psi }\\ c_{\theta }{s}_{\psi }& & s_{\varphi }{s}_{\theta }{s}_{\psi }+c_{\varphi }{c}_{\psi }& & c_{\varphi }{s}_{\theta }{s}_{\psi }-s_{\varphi }{c}_{\psi }\\ -s_{\theta }& & s_{\varphi }{c}_{\theta }& & c_{\varphi }{c}_{\theta }\end{array}\right)\end{array}$$

(6)

where $c_{\star }=cos(\star )$, $s_{\star }=sin(\star )$, $\star =\varphi ,\theta ,\psi$. The angular velocity dynamics are expressed using the matrix

$$M\left(\alpha \right)=\left(\begin{array}{ccccc}1& & s_{\varphi }{\mathrm{tg}}_{\theta }& & c_{\varphi }{\mathrm{tg}}_{\theta }\\ 0& & c_{\varphi }& & -s_{\varphi }\\ 0& & s_{\varphi }{\mathrm{sc}}_{\theta }& & c_{\varphi }{\mathrm{sc}}_{\theta }\end{array}\right)$$

with ${\mathrm{tg}}_{\star }=tan(\star )$, ${\mathrm{sc}}_{\star }=sec(\star )$, $\star =\varphi ,\theta ,\psi$. Under the assumption of small angles $\varphi$ (roll) and $\theta$ (pitch), which is acceptable in the case of a quadrotor performing non-aggressive maneuvers, this matrix can be approximated by the identity matrix, i.e., $M\left(\alpha \right)\simeq I_{3x3}$ [12]. The rolling torque ${\tau }_1$ is produced by the forces $F_2$ and $F_4$. Similarly, the pitching torque ${\tau }_2$ is produced by the forces $F_1$ and $F_3$. Due to the Newton’s third law, the propellers produce a yawing torque ${\tau }_3$ on the body of the quadrotor in the opposite direction of the propeller rotation. Moreover, ${\tau }_{gy}={\sum }_{i=1}^4{(-1)}^i{J}_{p,i}{\omega }_{p,i}\tilde{\omega }{ϵ}_3$ is the gyroscopic torque due to the propeller rotations, with $J_{p,i}$, $i=1,2,3,4$, the inertia moment of the $i_{\mathrm{th}}$ motor plus propeller, with respect to its rotation axis. Finally, $f_d$ and ${\tau }_d$ are the forces and torques due to the external disturbances, for example due to aerodynamic forces. Considering the nominal values $m^{\circ }$, $J_{p,i}^{\circ }=J_p^{\circ }$, $i=1,2,3,4$, $J^{\circ }=\mathrm{diag}\{J_1^{\circ },J_2^{\circ },J_3^{\circ }\}$, $f_p^{\circ }$, ${\tau }_p^{\circ }$, system (1) can be rewritten as

$$\begin{array}{cc}\hfill \dot{p}& =v\hfill \\ \hfill \dot{v}& =\frac{1}{m^{\circ }}{f}_p^{\circ }+\frac{1}{m}{f}_g+d\hfill & \\ \hfill \dot{\alpha }& =\omega \hfill \\ \hfill \dot{\omega }& =-J^{\circ }{}^{-1}\tilde{\omega }{J}^{\circ }\omega +J^{\circ }{}^{-1}{\tau }_p^{\circ }+J^{\circ }{}^{-1}{\tau }_{gy}^{\circ }+\delta \hfill \end{array}$$

(7)

where

$$\begin{array}{cc}\hfill d=\left(\begin{array}{c}{d}_x\\ d_y\\ d_z\end{array}\right)& =\left(\frac{1}{m}-\frac{1}{m^{\circ }}\right)f_p+\frac{1}{m}{f}_d+\frac{1}{m^{\circ }}\left(f_p-f_p^{\circ }\right)\hfill \\ \delta =\left(\begin{array}{c}{\delta }_1\\ {\delta }_2\\ {\delta }_3\end{array}\right)& =-\left(J^{-1}\tilde{\omega }J\omega -J^{\circ }{}^{-1}\tilde{\omega }{J}^{\circ }\omega \right)+\left(J^{-1}-J^{\circ }{}^{-1}\right){\tau }_p+\left(J^{-1}{\tau }_{gy}-J^{\circ }{}^{-1}{\tau }_{gy}^{\circ }\right)+J^{-1}{\tau }_d+\\ \hfill & +J^{\circ }{}^{-1}({\tau }_p-{\tau }_p^{\circ })\hfill \end{array}$$

(8)

are the perturbations acting on the drone, due to parametric uncertainties and external disturbances, and

$$\begin{array}{c}\hfill {\tau }_{gy}^{\circ }=\sum _{i=1}^4{(-1)}^i{J}_p^{\circ }{\omega }_{p,i}\tilde{\omega }{ϵ}_3=J_p^{\circ }{\omega }_p\left(\begin{array}{c}{\omega }_2\\ -{\omega }_1\\ 0\end{array}\right)\end{array}$$

(9)

is the nominal gyroscopic torque due to the propeller rotations, with ${\omega }_p=-{\omega }_{p,1}+{\omega }_{p,2}-{\omega }_{p,3}+{\omega }_{p,4}$ the so-called rotor relative speed. With these positions, mathematical model (7) of the quadrotor can be expressed as

$$\begin{array}{cc}\hfill \dot{x}& =v_x,\dot{y}=v_y,\dot{z}=v_z\hfill \\ \hfill {\dot{v}}_x& =\left(c_{\varphi }{s}_{\theta }{c}_{\psi }+s_{\varphi }{s}_{\psi }\right)\frac{u_p^{\circ }}{m^{\circ }}+d_x\hfill \\ \hfill {\dot{v}}_y& =\left(c_{\varphi }{s}_{\theta }{s}_{\psi }-s_{\varphi }{c}_{\psi }\right)\frac{u_p^{\circ }}{m^{\circ }}+d_y\hfill \\ \hfill {\dot{v}}_z& =-g+c_{\varphi }{c}_{\theta }\frac{u_p^{\circ }}{m^{\circ }}+d_z\hfill \\ \hfill \dot{\varphi }& ={\omega }_1,\dot{\theta }={\omega }_2,\dot{\psi }={\omega }_3\hfill \\ \hfill {\dot{\omega }}_1& =\frac{J_2^{\circ }-J_3^{\circ }}{J_1^{\circ }}{\omega }_2{\omega }_3+\frac{J_p^{\circ }}{J_1^{\circ }}{\omega }_p{\omega }_2+\frac{1}{J_1^{\circ }}{\tau }_1^{\circ }+{\delta }_1\hfill \\ \hfill {\dot{\omega }}_2& =\frac{J_3^{\circ }-J_1^{\circ }}{J_2^{\circ }}{\omega }_1{\omega }_3-\frac{J_p^{\circ }}{J_2^{\circ }}{\omega }_p{\omega }_1+\frac{1}{J_2^{\circ }}{\tau }_2^{\circ }+{\delta }_2\hfill \\ \hfill {\dot{\omega }}_3& =\frac{J_1^{\circ }-J_2^{\circ }}{J_3^{\circ }}{\omega }_1{\omega }_2+\frac{1}{J_3^{\circ }}{\tau }_3^{\circ }+{\delta }_3.\hfill \end{array}$$

(10)

Remark 1. 

Within the vectors d and δ, there are the variations of mass m, moments of inertia, thrust coefficient b and drag coefficient c linked to the propulsion force $F_p$, as well as possible external perturbations and unmodeled dynamics that are considered in $F_d$ and ${\tau }_d$, so as to incorporate the parametric variations, the unmodeled dynamics and external perturbations (as wind) that act both at the level of linear and angular velocity level (see Section 4.2).

Assumption A1. 

$f_p$ and ${\tau }_d$, and their derivatives, are bounded.

3. Design of Nonlinear Controller

The output vector to be controlled is composed by the 3D position and yaw rate and the control aim is to design a controller, in the presence of disturbances and parameters uncertainty, to ensure the asymptotic converge of the variables $\chi =(x,y,z,\psi )$ to some reference trajectories ${\chi }_r=(x_r,y_r,z_r,{\psi }_r)$. To this aim, let us consider the kinematic errors

$$e_{\chi }=\chi -{\chi }_r=\left(\begin{array}{c}x-x_r\\ y-y_r\\ z-z_r\\ \psi -{\psi }_r\end{array}\right)=\left(\begin{array}{c}{e}_x\\ e_y\\ e_z\\ e_{\psi }\end{array}\right)$$

whose dynamics are

$${\dot{e}}_{\chi }=\left(\begin{array}{c}{\dot{e}}_x\\ {\dot{e}}_y\\ {\dot{e}}_z\\ {\dot{e}}_{\psi }\end{array}\right)=\left(\begin{array}{c}{v}_x-{\dot{x}}_r\\ v_y-{\dot{y}}_r\\ v_z-{\dot{z}}_r\\ {\omega }_3-{\dot{\psi }}_r\end{array}\right)=\left(\begin{array}{c}{e}_{v_x}\\ e_{v_y}\\ e_{v_z}\\ e_{{\omega }_3}\end{array}\right)$$

and with

$$\begin{array}{cc}\hfill {\ddot{e}}_{\chi }& =\left(\begin{array}{c}\hfill {\ddot{e}}_x\\ \hfill {\ddot{e}}_y\\ \hfill {\ddot{e}}_z\\ \hfill {\ddot{e}}_{\psi }\end{array}\right)=\left(\begin{array}{c}\hfill {\dot{v}}_x-{\ddot{x}}_r\\ \hfill {\dot{v}}_y-{\ddot{y}}_r\\ \hfill {\dot{v}}_z-{\ddot{z}}_r\\ \hfill {\dot{\omega }}_3-{\ddot{\psi }}_r\end{array}\right)=\left(\begin{array}{c}\hfill {\dot{e}}_{v_x}\\ \hfill {\dot{e}}_{v_y}\\ \hfill {\dot{e}}_{v_z}\\ \hfill {\dot{e}}_{{\omega }_1}\end{array}\right)=\left(\begin{array}{c}\left(c_{\varphi }{s}_{\theta }{c}_{\psi }+s_{\varphi }{s}_{\psi }\right){\displaystyle \frac{u_p^{\circ }}{m^{\circ }}}+d_x-{\ddot{x}}_r\\ \left(c_{\varphi }{s}_{\theta }{s}_{\psi }-s_{\varphi }{c}_{\psi }\right){\displaystyle \frac{u_p^{\circ }}{m^{\circ }}}+d_y-{\ddot{y}}_r\\ -g+c_{\varphi }{c}_{\theta }{\displaystyle \frac{u_p^{\circ }}{m^{\circ }}}+d_z-{\ddot{z}}_r\\ {\displaystyle \frac{J_1^{\circ }-J_2^{\circ }}{J_3^{\circ }}}{\omega }_1{\omega }_2+{\displaystyle \frac{1}{J_3^{\circ }}}{\tau }_3^{\circ }+{\delta }_3-{\ddot{\psi }}_r\end{array}\right)=\hfill \\ & =\left(\begin{array}{c}\hfill \left(c_{\varphi }{s}_{{\theta }_r}{c}_{\psi }+s_{{\varphi }_r}{s}_{\psi }\right){\displaystyle \frac{u_p^{\circ }}{m^{\circ }}}+d_x-{\ddot{x}}_r\\ \left(c_{\varphi }{s}_{{\theta }_r}{s}_{\psi }-s_{{\varphi }_r}{c}_{\psi }\right){\displaystyle \frac{u_p^{\circ }}{m^{\circ }}}+d_y-{\ddot{y}}_r\\ -g+c_{\varphi }{c}_{\theta }{\displaystyle \frac{u_p^{\circ }}{m^{\circ }}}+d_z-{\ddot{z}}_r\\ {\displaystyle \frac{J_1^{\circ }-J_2^{\circ }}{J_3^{\circ }}}{\omega }_1{\omega }_2+{\displaystyle \frac{1}{J_3^{\circ }}}{\tau }_3^{\circ }+{\delta }_3-{\ddot{\psi }}_r\end{array}\right)+\left(\begin{array}{c}\left((s_{\theta }-s_{{\theta }_r})c_{\varphi }{c}_{\psi }+(s_{\varphi }-s_{{\varphi }_r})s_{\psi }\right){\displaystyle \frac{u_p^{\circ }}{m^{\circ }}}\\ \left((s_{\theta }-s_{{\theta }_r})c_{\varphi }{s}_{\psi }-(s_{\varphi }-s_{{\varphi }_r})c_{\psi }\right){\displaystyle \frac{u_p^{\circ }}{m^{\circ }}}\\ 0\\ 0\end{array}\right).\end{array}$$

Setting

$$\begin{array}{cc}\hfill \left(\begin{array}{c}\hfill s_{{\varphi }_r}\\ \hfill s_{{\theta }_r}\end{array}\right)& =\frac{1}{\sigma }\left(\begin{array}{cc}{c}_{\varphi }{c}_{\theta }{s}_{\psi }& -c_{\varphi }{c}_{\theta }{c}_{\psi }\\ c_{\theta }{c}_{\psi }& c_{\theta }{s}_{\psi }\end{array}\right)\left(\begin{array}{c}{\ddot{x}}_r-{\widehat{d}}_x+w_x\\ {\ddot{y}}_r-{\widehat{d}}_y+w_y\end{array}\right)\hfill \\ \hfill u_p^{\circ }& =\frac{m^{\circ }}{c_{\varphi }{c}_{\theta }}\left({\ddot{z}}_r+g-{\widehat{d}}_z+w_z\right)\hfill \\ \hfill {\tau }_3^{\circ }& =J_3^{\circ }\left({\ddot{\psi }}_r-{\displaystyle \frac{J_1^{\circ }-J_2^{\circ }}{J_3^{\circ }}}{\omega }_1{\omega }_2-{\widehat{\delta }}_3+w_{\psi }\right)\hfill \\ \hfill w_j& =-k_{j,0}{I}_{e_j}-k_{j,1}{e}_j-k_{j,2}{e}_{v_j},j=x,y,z\hfill \\ \hfill w_{\psi }& =-k_{\psi ,0}{I}_{e_{\psi }}-k_{\psi ,1}{e}_{\psi }-k_{\psi ,2}{e}_{{\omega }_3}\hfill \\ \hfill {\dot{I}}_{e_j}& =e_j,j=x,y,z,\psi \hfill \end{array}$$

(11)

where $\sigma ={\ddot{z}}_r+g-{\widehat{d}}_z+w_z$, $k_{j,0},k_{j,2}>0$ and $k_{j,1}>k_{j,0}/k_{j,2}$, $j=x,y,z,\psi$, one gets

$$\begin{array}{cc}\hfill {\ddot{e}}_{\chi }& =\left(\begin{array}{c}{w}_x\\ w_y\\ \\ w_z\\ w_{\psi }\end{array}\right)+\left(\begin{array}{c}{\mu }_x{\displaystyle \frac{u_p^{\circ }}{m^{\circ }}}\\ {\mu }_y{\displaystyle \frac{u_p^{\circ }}{m^{\circ }}}\\ 0\\ 0\end{array}\right)+\left(\begin{array}{c}{d}_x-{\widehat{d}}_x\\ d_y-{\widehat{d}}_y\\ \\ d_z-{\widehat{d}}_z\\ {\delta }_3-{\widehat{\delta }}_3\end{array}\right),\hfill \end{array}$$

where

$${\mu }_x=[(s_{\theta }-s_{{\theta }_r})c_{\varphi }{c}_{\psi }+(s_{\varphi }-s_{{\varphi }_r})s_{\psi }],$$

$${\mu }_y=[(s_{\theta }-s_{{\theta }_r})c_{\varphi }{s}_{\psi }-(s_{\varphi }-s_{{\varphi }_r})c_{\psi }].$$

As far as the errors $s_{\varphi }-s_{{\varphi }_r}$, $s_{\theta }-s_{{\theta }_r}$ are concerned, let us consider the roll error $e_{\varphi }=\varphi -{\varphi }_r$ and the pitch error $e_{\theta }=\theta -{\theta }_r$, with ${\varphi }_r=arcsin{s}_{{\varphi }_r}$, ${\theta }_r=arcsin{s}_{{\theta }_r}$. Their dynamics are

$$\begin{array}{cc}\hfill {\dot{e}}_{\varphi }& ={\omega }_1-{\dot{\varphi }}_r=e_{{\omega }_1},{\dot{e}}_{\theta }={\omega }_2-{\dot{\theta }}_r=e_{{\omega }_2}\hfill \\ \hfill {\dot{e}}_{{\omega }_1}& =\frac{J_2^{\circ }-J_3^{\circ }}{J_1^{\circ }}{\omega }_2{\omega }_3+\frac{J_p^{\circ }}{J_1^{\circ }}{\omega }_p{\omega }_2+\frac{1}{J_1^{\circ }}{\tau }_1^{\circ }+{\delta }_1-{\ddot{\varphi }}_r\hfill \\ \hfill {\dot{e}}_{{\omega }_2}& =\frac{J_3^{\circ }-J_1^{\circ }}{J_2^{\circ }}{\omega }_1{\omega }_3-\frac{J_p^{\circ }}{J_2^{\circ }}{\omega }_p{\omega }_1+\frac{1}{J_2^{\circ }}{\tau }_2^{\circ }+{\delta }_2-{\ddot{\theta }}_r\hfill \end{array}$$

and the controls

$$\begin{array}{cc}\hfill {\tau }_1^{\circ }& =J_1^{\circ }{\ddot{\varphi }}_r-(J_2^{\circ }-J_3^{\circ }){\omega }_2{\omega }_3-J_p^{\circ }{\omega }_p{\omega }_2-J_1^{\circ }{\widehat{\delta }}_1+J_1^{\circ }{w}_{\varphi }\hfill \\ \hfill {\tau }_2^{\circ }& =J_2^{\circ }{\ddot{\theta }}_r-(J_3^{\circ }-J_1^{\circ }){\omega }_1{\omega }_3+J_p^{\circ }{\omega }_p{\omega }_1-J_2^{\circ }{\widehat{\delta }}_2+J_2^{\circ }{w}_{\theta }\hfill \\ \hfill w_{\varphi }& =-k_{\varphi ,0}{I}_{e_{\varphi }}-k_{\varphi ,1}{e}_{\varphi }-k_{\varphi ,2}{e}_{{\omega }_1}\hfill \\ \hfill w_{\theta }& =-k_{\theta ,0}{I}_{e_{\theta }}-k_{\theta ,1}{e}_{\theta }-k_{\theta ,2}{e}_{{\omega }_2}\hfill \end{array}$$

(12)

with $k_{j,0},k_{j,2}>0$, $k_{j,1}>k_{j,2}/k_{j,0}$, $j=\varphi ,\theta$. Therefore, one obtains

$$\begin{array}{cc}\hfill \left(\begin{array}{c}{\ddot{e}}_x\\ {\ddot{e}}_y\\ {\ddot{e}}_z\end{array}\right)& =\left(\begin{array}{c}{w}_x\\ w_y\\ w_z\end{array}\right)+\left(\begin{array}{c}{\mu }_x{\displaystyle \frac{u_p^{\circ }}{m^{\circ }}}\\ {\mu }_y{\displaystyle \frac{u_p^{\circ }}{m^{\circ }}}\\ 0\end{array}\right)+\left(\begin{array}{c}{d}_x-{\widehat{d}}_x\\ d_y-{\widehat{d}}_y\\ d_z-{\widehat{d}}_z\end{array}\right)\hfill \end{array}$$

(13a)

$$\begin{array}{cc}\hfill \left(\begin{array}{c}{\ddot{e}}_{\varphi }\\ {\ddot{e}}_{\theta }\\ {\ddot{e}}_{\psi }\end{array}\right)& =\left(\begin{array}{c}{w}_{\varphi }\\ w_{\theta }\\ w_{\psi }\end{array}\right)+\left(\begin{array}{c}{\delta }_1-{\widehat{\delta }}_1\\ {\delta }_2-{\widehat{\delta }}_2\\ {\delta }_3-{\widehat{\delta }}_3\end{array}\right)\end{array}$$

(13b)

Remark 2. 

For the computation of the first and second derivatives of the signals ϕ and ψ, since it is not possible to have a direct analytical expression, they will be implemented using the high-order sliding mode differentiator (HOSMD) designed in [13].

Disturbances Estimator

To achieve robustness of the closed-loop system with respect to external disturbances, a robust disturbance estimator is designed. In what follows, given the measured variables the HOSM estimator wants to converge in finite time to disturbances values, finally, based on [14], one considers,

$$\begin{array}{cc}\hfill {\dot{\nu }}_x& =\left(c_{\varphi }{s}_{\theta }{c}_{\psi }+s_{\varphi }{s}_{\psi }\right){\displaystyle \frac{u_p^{\circ }}{m^{\circ }}}-{\lambda }_{x,1}{⌊{\nu }_x-v_x⌉}^{1/2}+{\widehat{d}}_x\hfill \\ \hfill {\dot{\nu }}_y& =\left(c_{\varphi }{s}_{\theta }{s}_{\psi }-s_{\varphi }{c}_{\psi }\right){\displaystyle \frac{u_p^{\circ }}{m^{\circ }}}-{\lambda }_{y,1}{⌊{\nu }_y-v_y⌉}^{1/2}+{\widehat{d}}_y\hfill \\ \hfill {\dot{\nu }}_z& =-g+c_{\varphi }{c}_{\theta }\frac{u_p^{\circ }}{m^{\circ }}-{\lambda }_{z,1}{⌊{\nu }_z-v_z⌉}^{1/2}+{\widehat{d}}_z\hfill \\ \hfill {\dot{\nu }}_1& =\frac{J_2^{\circ }-J_3^{\circ }}{J_1^{\circ }}{\omega }_2{\omega }_3+\frac{J_p^{\circ }}{J_1^{\circ }}{\omega }_p{\omega }_2+\frac{1}{J_1^{\circ }}{\tau }_1^{\circ }-{\lambda }_{\varphi ,1}{⌊{\nu }_1-{\omega }_1⌉}^{1/2}+{\widehat{\delta }}_1\hfill \\ \hfill {\dot{\nu }}_2& =\frac{J_3^{\circ }-J_1^{\circ }}{J_2^{\circ }}{\omega }_1{\omega }_3-\frac{J_p^{\circ }}{J_2^{\circ }}{\omega }_p{\omega }_1+\frac{1}{J_2^{\circ }}{\tau }_2^{\circ }-{\lambda }_{\theta ,1}{⌊{\nu }_2-{\omega }_2⌉}^{1/2}+{\widehat{\delta }}_2\hfill \\ \hfill {\dot{\nu }}_3& =\frac{J_1^{\circ }-J_2^{\circ }}{J_3^{\circ }}{\omega }_1{\omega }_2+\frac{1}{J_3^{\circ }}{\tau }_3^{\circ }-{\lambda }_{\psi ,1}{⌊{\nu }_3-{\omega }_3⌉}^{1/2}+{\widehat{\delta }}_3\hfill \\ \hfill {\dot{\widehat{d}}}_j& =-{\lambda }_{k,2}{⌊{\nu }_j-v_j⌉}^0\hfill \\ \hfill {\dot{\widehat{\delta }}}_i& =-{\lambda }_{k,2}{⌊{\nu }_i-{\omega }_i⌉}^0\hfill \end{array}$$

(14)

which ensure the convergence in finite time of the estimation errors $d_j-{\widehat{d}}_j$, $j=x,y,z$, and ${\delta }_i-{\widehat{\delta }}_i$, $i=1,2,3$, for some ${\lambda }_{k,1},{\lambda }_{k,2}$ with $k=x,y,z,\varphi ,\theta ,\psi$, where ${⌊\sigma ⌉}^{1/2}={\left|\sigma \right|}^{1/2}\mathrm{sign}\left(\sigma \right)$, ${⌊\sigma ⌉}^0=\mathrm{sign}\left(\sigma \right)$ for any scalar signal $\sigma$.

Theorem 1. 

Given the reference ${\chi }_r$, bounded along with its derivatives ${\dot{\chi }}_r$, ${\ddot{\chi }}_r$, the controller (11), (12), with the gains ${\lambda }_{\eta ,1},{\lambda }_{\eta ,2}>0$ with $\eta =x,y,z,\varphi ,\theta ,\psi$, ensures that the error $e_{\chi }$ tends to zero asymptotically and the perturbation estimations ${\widehat{d}}_{\eta }$, in (14) tends to $d_{\eta }$ in finite time, with $\eta =x,y,z,\varphi ,\theta ,\psi$.

Proof. 

Since the subsystems (13a) and (13b), with ${\left(d_x-{\widehat{d}}_x{d}_y-{\widehat{d}}_y{d}_z-{\widehat{d}}_z{\delta }_1-{\widehat{\delta }}_1{\delta }_2-{\widehat{\delta }}_2{\delta }_3-{\widehat{\delta }}_3\right)}^T$ as input, is input-to-state stable (in fact, the dynamics matrix of each subsystem is Hurwitz, and the ISS Lyapunov function is quadratic) and the system, for $\eta =x,y,z,\varphi ,\theta ,\psi$,

$$\begin{array}{cc}\hfill {\xi }_{\eta }& ={\nu }_{\eta }-v_{\eta },{\tilde{d}}_{\eta }=d_{\eta }-{\widehat{d}}_{\eta },\hfill \\ \hfill {\dot{\xi }}_{\eta }& =-{\lambda }_{\eta ,1}{⌊{\xi }_{\eta }⌉}^{1/2}+{\tilde{d}}_{\eta }\hfill \\ \hfill {\dot{\tilde{d}}}_{\eta }& =-{\lambda }_{\eta ,2}{⌊{\xi }_{\eta }⌉}^0\hfill \end{array}$$

converges to the origin in finite time, and the systems (13a) and (13b) is globally uniformly exponentially stable. In fact, considering the part relating to the position on the x-axis ${\chi }_x$ (the demonstrations are analogous for the other variables taken into consideration), it is possible to consider the following Lyapunov function

$$V_{e,x}=\frac{1}{2}{\chi }_x^T{P}_{e,x}{\chi }_x,{\chi }_x=\left(\begin{array}{c}{\chi }_{x,1}\\ {\chi }_{x,2}\end{array}\right)=\left(\begin{array}{c}{⌊{\nu }_x-v_x⌉}^{1/2}\\ d_x-{\widehat{d}}_x\end{array}\right)$$

with

$$P_{e,x}=P_{e,x}^T=\left(\begin{array}{ccc}{\lambda }_{x,1}^2+4{\lambda }_{x,2}& & -{\lambda }_{x,1}\\ -{\lambda }_{x,1}& & 2\end{array}\right)>0.$$

since ${\lambda }_{x,2}>0$. Notice that ${\lambda }_{min}^{P_{e,x}}\parallel {\chi }_x{\parallel }_2^2/2\le V_{e,x}\le {\lambda }_{max}^{P_{e,x}}{\parallel {\chi }_x\parallel }_2^2/2$, with ${\lambda }_{min}^{P_{e,x}}$, ${\lambda }_{max}^{P_{e,x}}$ the minimum and maximum eigenvalues of $P_{e,x}$. To show the convergence, one computes the derivative of the Lyapunov candidate

$${\dot{V}}_{e,x}=-\frac{1}{2{\chi }_{x,1}}{\chi }_x^T{P}_{e,x}{A}_x{\chi }_x=-\frac{{\lambda }_{x,1}}{2{\chi }_{x,1}}{\chi }_x^T{Q}_{e,x}{\chi }_x$$

where

$$Q_{e,x}=\left(\begin{array}{ccc}({\lambda }_{x,1}^2+2{\lambda }_{x,2})& & -{\lambda }_{x,1}\\ -{\lambda }_{x,1}& & 1\end{array}\right)>0$$

$$\begin{array}{cc}\hfill P_{e,x}{A}_x& =\left(\begin{array}{cc}{\lambda }_{x,1}^2+4{\lambda }_{x,2}& -{\lambda }_{x,1}\\ -{\lambda }_{x,1}& 2\end{array}\right)\left(\begin{array}{cc}{\lambda }_{x,1}& -1\\ 2{\lambda }_{x,2}& 0\end{array}\right)={\lambda }_{x,1}{Q}_{e,x}>0\hfill \end{array}$$

since ${\lambda }_{x,1},{\lambda }_{x,2}>0$. Since $\parallel {\chi }_x{\parallel }_2^2=\left|{\xi }_x\right|+{\tilde{d}}_x^2$, clearly $|{\xi }_x{|}^{1/2}\le {\parallel {\chi }_x\parallel }_2\le \sqrt{2}{V}_{e,x}^{1/2}/\sqrt{{\lambda }_{min}^{P_{e,x}}}$, so that

$$-{\displaystyle \frac{1}{|{\xi }_x{|}^{1/2}}}\le -\frac{\sqrt{{\lambda }_{min}^{P_{e,x}}}}{\sqrt{2}}\frac{1}{V_{e,x}^{1/2}}$$

and

$$\begin{array}{cc}\hfill {\dot{V}}_{e,x}& =-\frac{{\lambda }_{x,1}}{2|{\xi }_x{|}^{1/2}}{\chi }_x^T{Q}_{e,x}{\chi }_x\le -{\lambda }_{min}^{Q_{e,x}}{\parallel {\chi }_x\parallel }_2^2\le -\frac{{\lambda }_{min}^{Q_{e,x}}}{{\lambda }_{max}^{P_{e,x}}}2V_{e,x}=-{\gamma }_{min}^{e,x}{V}_{e,x}^{1/2}\hfill \end{array}$$

$${\gamma }_{min}^{e,x}={\lambda }_{x,1}\frac{\sqrt{{\lambda }_{min}^{P_{e,x}}}}{\sqrt{2}{\lambda }_{max}^{P_{e,x}}}{\lambda }_{min}^{Q_{e,x}}$$

where

$${\lambda }_{min}^{Q_{e,x}}=\frac{{\lambda }_{x,1}^2+2{\lambda }_{x,2}+1}{2}-\frac{1}{2}\sqrt{{({\lambda }_{x,1}^2+2{\lambda }_{x,2}+1)}^2-8{\lambda }_{x,2}}.$$

This implies that ${\xi }_x$, ${\tilde{d}}_x$ converge to zero in finite time for ${\lambda }_{x,1},{\lambda }_{x,2}>0$, since

$$\begin{array}{cc}\hfill \frac{\sqrt{{\lambda }_{min}^{P_{e,x}}}}{\sqrt{2}}{\parallel {\chi }_x\left(t\right)\parallel }_2& \le V_{e,x}^{1/2}\left(t_0\right)-\frac{{\gamma }_{min}^{e,x}}{2}(t-t_0)\le \frac{\sqrt{{\lambda }_{max}^{P_{e,x}}}}{\sqrt{2}}{\parallel {\chi }_x\left(t_0\right)\parallel }_2-\frac{{\gamma }_{min}^{e,x}}{2}(t-t_0)=0\hfill \end{array}$$

for $\overline{t}=t_0+\sqrt{2{\lambda }_{max}^{P_{e,x}}}{\parallel {\chi }_x\left(t_0\right)\parallel }_2/{\gamma }_{min}^{e,x}$. □

The control scheme is divided into an inner and outer loop design. The inner loop design is responsible for attitude and altitude control, while the outer loop is responsible for controlling lateral and longitudinal motion. The desired trajectories for translation and heading movements are provided externally, while the roll and pitch trajectories are generated from the outer loop. The functional scheme of the proposed controller is the one shown in Figure 2.

4. Simulation and Results Analysis

The performance of the proposed controller is discussed in this section. First of all, the simulation environment used is briefly described in Section 4.1; then, the main wind models used in aerospace research are briefly recalled in Section 4.2, and specifically, the one used in the simulations of this paper is described. Finally, the results obtained by the robust controller in a generic flight mission are shown and compared in Section 4.3 with the non-robust controller in the presence of disturbances and parameters uncertainty present in the physical model of the chosen simulation environment.

4.1. Co-Simulation Environment with Matlab-Simulink and Simcenter Amesim with FPGA Implementation

The software Simcenter Amesim (see [15]), a multi-physics system simulation tool developed by Siemens Digital Industries Software, is used to model the different subsystems of the UAV: batteries, propulsion chain and flight dynamics. In general, the tool provides off-the-shelf components available in libraries, covering several physical (fluids, mechanical, electrical, thermal…) and application (aerospace, automotive, gas turbines…) domains. Those subsystem plant models have been favorably compared against experimental data provided by the manufacturers first. The use of this tool was therefore aimed at validating the controller and the estimators, which in the first instance were developed and tested only in the Matlab–Simulink environment. The overall plant model was then integrated in a co-simulation framework capable of modeling drone’s navigation sensors, mission environment and algorithms to simulate the drone’s behavior under different scenarios. The model present in Simcenter Amesim is complex with respect to the reference mathematical model considered for the design of the controller presented in Section 2; therefore, it contains both unmodeled dynamics and parameters that are time-varying and is therefore a valid test bed to evaluate controller’s performance. In the present research, the UAV model considered is seen as a black box model, and the parameters used in the model-based controller have been preliminary identified with Least Square Method. For example, at the modeling level, the aerodynamic coefficients such as the thrust and drag coefficient play an important role in the design of the controller. In fact, several studies have shown that they are variable based on the operating conditions of the drone. In both the physical models used in the advanced, realistic simulator used in this paper (see [15]) and in experimental studies, they are determined by multivariable functions (see [16,17]) such as orientation, angular and total velocity as well as the geometric characteristics themselves. For these reasons, they are parameters that can be identified in a static way, and this is the practice that can be followed, but they are actually time-varying during the simulation. For this reason, it is necessary to take them into account in the design phase of the controller in order to make it robust, and in our case, all these variations in propulsion level are also considered in the disturbance term. The drone simulation parameters used in the model-based design controller are presented in

Table 1. Nominal quadrotor parameters.

m	Mass of the Airframe	5 kg
l	Distance of CoG to the rotor shaft	0.3 m
$J_1$	Inertia in the x-axis	0.011521 kg $·$ m${}^2$
$J_2$	Inertia in the y-axis	0.0362132 kg $·$ m${}^2$
$J_3$	Inertia in the z-axis	0.029142 $·$ m${}^2$
$J_p$	Inertia of the propellers	0.0003 $·$ m${}^2$
g	Gravity acceleration	9.81 m $·$ s${}^{-2}$
b	Trust factor	4.5625 $\times$ 10${}^{-5}$ N $·$ s${}^2$$·$ rad${}^{-2}$
c	Drag factor	1.375 $\times$ 10${}^{-5}$ N $·$ s${}^2$$·$ rad${}^{-2}$ $·$ m

.

In addition to the use of the physical model of the Siemens drone with identical performance to the real one, a perfect candidate for this application is the FPGA implementation, since it can provide strict real–time execution, high power efficiency, small size, and high reliability, at the cost of a reduced dynamic range. Specifically, the Intel Cyclone 10 LP FPGA was available and corresponded to what was needed; then, it has been used as a reference, and it comes on a minimal board, which is useful for actual power consumption tests and fast prototyping, allowing installation on small UAVs for real-world tests. Without going into the details of the FPGA implementation, Figure 3 shows the structure of the implemented controller, which is only part of a complete system architecture. During the implementation on FPGA, the typical problems related to the operations of trigonometric function, their inverse, reciprocal or division, square root, computation of the derivatives, discrete approximation, variables normalization, as well as to the fixed-point arithmetic of the controller and input variables have been addressed. The simulations presented were carried out in a co-simulation environment between Matlab and Simcenter Amesim in order to have a complex model on which to test the controller and were implemented in FPGA at the end.

4.2. Environmental Disturbances Acting on the Quadrotor: Wind Modeling

The categories, features, and mathematical models of the wind, which have a big influence on UAV in the low-altitude condition, can be classified in four groups: constant wind, turbulent flow, many kinds of wind shear, and the propeller vortex (see [18]). The attitude, speed, and position of UAV are sensitive to wind disturbance, and consequently, the wind is seen as a parameter that has the bigger influence on efficiency. The main features of each wind category are shortly evoked below:

• Constant wind deals with the average wind speed in a specified environment. Constant wind does not exist in nature, because it is only the reference value of the wind speed in a given environment. To perform a more realistic simulation test, other types of wind need to be considered. Generally through statistical data, it is possible to obtain its value in the different temporal or spatial conditions. The numerous relevant websites that log wind data are supportive in this regard. Usually, constant wind is used in simulation tests of UAVs at first, but it is unable to reproduce the flight environment correctly.
• Turbulent flow is a continuous random change, which is always conducted by constant wind. The source of turbulent flow is coupled to many elements, such as wind shear, heat exchange, topographical factors, etc. (see [19]). Stochastic Process Theory is commonly used to describe atmospheric turbulence in engineering applications. Turbulent flow models incorporate the Dryden model and the Von Karman model, both of which depend on a vast amount of measurements and statistics (see [20]). The difference is that Dryden’s model finds the correlation function of turbulence before giving rise to the spectral function, while Von Karman, on the contrary, first establishes the spectral function and then deduces the correlation function of turbulence. From the aforementioned literature, it emerges that between the two models, there is no relevant difference for the application use in the simulation phase; therefore, both can be used within engineering problems designed to simulate the wind.
• Wind shear The American National Research Council defines the degree of wind shear as the ratio of the difference between two wind vectors at two points and the distance between those two points. There are many sources of wind shear. In UAV flight, common low-altitude wind shear contains frontal wind shear, topography-induced wind shear, microdownburst and night jet wind shear. The frequency of the wind shear is also not as large as the turbulent flow, it changes every few seconds as its maximum frequency. The speed of the wind shear is a function of the spatial position, and its model is relatively simple. In the simulation test, the abrupt change in constant wind speed is used to indicate the occurrence of wind shear. Wind shear is a discrete or deterministic wind speed, which often occurs in a very short time, and it is a strong atmospheric perturbation. The wind shear model can be divided into several categories based on its profile geometry, including the rectangular model, the trapezoidal model, and the ’1-consine’ model.
• Propeller vortex is related to large manned aircraft. There are wake vortex caused by aircraft wings that significantly affect other aircraft. As for UAVs, due to their slower flight speed, the propeller vortex will no longer play an important role, while here, another vortex caused by the propellers could be considered, just as the name suggests, but it can be kept in consideration in the modeling phase of the system to be controlled.

The Principle of Wind Affects the UAV

All kind of wind modify the speed of the UAV, as a result altering the angle of attack and the angle of drift, which influences the values of the aerodynamic derivatives and eventually causes disturbances to the outcome force and torque of the UAV. In general, the most common approach to relate the consequences of winds on UAVs is from the point of view of speed. The wind models cited above explain only the distribution of velocity in space and time. Unfortunately, under most conditions, the induced speed of a certain wind field is not uniform, such as the wind gust field and the vortex field of the propeller, which denotes that the speed is not a constant value and changes together with the different parts of the UAV. To figure out this problem, in [21], we proposed an averaging way to calculate the actual wind speed and the wind gradient in order to derive the relationships between the variation of the air speed and the actual wind as well as between the variation of the body angular velocities and wind gradient. On the other hand, from the point of view of force, the effects of the non-uniform part of the flow field are considered an external force on the UAV. It is a well-known fact that interaction forces are the essential reason for the change in the kinetic states of objects, and the force caused by air is called resistance. In the simulation results presented, the Dryden wind turbulence model is used. The Dryden model considers the linear and angular velocity components of continuous gusts as spatially varying stochastic processes and specifies each component’s power spectral density. The Dryden wind turbulence model is characterized by rational power spectral densities, so exact filters can be designed that take white noise inputs and output stochastic processes with the Dryden gusts’ power spectral densities. The trends of the wind disturbances considered in the subsequent simulation are shown in Figure 4, dividing the contributions with respect to the three reference axes by linear and angular speed. For more mathematical details on the considered model, it is possible to refer to [22].

Table 2. Dryden Wind Turbulence Model Parameters.

Model	Continuos Dryden (+q −r): Use continuous
Type	representation of Dryden velocity
	spectra with positive vertical
	and negative lateral angular rates spectra.
Altitude	z
UAV	$\sqrt{v_x^2+v_y^2+v_z^2}$
Velocity
Measured	15 m/s.
wind speed	Measured wind speed at a height of 6 m
speed
Wind	0
direction	Measured wind direction at a height of 6 m
	(degrees clockwise from north)
	to aid in transforming the low-altitude
	turbulence model into a body coordinates.
Turbulence	${10}^{-2}$ -Light
intensity	Probability of the turbulence intensity being exceeded
Wingspan	10 m
	required in the turbulence calculation on the angular rates.
Band limited	0.01 s
noise sample	Noise sample time at which
time	the unit variance white noise signal is generated
Noise seeds	(23,341 23,342 23,343 23,344)
(ug vg wg pg)	Random noise seeds, specified as a four-element vector one for each
	of the three velocity components
	and one for the roll rate.

shows the main parameters and the values used for the generation of disturbances with reference to the Simulink block called “Dryden Wind Turbulence Model (Continuous)”. The gains value of the HOSM controller and of the disturbances estimators used in the simulations shown are reported in the

Table 3. Parameters of HOSM Controller and Estimators.

${\mathit{k}}_{\mathit{x},0}=3$	${\mathit{k}}_{\mathit{x},1}=10$	${\mathit{k}}_{\mathit{x},2}=3$	${\mathit{k}}_{\mathit{y},0}=3$	${\mathit{k}}_{\mathit{y},1}=10$	${\mathit{k}}_{\mathit{y},2}=2$
$k_{z,0}=3$	$k_{z,1}=10$	$k_{z,2}=2$	$k_{\varphi ,0}=1.5$	$k_{\varphi ,1}=5$	$k_{\varphi ,2}=1$
$k_{\theta ,0}=1.5$	$k_{\theta ,1}=5$	$k_{\theta ,2}=1$	$k_{\psi ,0}=2.7$	$k_{\psi ,1}=9$	$k_{\psi ,2}=1.8$
${\lambda }_{x,1}=10$	${\lambda }_{x,2}=3$	${\lambda }_{y,1}=10$	${\lambda }_{y,2}=3$	${\lambda }_{z,1}=10$	${\lambda }_{z,2}=3$
${\lambda }_{\varphi ,1}=2$	${\lambda }_{\varphi ,2}=1$	${\lambda }_{\theta ,1}=2$	${\lambda }_{\theta ,2}=1$	${\lambda }_{\psi ,1}=2$	${\lambda }_{\psi ,2}=1$

. We want to emphasize that controller parameters give a formal solution of the stated control problem. Nevertheless, in practice, one often requires to calibrate the convergence rate, either to slow it down (relaxing the burden on the actuators) or to accelerate it in order to meet some system requirements. Consider, in this context, redundantly augmenting the magnitude parameter does not accelerate the convergence but rather only increases the chattering, while its reduction may lead to a loss of convergence.

4.3. Simulation Results

In the test mission shown, we consider a rest-to-rest maneuver of the quadrotor with an initial state $\left[xyz\psi \right]=\left[0000\right]$ and a final state $\left[xyz\psi \right]=[4515\pi /4]$. In fact, we assumed zero linear and angular velocity of the UAV at the boundary states, with zero body orientation at time $t=0$. The reference generator of the trajectories to be tracked is generated in such a way to have minimum accelerations and with a time suitable for carrying out the mission (see [23] for more details). Figure 5, Figure 6, Figure 7, Figure 8, Figure 9, Figure 10, Figure 11 and Figure 12 are representative of the simulation considered in the cosimulation environment with FPGA implementation described in (Section 4.1): they compare and show the tracked characteristic state variables between the robust controller with disturbances compensation and the non-robust controller and the errors, the control inputs, the disturbances $d_x$, $d_y$, $d_z$, ${\delta }_3$ and their estimation ${\widehat{d}}_x$, ${\widehat{d}}_y$, ${\widehat{d}}_z$, ${\widehat{\delta }}_3$; the effectiveness of the robust controller compared to the non-robust one and good estimation of the disturbances is evident. In fact, if the wind disturbance is not compensated, it takes the drone away from the trajectory to follow and leads it to dangerous states in terms of dynamics and greater energy consumption due to the removal from the desired trajectory. In particular, Figure 5, Figure 6 and Figure 7 show the trends of the longitudinal (x), latitudinal (y), and altitudinal (z) positions, respectively, and their errors: in order to make evident the tracking difference between the robust and non-robust controller, their trends have been separated as well as their errors. The tracking errors of the non-robust controller are always large due to disturbances and parametric variations, but in the final part of the simulation, it cannot even keep the drone stable. With the same philosophy at the graphic level, Figure 8 shows the trend of yaw rate for the two controllers: the results are similar; therefore, with better performance than the robust controller, it should be emphasized that unlike the previous figures in this case, probably due to the wind and the unknown time-varying dynamics of some coefficients present in the simulator, the error reaches values of 12 degrees while keeping it always bounded.

Figure 9 shows the control inputs of robust control expressed as angular velocities of the four propellers. Modern-day ESC protocols can communicate at speeds of 37.5 kHz or greater (in general, it is between 30 and 60 kHz because of the relatively low inductance of the high-speed motors and the potential interference to sensor boards), see [24]. The frequency of the actuation signals of the robust controller is $100\mathrm{Hz}$ and thus abundantly meets the requirements of modern commercial drones. The control inputs in the robust case have more rapid variations, while in the non-robust case, they assume a smooth trend: this is due to the terms of the estimators trying to compensate for what they see as possible parametric disturbance or uncertainty. Figure 10 shows the three-dimensional spatial trajectories of the two controlled drones in order to be able to visually quantify the error that occurs when the controller is not made robust with respect to disturbances. It is evident that the error that non-robust control makes at the 3D trajectory level is high, and therefore, there is an increase in consumption and the risk of drone instability. Finally, the disturbances $d_x$-$d_y$-$d_z$-${\delta }_3$ and their estimates used in robust control, together with their respective errors, are shown in Figure 11 and Figure 12. The estimates are good for all disturbances: we note the presence of noise and a slightly larger error in ${\delta }_3$, which is an event that causes errors of some degree in the tracking of the yaw rate, error due to the complex simulator considering that it has variable parameters over time, dynamics not considered in our model, and physical models based on experimental data; in fact, in the simulations not shown due to the lack of space made with the plant in Matlab, this does not happen.

To compare the performance of two proposed controllers in the considered mission, the Mean Absolute Error (MAE), Mean Square Error (MSE), Integral Squared Error (ISE), Integral Absolute Error (IAE), Integral Time-Weighted Absolute Error (ITAE) and Integral Time Square Error (ITSE) measures have been considered as comparative indices for the tracking along x, y, z and $\psi$. These performance indices are defined, for the variable x (the argument is analogous for the other variables) with respect to the reference $x_{ref}$ to be tracked, during the N sampling instants, as follows

$$\mathrm{MAE}=\frac{{\sum }_{i=1}^N[x_{ref,i}-x_i]}{N}$$

$$\mathrm{MSE}=\frac{1}{N}\sum _{i=1}^N{(x_{ref,i}-x_i)}^2$$

$$\mathrm{ISE}=\int {(x_{ref}-x)}^{2}dt$$

$$\mathrm{IAE}=\int \mid x_{ref}-x\mid dt$$

$$\mathrm{ITAE}=\int (t\mid x_{ref}-x\mid )dt$$

where $x_{ref,i}$ and $x_i$ are, respectively, the values of the variable $x_{ref}$ and x at sampling time i, and ISE, IAE and ITAE integrate over time. The results are reported in

Table 4. Performance indices for the x, y, z and $\psi$ tracking errors.

Controller	Controlled Variable	MAE	MSE	ISE	IAE	ITAE	ITSE
Non-robust	x	0.4312	1.3198	12.92	4.2777	35.9148	125.34
Robust	x	0.018	0.000063	0.0063	0.1798	0.9743	0.0339
Non-robust	y	0.3867	1.3198	7.0007	3.8442	32.3692	66.1772
Robust	y	0.036	0.0019	0.0196	0.3604	2.2055	0.127
Non-robust	z	0.216	0.1321	1.2969	2.1505	13.9362	11.0271
Robust	z	0.0124	0.000021	0.0021	0.1238	0.5591	0.0089
Non-robust	$\psi$	19.3308	1492.2	1475.3	1.9254	1655.83	138,009
Robust	$\psi$	3.0647	17.5	175.7	30.6552	1.3336	629

, which shows the difference for each performance index of the tracked variables between the non-robust and robust controller.

Finally, the control effort for each of the two controllers is also reported in

Table 5. Actuation effort of robust and non-robust controllers as propeller angular velocity square integral along all the mission time.

Controller	Non-Robust	Robust
${\omega }_{p,1}$	2.535937464 × ${10}^8$	2.522437042 × ${10}^8$
${\omega }_{p,2}$	2.534601465 × ${10}^8$	2.512100166 × ${10}^8$
${\omega }_{p,3}$	2.533127333 × ${10}^8$	2.514066778 × ${10}^8$
${\omega }_{p,4}$	2.534819281 × ${10}^8$	2.515695169 × ${10}^8$
Total	1.013848554 × ${10}^9$	1.006429915 × ${10}^9$

. It is computed as the time integral of the square of the angular velocity for each propeller. The effort of non-robust control is greater, but above all, it has not fulfilled the mission it was supposed to carry out, and therefore, additional time and energy will be required. In conclusion, Table 4 and Table 5 point out, respectively, that the robust controller has better performance in terms of tracking error and lower total actuation effort.

Future work will deal with a comparison with another control strategy such as a neuroadaptive learning algorithm for constrained nonlinear systems with disturbance rejection [25] and asymptotic tracking with novel integral robust schemes for mismatched uncertain nonlinear systems [26].

5. Conclusions and Future Works

This paper presents a dynamic controller for a quadrotor where the user has partial knowledge of the model and in the presence of environmental disturbances such as wind along the linear and angular velocity. The proposed PI-like controller is able to compensate these issues by means of the use of HOSM estimators. The performance of such a controller has been tested in a co-simulation integrated environment between Matlab-Simulink and Simcenter Amesim with FPGA implementation. The results obtained by the robust controller are better than a basic version of the control both in terms of desired trajectories tracking and of actuation effort.

Future investigations will concern the development of theoretical results concerning the digital implementation of the proposed controller by exploiting the methodologies provided, for instance, in [27,28], experimental simulation and multi-agent (quadrotor) control systems in the presence of delay, quantization, and transmission errors (see [29]).