Aerodynamic Effects Compensation on Multi-Rotor UAVs Based on a Neural Network Control Allocation Approach

I. Introduction

IN the last decade, there has been a substantial growth in the usage of unmanned aerial vehicles (UAV) in a series of applications [1]. They are widely used to perform difficult or dangerous missions, military ones [2], or to make civilian tasks more efficient, like package delivery. The research of modeling and control techniques for this type of vehicles, especially the quadrotors, has also grown considerably [1]. A quadrotor dynamic model is usually an under-actuated nonlinear system [3]. Quadrotor UAVs are considered classic unstable dynamical systems [4], hence, the study of suitable algorithms to improve stability and reduce the influence of disturbances is indispensable.

A control allocation algorithm aims at computing a control input that ensures that the commanded virtual control is produced by the set of effectors and actuators at all times. This objective would fail to be met if the commands require forces beyond the capabilities of the actuators due to saturation or other physical limitations such as aerodynamic effects [5].

The use of control allocation offers the advantage of a modular design where the motion control algorithm can be designed without detailed knowledge about the actuators. Issues such as input saturation, rate constraints, actuator fault tolerance, power efficiency, etc., could be managed within the control allocation algorithm [5]. The resulting control laws specify which total control effect to produce but not how to produce it, and hence the need for control allocation arises [6].

Conventional multirotor UAVs have the rotors in parallel directions to collectively counteract gravity and optimize flight time. Consequently, these aircrafts have under-actuated dynamics due to the coupling between the horizontal translational and rotational dynamics [7]. As a result, the control of a quadrotor tends to present a high sensibility when subjected to an input constraint [8]. Therefore, control allocation is necessary for setting suitable commands for each rotor, so that it generates the desired control effect, that is, the forces and torques, while presenting an acceptable performance and respecting the constraints imposed by the minimal and maximal angular velocities and hence the maximum and minimum forces and torques [8].

In the classic quadcopter modeling used for most control techniques, the control allocation is often performed through the control effectiveness matrix, usually named mixing matrix or allocation matrix. It provides an invertible relation between mechanical conjugates and motor rotation speeds. This matrix results from a model simplification that neglects important aerodynamic effects that influence thrust forces and torque behavior, especially regarding small aircraft or aggressive maneuvering [9]. As a matter of fact, this negligence leaves the quadcopter subject to saturation, and the torque and thrust commands sent by the controller are not properly generated by the actuators, because the motors cannot achieve the necessary rotation speed.

In this sense, a new control allocation method was proposed to consider aerodynamic effects, most notably those regarding blade flapping and inflow effects as described in [10], based on a function fitting neural network to replace the classic mixing matrix. This method allows us to improve the controller performance without the need for directly using aerodynamic effects equations in the control algorithm or the allocation problem. Furthermore, network training is performed offline, and the in-flight control allocation depends only on the virtual control inputs and some state-space variables, which will be detailed further in this paper.

While it is true that there are several results regarding UAV control in the literature, we did not find works approaching the problem of the error inserted by the disregard of aerodynamic effects in the control allocation, which is the target of this paper. Even the robust controllers available in the literature and mentioned in the next section, still use the classic mixing matrix, forcibly inserting an error in the rotation speeds established at the motors and hence in the control performance.

Thus, the main contributions of this paper are:

1) We propose the use of a classic neural network as a novel control allocation approach that takes into account the aerodynamic effects.

2) The proposed method can be used to improve the performance of quadrotor controllers in general, since it removes the error introduced by the classic allocation matrix.

3) Our experimental results verified that neglecting aerodynamic effects in the control allocation has a direct negative impact in the flight performance;

4) Several aspects of the system were improved with the new control allocation, including the decrease in the allocation error, the increase in the controller performance, and a better overall signal quality.

II. Related Work

Dynamic models with minimum complexity for quadrotors are well documented in the literature [11]. These models generally use six-degree-of-freedom (6-DOF) rigid-body equations and a highly simplified rotor dynamics as a function of the square of motor speed. The treatment of quadcopter dynamics around steady-state conditions has often ignored rotorcraft aerodynamic effects due to its complex physical modeling [12]. In certain circumstances similar to hover flight, these effects could be neglected. Nevertheless, they can prove to be very relevant during agile maneuvers [9] or in the case of small aircrafts [13], [14]. Davis and Pounds [13] sought to understand and exploit the aerodynamic effects acting on a small-scale quadrotor to regulate relative position, attitude, and velocity without extrinsic sensors. Ning et al. [14] investigated the aerodynamics and aeroacoustics of propellers for small unmanned aerial vehicles. Powers et al. [15] demonstrated that the dynamic response and performance of a micro UAV are greatly influenced by its aerodynamics, and the introduction of these effects into controllers is also discussed. They show that the controller performance can be improved when considering aerodynamic effects into the model.

Pounds et al. [16], [17] have included some aerodynamic effects modeling in the rotor dynamics equations to enhance model accuracy. The authors adopted a generalized flapping dynamics structure, based on helicopter theories, to model the rotor flapping phenomenon. Bristeau et al. [18] studied the aerodynamic role of the propellers. Hoffman et al. [19] and Huang et al. [9] modeled aerodynamic effects associated with blade flapping, airflow distribution, and thrust generation, to derive a model that can cover extreme maneuvers.

Cai et al. [11] introduced a frequency-domain modeling methodology that can be applied to various multi-rotor aerial vehicles. Niemiec and Gandhi [20] created a dynamic simulation for a 2 kg quadcopter, where several inflow models were used to analyze the behavior of the aircraft in trim, in hover, and forward flight. Kantue and Pedro [12], [21] used a radial basis functions neural network (RBF-NN) for model identification of an unmanned quadcopter in accelerated flight to estimate rotor dynamics parameters from an unknown flapping dynamics model.

The works presented above are mainly about modeling the aerodynamic effects that impact multirotor aircraft systems. However, the literature is still lacking in regard to techniques to actually take these effects into account in the quadrotor system, especially when it comes to the control allocation, which none of the previous works mention.

Omari et al. [22] present the design and evaluation of a nonlinear control scheme for a six-rotor helicopter that takes the first-order drag effects into account explicitly, and a hierarchical nonlinear controller is designed to actively compensate for the nonlinear effects caused by these drag forces. Note that, even though a dynamic model including the blade flapping and induced drag forces is presented in this paper, the corresponding traditional mixing matrix for a six-rotor helicopter is still employed by the authors. Even if the controller takes into account the drag forces, the control allocation completely neglects their existence, jeopardizing some of the controller’s performance when converting its commands to the actuators’ rotation speeds.

Chebbi et al. [23] use a polynomial regression to obtain the control allocation of a coaxial multirotor. However, the problem the authors are trying to solve in this work differs from the one we attack in our paper. They try to estimate a model that will solve the interference between the propellers in a pair of coaxial motors. Again, even though they estimate a thorough model of the actuators’ nonlinear effects, they neglect them in the control allocation.

Johansen and Fossen [5] and Alves [8] present the fundamentals of control allocation theory. A survey regarding different control allocation algorithms and applications to aerospace, maritime, automotive, and other application areas is presented in [5], while the applications of several methods of control allocation to an experiment with a differential-drive mobile robot and simulations of a quadrotor UAV and a Remotely Operated Underwater Vehicle are presented in [8].

Most recently, the interest in control allocation applied to multirotor vehicles has increased, given the rise in the usage of this type of aircraft. Kirchengast et al. [24] showed the comparison of different control allocation techniques, as well as Alves [8], but this time they are applied to a laboratory setup which emulates the rotational dynamics of a quadrotor, though not to a real quadcopter. Monteiro et al. [25] also analyzed control allocation algorithms to countermeasure saturation in the actuators, however, they generalized the classic direct control allocation (DCA) to suggest new methods, called partial control allocation (PCA) [25].

Doman and Oppenheimer [26] stated that forces and moments produced by the vehicle’s aerodynamic control surface are often nonlinear functions and that this limits most control allocation algorithms applications, since they are based on the assumption that the control variable rates are linear functions. The authors then presented a mixed optimization control allocation scheme that minimizes the norm of the error and the deviation of the control input from a preferred value, to attack the inaccuracies introduced by these linear assumptions.

In addition, Smeur et al. [27] dealt with the problem of actuator saturation for Incremental nonlinear dynamic inversion (INDI) controlled flying vehicles. The authors integrated the Weighted Least Squares optimization-based control allocation algorithm into INDI, which allows for prioritization among roll, pitch, yaw, and thrust. The conventional application area for control allocation algorithms is over-actuated mechanical systems, however, the principles of control allocation are general and can be applied to under-actuated systems subject to saturation and nonlinear effects.

The main technique most control allocation algorithms have in common is solving a constrained optimization problem [28], and this process is usually performed on-line [29], requiring computational resources from embedded systems [5]. Even though their computational complexity is mostly within the capabilities of today’s computer technology, some works highlight the advantages of requiring a low computational power on-line [30] such as the method we propose in this paper.

The control allocation theory establishes that it can be used to treat actuator saturation problems or other disturbances which may affect their optimal functioning [5]. However, at the moment of this research, there are no works proposing to specifically compensate aerodynamic effects inside the control allocation, replacing the traditional control effectiveness matrix and thus mitigating the allocation error.

Neural networks have been receiving wide attention in the control field, being proposed especially for identification and control of nonlinear dynamic plants [31]. A common application of neural networks to controllers is as an adaptive element to dynamic inversion plants, as stated by Cao et al. [32], Jiang et al. [33], and Wang et al. [31]. Likewise, Ansari et al. [34], [35] employed a baseline controller with a dynamic inversion in the inner loop for attitude tracking and body rates stabilization, with the addition of a Radial Basis Function Neural Network to estimate the unknown attitude dynamics of the quadrotor.

According to Achermann et al. [36], small UAVs are particularly susceptible to wind drag because of their relatively low mass and their operation at low altitudes where the wind environment can be more complex. The authors proposed an approach for predicting high-resolution wind fields using a deep convolutional neural network to generate 3D wind estimates. In a similar topic, Shi et al. [37] presented a deep-learning-based robust nonlinear controller called Neural-Lander that improves control performance of a quadcopter during landing, since near-ground trajectory control is difficult for multi-rotor drones, due to the aerodynamic effects caused by interactions between multi-rotor airflow and the environment. Verberne and Moncayo [38] described the design of an extension of the Nonlinear Dynamic Inversion control architecture for wind disturbance rejection in quadrotor UAVs, followed by the application of adaptive artificial neural networks to correct inversion errors caused by wind disturbance.

Neural networks have also been used for fault-tolerant control design in quadcopters, as in the works of Abbaspour et al. [39] and Song et al. [40]. The former used a neural network adaptive structure for fault detection and isolation combined with a nonlinear dynamic inversion technique to compensate for faults in the actuators, while the latter combined the radial basis function neural networks with virtual parameter estimating algorithms to develop an indirect neural network-based adaptive fault-tolerant control scheme. Finally, Reddinger and Gandhi [41] trained a predictive neural network to estimate power as a function of the redundant control settings for a range of flight speeds for a compound helicopter. This neural network is used as a surrogate model for a gradient-based optimization to allocate the redundant control settings with a minimized power requirement.

Note that even though neural networks are used to identify aerodynamic models, the application of a neural network to perform the control allocation of a quadrotor control system is still uncommon in the literature, specially considering aerodynamic effects, which is an application our research did not find in the most relevant works in this area of study.

III. Modeling and Control of a Quadcopter

The traditional quadrotor model assumes that the aircraft is built as a rigid structure, usually in form of X, with four arms, as shown in Fig. 1, where the axis convention used in this paper is presented. The angles use the Tait-Bryan convention, where the movement of rotation around the $x$-axis is called roll and is represented in this paper by the angle $\phi$. The roll motion induces a linear displacement of the aircraft to the sides. The rotation around the $y$-axis is called pitch, represented by the angle $\theta$, and causes a linear displacement in the forward or backward direction. The rotation around the z-axis is called yaw, represented by the angle $\psi$. The torques that induce each one of these motions are here called $\tau_\phi$, $\tau_\theta$, $\tau_\psi$, respectively. The force that enables the drone to move upwards is the $thrust$ $T$. Ideally, it would only have one vertical component in the z-axis.

Figure  1.  In the front, we show the axis configuration in the local frame, as well as the rotation sense for each rotor and the roll ($\phi$), pitch ($\theta$), yaw ($\psi$) notations used in the quadcopter model. In the back, there is the axis configuration in the inertial frame.

DownLoad: Full-Size Img PowerPoint

Since our system is under-actuated, a control allocation algorithm is needed to send the appropriate commands to the actuators. The controller commands virtual control inputs [5] for the roll, pitch, and yaw torques, as well as the thrust force. The control allocation then transforms these mechanical conjugates into angular velocities commands for the four motors.

Table I shows the influence of the angular velocities of the motors in the resulting motion of the aircraft. Here, $\omega_i$, $i = 1, 2, 3, 4$, is the angular velocity of each rotor. An upward arrow means an increase in the velocity of a specific motor, while a downward arrow corresponds to a decrease in velocity. Note that the arrows can be inverted for a change of direction of motion.

Table  I.  Rotor Speed Influence on Quadcopter Motion

Type of motion	$\omega_1$	$\omega_2$	$\omega_3$	$\omega_4$
Thrust ($T_{\textit{z}}$)	$\uparrow$	$\uparrow$	$\uparrow$	$\uparrow$
Roll ($\phi$)	$\downarrow$	$\uparrow$	$\uparrow$	$\downarrow$
Pitch ($\theta$)	$\downarrow$	$\downarrow$	$\uparrow$	$\uparrow$
Yaw ($\psi$)	$\downarrow$	$\uparrow$	$\downarrow$	$\uparrow$

 | Show Table

DownLoad: CSV

A Method Overview and Notation

In our control allocation, the objective is to produce actuator commands that will generate resultant torques and thrusts on the UAV body as close as possible to the requested control signals, also called virtual control inputs. In this sense, we provide here an overview of the notations used throughout this paper, for better understanding of the upcoming sections.

Definition 1: The virtual control inputs $[\mathbf{T}_c, {{\boldsymbol\tau}}_c] = [ T_c \;\; \tau_{\phi_c} \;\;\tau_{\theta_c} \;\; \tau_{\psi_c} ]^T \in \mathbb{R}^4$ are the commanded signals at the output of the controller, where $\mathbf{T}$ represents the thrust force, and ${\boldsymbol\tau}$, the torques.

Remark 1: In each of the four actuators $i$, $i = 1, 2, 3, 4$, there are the body-frame thrusts, $\mathbf{T}^B_i \in \mathbb{R}^3$, and torques, ${\boldsymbol\tau}_i \in \mathbb{R}^3$, as well as their components $[T_{x_i} \;\; T_{y_i} \;\; T_{{\textit{z}}_i} ]^T$ and $[ \tau_{\phi_i} \;\; \tau_{\theta_i} \;\; \tau_{\psi_i} ]^T$, respectively.

Definition 2: $\mathbf{u} = [ \omega_1 \;\;\omega_2 \;\; \omega_3 \;\; \omega_4]^T \in \mathbb{R}^4$ is the vector of angular velocities sent as actuator commands by the control allocator. These velocities should produce the desired resultant thrust and torques.

Definition 3: The vector $[\mathbf{T}, {{\boldsymbol\tau}}] = [T \;\;\tau_{\phi} \;\; \tau_{\theta}\;\; \tau_{\psi} ]^T \in \mathbb{R}^4$ represents the real resultant mechanical conjugates established at the UAV body. In this paper, they were calculated by means of experimental data and dynamic model equations. Again, $\mathbf{T}$ represents the thrust force, and ${\boldsymbol\tau}$, the torques.

Assumption 1: If $F$ is a virtual function that describes the transformation of the actuator commands $\omega_i$ to real thrust and torques, and $\widetilde{F}$ is its model, we have

Hence, the desired control allocation must perform the inverse operation, using the virtual control inputs

where $\widetilde{F}^\ast$ represents the function which performs the inverse operation of $\widetilde{F}$.

Using this notation, the control allocation algorithm goal is to achieve $\mathbf{T} = \mathbf{T}_c$ and ${\boldsymbol\tau} = {\boldsymbol\tau}_c$. Having $\begin{bmatrix} \mathbf{T}\\ {\boldsymbol\tau} \end{bmatrix} = F(\widetilde{F}^\ast(\mathbf{T}_c, {{\boldsymbol\tau}}_c))$, the composition of the inverse mappings with the real system behaviour functions must be equal to identity [23]. See Fig. 2.

Figure  2.  Data flow from the virtual control inputs $\mathbf{T}_c, {{\boldsymbol\tau}}_c$ to the real mechanical conjugates applied on the UAV body $\mathbf{T}, {{\boldsymbol\tau}}$. $\widetilde{F}^\ast$ can represent the control allocation block.

DownLoad: Full-Size Img PowerPoint

Remark 2: $\widetilde{F}^\ast$ can be interpreted as the control allocator. In this paper, it can refer to the classic mixing matrix or the proposed neural network, and we will compare the performance of both methods.

Remark 3: $\widetilde{F}$ represents the quadcopter dynamic model considering aerodynamic effects such as blade flapping and will be further detailed in Sections III-B and III-C of this paper.

B Quadcopter Model

The dynamic and kinematic behaviors of the UAV can be described by the following equations [42]:

where $g$ represents the gravitational acceleration, $m$ is the total mass of the vehicle, $\mathbf{T}^B$ is the total thrust acting on the center of mass, $\mathbf{J}$ is the inertia matrix, ${\boldsymbol\tau} = [\tau_\phi \;\; \tau_\theta \;\; \tau_\psi]^T$ refers to the vector of torques about the three axes, and $\mathbf{F_d}$ represents the disturbances caused by drag forces. The state variables are:

1) The linear position vector $\mathbf{r}(t) = [ x(t) \;\; y(t) \;\; {\textit{z}}(t) ]^T \in \mathbb{R}^3$;

2) The angular position vector $\boldsymbol\Psi(t) = [ \phi (t) \;\; \theta (t) \;\; \psi (t) ]^T \in \mathbb{R}^3$;

3) The linear velocity vector $\mathbf{v}(t) = [ v_x (t) \;\; v_y (t) \;\; v_{\textit{z}} (t) ]^T \in \mathbb{R}^3$;

4) The angular velocity vector in the body-frame ${\boldsymbol\Omega}(t) = [ P(t) \;\; Q(t) \;\; R(t) ]^T \in \mathbb{R}^3$.

Note that the linear position $\mathbf{r}$, the angular position $\Psi$, and the velocity of the center of mass $\mathbf{v}$ are represented in the inertial frame, while the vehicle’s angular velocity ${\boldsymbol\Omega}$ is represented in the local frame.

Furthermore, $\mathbf{R_{il}}$ is the $ZYX$ rotation matrix from local to inertial coordinates, as shown in

where $c(\cdot)$ is the cosine and $s(\cdot)$ is the sine of the angle.

$\mathbf{S}(\boldsymbol\Psi)$ is the coupling matrix, expressing an interconnection in the system’s angle change, and it is given by

C Aerodynamic Effects

There are several aerodynamic and gyroscopic effects that can modify the force model introduced in (3)–(6), the most important of them being the ones called blade flapping and induced drag, which are of significant importance in understanding the natural stability of quadrotors. These effects are particularly relevant because they induce forces in the $x$-$y$ rotor plane of the quadrotor, precisely the under-actuated directions in the dynamics that cannot be easily dominated by high gain control [43].

While disturbances caused by drag forces are already included in the basic model in (4) as an additive term, the inclusion of blade flapping equations is not so straight-forward. This phenomenon is caused by the flapping effect generated in the blades of the propellers by horizontal translation of the rotor in the air. During the translation of the rotor, one blade (called the advancing blade) will move in the direction of translation, while the other (the retreating blade) will move opposite to the direction of motion. The net velocity of the advancing blade will increase, and that of the retreating blade will decrease with respect to the air. This creates a disparity of lift between the advancing side of the rotor and the retreating side of the rotor (see Fig. 3) [17].

Figure  3.  Perspective views of a rotor during flight. Longitudinal and lateral blade flapping angles are $a_{1,s}$ and $b_{1,s}$, respectively, and $\Psi$ here is the blade rotation azimuth angle.

DownLoad: Full-Size Img PowerPoint

According to Pounds et al. [17], a rotor blade balances the lift imbalance (aerodynamic force), centripetal force, and blade weight force by allowing free rotation around a hinge or flexure, and this is known as rotor or blade flapping. The advancing blade will rise, the retreating blade will fall, as in Fig. 3, and therefore the rotor plane will tilt in an effort to meet a new equilibrium. Increasing translational velocity leads to increasing rotor tilt. When the rotor plane tilts, the thrust force generated by the blades is correspondingly tilted, resulting in a horizontal component that opposes the direction of translation.

Allowing the rotor blade to bend is very important in the mechanical design of the quadcopter. A rotor that is too rigid can cause transmission of these aerodynamic forces directly through to the rotor hub and may result in a mechanical failure of the motor mounting or the airframe itself [43]. Leishman [44] also states that if the rotor disc is perfectly rigid, the helicopter will flip, as was the case in early experiments.

As aforementioned, rotor flapping causes the rotor plane to tilt and hence the thrust force will also be tilted and will have a horizontal component. The total thrust $\mathbf{T}^B$ that appears in (4) results from the sum of the thrusts $\mathbf{T}^B_i$ produced by each of the four propellers $i = 1, 2, 3, 4$ [10]

with $T_i = \alpha_T \omega_i^2$, where $\alpha_T$ is called the thrust coefficient, $b_{1,s,i}$ and $a_{1,s,i}$ are detailed in Table II, as described by Riether [10].

Table  II.  Aerodynamic Effects on Total Thrust $\mathbf{T}^B_i$

Relative air speed at propeller $i$: $\mathbf{v}^B_{ra,i}$	${\bf\Omega} \times \mathbf{r}^B_{prop,i} + \mathbf{R}_a^{-1}\mathbf{v}$
Planar components: $\mu_i$	$\sqrt{{v^B_{x,ra,i}}^2 + {v^B_{y,ra,i}}^2}/ \left | \omega_i\frac{d}{2} \right |$
Non-dimensionalized normal inflow: $l_i$	$|v^B_{{\textit{z}},ra,i}| / \left | \omega_i\frac{d}{2} \right |$
Sideslip azimuth relative to x-axis $x$: $j_i$	$\arctan\frac{v^B_{x,ra,i}}{v^B_{y,ra,i}}$
Sideslip rotation matrix: $\mathbf{J}_{ss,i}$	$\begin{bmatrix} \cos(j_i)&-\sin(j_i)\\ \sin(j_i)&\cos(j_i)\end{bmatrix}$
Longitudinal flapping: $\beta_{lf,i}$	$\mathbf{J}^T_{ss,i}\begin{bmatrix} ((\frac{8}{3}\theta_{b0} + 2\theta_{b1}) - 2l_i)/(1/\mu_i - \mu_i/2)\\0\end{bmatrix}$
$a_{1,s,i}$	$\beta_{x,lf,i} - 16Q/(\gamma_{aero}|\omega_i|)$
$b_{1,s,i}$	$\beta_{y,lf,i} - 16P/(\gamma_{aero}|\omega_i|)$

 | Show Table

DownLoad: CSV

In this table, $\mathbf{r}^B_{prop,i}$ refers to the propeller $i$’s position with respect to the center of mass, $\theta_{b0}$ is the equivalent blade pitch at the rotor axis, and $\theta_{b1}$ is the washout of a linear twist blade [17].

The constant $\gamma_{aero} = {\rho a_0 c R}/{I_b}$ is the non-dimensional Lock Number [19], [44], which can be viewed as a measure of the ratio of aerodynamic forces to inertial forces. Here, $I_b$ is the moment of inertia of the blade about the hinge, $c$ is the chord of the blade, $R = {d}/{2}$ is the rotor radius, $\rho$ is the density of air, and $a_0$ represents the slope of the lift curve per radian [19], [17]. For a typical helicopter rotor, the value of the Lock Number varies from $5$ to $10$ [44]. In the case of our minidrone, it is $0.6051$, as detailed in Table III.

Table  III.  Constants Shown in Table II

Propeller $i$’s positions: $\mathbf{r}^B_{prop,i}$	$\begin{bmatrix} \pm 0.0441\;\;\;\pm 0.0441\;\;\;-0.0159 \end{bmatrix}^T$
Rotor radius: $d/2$	$0.0330$
Lock number: $\gamma_{aero}$	$0.6051$
Blade root angle: $\theta_{b0}$	$0.2548$
Blade twist angle: $\theta_{b1}$	−0.1361

 | Show Table

DownLoad: CSV

Finally, the vector of torques ${\boldsymbol\tau} = [\tau_\phi \;\; \tau_\theta \;\; \tau_\psi]^T$ acting on the body-frame of the drone can be calculated as

where $Q_i = \alpha_Q \omega_i|\omega_i|$ is the propeller-drag-induced torque [10]. Note that, in (13), $\mathbf{T}^B_i$ is as in (12) and is dependent on aerodynamic effects, since all the terms multiplied by the Lock Number $\gamma_{aero}$ come from aerodynamics [44], as is the case with $b_{1,s,i}$ and $a_{1,s,i}$. Hence, the vector of torques ${\boldsymbol\tau}$ will likewise be influenced by the aerodynamic effects. Appendix presents a more thorough calculation of the thrust rotation caused by aerodynamic effects, which also provides an alternative manner to visualize the thrust equation.

These aerodynamic effects are neglected in most models, assuming $a_{1,s,i} = b_{1,s,i} = 0$ and ignoring the ratio of aerodynamic forces to inertial forces. Therefore, in those cases, the thrust vectors are supposed to be fully aligned with the body-frame z-axis and do not have horizontal components caused by blade flapping. The angles $a_{1,s,i}$ and $b_{1,s,i}$ cannot be zero, as explained above, after all the rotor bend is essential to the quadrotor design. Since increasing translational velocities increases rotor tilt, aggressive maneuvering is specifically subject to the consequences of ignoring rotor flapping, as well as small aircrafts, as referenced in Section I. In this paper, we aim to compensate for the described aerodynamic effects using control allocation, without the need to include them in the control algorithm itself.

Overlooking these effects results in $\mathbf{T}^B_i \approx [0\;\; 0\;\; T_i ]^T$ and

where $T \approx \sum T_i$. With these assumptions, the plant input is then simplified from six dimensions $[ \mathbf{T}_c^B \;\; \tau_{\phi_c} \;\; \tau_{\theta_c} \;\; \tau_{\psi_c}]^T$ to four dimensions $[ T_c \;\; \tau_{\phi_c} \;\; \tau_{\theta_c} \;\; \tau_{\psi_c} ]^T$. Hence the above mentioned equations can be simplified as

Definition 4 (Mixing Matrix Control Allocation (MMCA)): The control allocation performed by means of a matrix that provides an invertible relation between squared motor speeds and the plant input $[ T_c \;\; \tau_{\phi_c} \;\;\tau_{\theta_c} \;\; \tau_{\psi_c} ]^T$, as in (15), provided that the matrix is non-singular and the angular velocities $\omega_i$ are not saturated.

If we write (15) as

the angular velocities would then be determined as

where $\mathbf{B}$ is the aforementioned control effectiveness matrix, mixing matrix, or allocation matrix.

Remark 4: Assuming the angular velocities $\omega_i$ are not saturated is actually a very bold assumption when dealing with quadrotor systems [8]. Also, with (11) and (13), obtaining the analytical inversion $\widetilde{F}^\ast$ that transforms the virtual thrust and torques $[\mathbf{T}_c, {{\boldsymbol\tau}}_c]$ into the angular velocities $\omega_i$ would be much more complex than in (17). That is the reason why this simplification disregarding aerodynamic effects is prevailing in control models found in the state-of-the-art literature.

The new control allocation technique explained further along in this paper proposes precisely to perform the inversion of (11) and (13), by means of a function fitting neural network that models $\widetilde{F}^\ast$, to take into consideration the aerodynamic forces that influence the quadcopter behavior and therefore command the correct angular velocities to the actuators.

Additionally, expanding (13), we obtain the following:

which gives us, for each actuator $i$,

In (19)–(21), we make explicit the influence of the thrust components on $x$ and $y$-axes on the torques that perform the roll, pitch, and yaw movements, $\tau_\phi$, $\tau_\theta$, and $\tau_\psi$, respectively. These components are neglected when using the traditional MMCA, which considers that the thrust only has the z-axis component. We state once more that, due to aerodynamic effects such as blade flapping, however, thrust and moments are not necessarily orthogonal to the body axis [19], so these $x$-axis and $y$-axis components do exist. In Section V, we show that these components $T_x$ and $T_y$ are not of negligible order when compared to the order of the system torques they have an influence on.

D Quadrotor Control Strategy

The block diagram for the position and attitude control system is shown in Fig. 4. It is composed of an external PD trajectory control loop and an internal PID loop that controls the attitude which must be achieved so that the desired path is followed. Altitude or height control is also made through a PID controller.

Figure  4.  Block diagram for the quadrotor control strategy.

DownLoad: Full-Size Img PowerPoint

The Trajectory Generator block sends the path commands $(x,y)$ to the Position Controller. This controller outputs the roll ($\phi$) and pitch ($\theta$) commands to the Attitude Controller. The altitude (${\textit{z}}$) command comes directly from the Trajectory Generator, and the yaw ($\psi$) controller is independent. The respective adjustments in the angular velocities of the motors to achieve the desired behavior in $\phi$, $\theta$, $\psi$ and ${\textit{z}}$ will then be sent to the quadcopter through the Control Allocation. The sensing system must be adequate for good state feedback and, thereafter, an accurate calculation of the error.

We have the following expressions for the attitude PID controller, where $e_{\phi}$ and $e_{\theta}$ represent the roll and pitch errors, respectively; $P$ is the angular velocity around the $x$-axis in the body-frame, and $Q$ is the angular velocity around the $y$-axis in the body-frame:

As aforementioned, the attitude controller commands are generated by the PD trajectory controller, shown in (24) and (25). Note that we need to perform a rotation around the z-axis ($\psi$) during the error calculation, because the body and inertial frames are not necessarily aligned in the same direction.

The yaw control is made separately, as in (26), where $e_{\psi}$ represents the yaw error and $R$ is the angular velocity around the z-axis in the body-frame

Finally, the height control is performed through a PID, as shown below:

Since we aim to compare two different control allocation methods in this paper, the following assumption is naturally required.

Assumption 2: For the same controller design, simply replacing the traditional mixing matrix control allocation with the neural network control allocation here proposed does not harm the control project.

In order to verify this assumption, we designed a simple controller for the quadrotor system with the mixing matrix and tested to see if the project criteria were maintained when replacing it with the novel neural network control allocation we propose further ahead in this paper.

We set the project criteria of 1 s rising time ($t_r$), which is the time when the output first reaches the reference, and $3$ s settling time ($t_s$), defined as the time when the output is within a 1% error when compared to the steady state. With this, we specified the corresponding proportional ($K_p$), integral ($K_i$) and derivative ($K_d$) attitude controller gains, shown in Table IV.

Table  IV.  Step Response Characteristics

	Rise time ($t_r$)	Settling time ($t_s$)	Overshoot
Mixing Matrix	$0.605$ s	$3$ s	$14.9\%$
Neural Network	$0.6$ s	$3$ s	$14.6\%$
$K_p$		$0.0092$
$K_i$		$0.0097$
$K_d$		$0.003$

 | Show Table

DownLoad: CSV

The controllers were implemented in discrete-time, with the Bogacki-Shampine method, which is a requirement for it to be deployed to the hardware for real robot experiments afterwards. The step size was of $5$ ms, which was chosen to match the sampling time of the quadrotor during flight experiments. Other parameters necessary for the simulation, such as constants that model the aircraft, were obtained from the documentation provided by the manufacturer.

The obtained step response behaviors are presented in Fig. 5 and detailed in Table IV, for the traditional mixing matrix control allocation and the novel neural network control allocation. We can infer, comparing both step responses, that Assumption 2 is in fact true, as the difference between them is negligible. Therefore, it is possible to use the new control allocation approach and apply it to more elaborate maneuvers to observe the results. The PD-PID controllers were used here merely as a tool to perform the experiments and simulations, and are not the focus of the contributions of this paper.

Figure  5.  Roll step responses for the control system using the traditional mixing matrix (MMCA), on the top, and the proposed network control allocation (NNCA), on the bottom.

DownLoad: Full-Size Img PowerPoint

IV. Neural Network Control Allocation

Generally speaking, a control allocation solves an undetermined, typically constrained, system of equations [29]. The input to the control allocation block is the total control effect to be produced, or the virtual control input $\tau_c(t) \in \mathbb{R}^k$. Its output, on the other hand, is the achieved actuator command $u(t) \in \mathbb{R}^m$, where $m > k$ if the system is over-actuated, or $m < k$ if the system is under-actuated. That is, with a given $\tau_c(t)$, we search $u(t)$ such as

where $g : \mathbb{R}^m \rightarrow \mathbb{R}^k$ is the mapping from the achieved actuator commands to the virtual control inputs in the system to be controlled [6]. In the control allocation literature, (28) usually becomes

with B being called the control effectiveness matrix of dimensions $k \times m$ and rank $k$. Fig. 6 illustrates the control allocation process. It must transform the virtual inputs $\tau_c$ into commands $u$ sent to the actuators such that they generate the allocated inputs $\tau$ = $\tau_c$.

Figure  6.  Control system structure including control allocation. The motion control specifies a virtual control effort $\tau_c$ that must be produced, based on the reference signal r and the state feedback $x$. The control allocation distributes this control demand among the actuators, which generate the actual allocated control effort $\tau$. If the allocation is successful, then $\tau_c = \tau$

DownLoad: Full-Size Img PowerPoint

It is important to stress that, when it comes to the process of control allocation, the terminology used depends on the application in mind, so the control allocator can also be known as control mixer, control blender, control selector, or control distributor [6]. Therefore, the control effectiveness matrix can be referenced by several notations. In quadcopter applications, it is usually called mixing matrix, referencing the control mixer; the control inputs are the thrusts and torques, and the actuator commands are the angular velocities $\omega_i$ for each rotor.

The traditional formulation of the linear control allocation problem is the system (30), which consists of (28) with the application of the actuator constraints [5], [6].

There are three possible outcomes when solving (30):

1) There is an infinite number of solutions;

2) There is one unique solution;

3) No solution exists.

Many proposed control allocation methods correspond to different ways of computing the solution for a certain allocation objective, which most of the times is to avoid actuator saturation, rather than different objectives. The majority of well-established control allocation algorithms depend on numerical optimization to solve the control allocation problem. In this work, we briefly mention and describe the principal allocation methods according to Härkegård [6].

Optimization-based control allocation methods rely on a pragmatic interpretation of the control allocation problem, which states: for the system (30), determine a feasible control input so that this system is satisfied. If there are several solutions, pick the optimal one. In direct control allocation, the choice of control input is made unique by geometric reasoning. In daisy chain control allocation, the allocation suite is divided into groups, which are successively employed to generate the total control effort. All of these alternatives can be solved using numerical algorithms.

The main problem approached when using control allocation in the literature is the saturation of actuators in over-actuated systems. The principles of control allocation are however general and can be used to optimize other disturbances in the systems, from power or fuel consumption to aerodynamic effects [5], which is the case of this paper.

A Neural Network Design

The problem of allocating the proper angular velocities to the rotors, which will generate the commanded virtual thrust and torques, is solved with an analytical inversion of (11) and (13), which transforms the angular velocities into thrust and torques, $\omega_i \rightarrow [\mathbf{T}, {\boldsymbol\tau}]$.

We need therefore a solution to perform the opposite operation, $\widetilde{F}^\ast$ (as described in Section III-A), from the virtual control inputs to the rotor velocities, $[\mathbf{T}_c, {\boldsymbol\tau}_c] \rightarrow \omega_i$, considering the aerodynamic effects introduced in Section III. This inverse transformation would derail the invertible relation of the classic mixing matrix and hinder its usage. In this work, we propose the use of a neural network to perform the function fitting for finding $\widetilde{F}^\ast$, and use it as the control allocation.

Definition 5 (Neural Network Control Allocation (NNCA)): The control allocation performed by means of a function fitting neural network that provides a transformation $\widetilde{F}^\ast$ of the plant input $[\mathbf{T}_c, {{\boldsymbol\tau}}_c] = [T_{{\textit{z}}_c} \;\; \tau_{\phi_c} \;\; \tau_{\theta_c} \;\; \tau_{\psi_c} ]^T$ to the actuator commands $\mathbf{u} = [ \omega_1 \;\; \omega_2 \;\; \omega_3 \;\; \omega_4 ]^T$, taking into account aerodynamic effects as described in (11) and (13).

That is, as in Section III-A,

Function fitting is the process of training a neural network on a set of inputs with the aim of producing an associated set of target outputs. After the network has been created with the desired hidden layers and the training algorithm, the neural network must be trained using a set of training data. Once the neural network has fitted the data, it forms a generalization of the input-output relationship. The trained network can therefore be used to generate outputs for a new set of inputs that was not used in training. Neural networks are good at fitting functions. There is proof that a fairly simple neural network can fit any practical function [45], [46].

The proposed function fitting neural network is composed of a two-layer feed-forward network with sigmoid hidden neurons and linear output neurons. This structure was chosen since it can fit multi-dimensional mapping problems arbitrarily well, provided consistent data and enough neurons in its hidden layer are given. In this case, 20 neurons were used for the first step of the training, and 10 neurons for the second one, as described below.

The neural network control allocation approach has a two-step training process in order to take into account the aerodynamic effects, as well as achieving desirable training performance and mean squared error (MSE) order. The second step of the training process is independent from the mixing matrix.

In the mathematical model presented in Section III-C, we do not represent uncertainties caused by nonlinear terms. Nonetheless, the curve fitting provided by the neural network imposes an approximation error influenced by uncertain nonlinear terms.

With the purpose of minimizing the effect of these uncertainties, we aimed for a certain MSE order based on recent works in the literature regarding curve-fitting using neural networks [47]-[49]. To achieve this goal, we used our two-step training process, as will be described next. The references found in literature used radial basis function (RBF) networks for the curve-fitting. Considering that, we first trained RBF neural networks using the same mentioned two-step training method. Our objective was to compare the MSE of these RBF results to the MSE of the network structure we are using and see if they agree with the results found in the literature [47]-[49]. The complete neural network training is described next and summarized in Table V.

Table  V.  MSE for the Neural Network Training Process

		1st training: Bayesian-regularization ($20$ neurons, $64\;004$ samples)	2nd training: Levenberg-Marquadt ($10$ neurons, $34\;136$ samples)
Sample division	Dataset	Mean-squared error	Mean-squared error
$70\%$	Train	$3.28812 \times 10^{-2}$	$4.4857 \times 10^{-4}$
$15\%$	Validation	–	$5.98463 \times 10^{-4}$
$15\%$	Test	$4.09559 \times 10^{-2}$	$4.50975\times 10^{-4}$

 | Show Table

DownLoad: CSV

For the RBF networks, the training process was the same as below, with the exception of the training algorithm, which was performed using the MATLAB function newrbe(). In the first training, the RBF network obtained an MSE of $\mathbf{4.2 \times 10^{-2}}$. In the second training, the RBF network achieved an MSE of $\mathbf{6.0297 \times 10^{-4}}$. These error orders are expected, according to the literature [47]-[49]. If we compare these results with the ones of our proposed network, detailed in Table V, it is possible to conclude that not only are they within the expected MSE order of previous works, but our results are slightly better than the MSE obtained with the RBF networks and the same training data.

● First Training

1) An open-loop simulation of the quadrotor system without any controller was performed, using pseudo-random binary sequence (PRBS) pulses with amplitude $0.2$ as mechanical conjugates ($\mathbf{T}_c, {\boldsymbol\tau}_c$) to collect the corresponding angular velocities outputs ($\omega_i$). Other variables ($\mathbf{v}$ and ${\boldsymbol\Omega}$) were also collected for training, as detailed in Fig. 7.

Figure  7.  Block diagrams representing the variables used for neural network training and therefore for neural network control allocation.

DownLoad: Full-Size Img PowerPoint

2) The collected angular velocities data were applied to the aerodynamic effects equations (11) and (13) aiming to collect the real mechanical conjugates ($\mathbf{T}, {\boldsymbol\tau}$) established at the aircraft.

3) The input-output data ($\mathbf{T}, {\boldsymbol\tau}, {\boldsymbol\Omega}, \mathbf{v} \rightarrow \omega_1, \omega_2, \omega_3, \omega_4$) were then used for the first training step, using the Bayesian-regularization algorithm.

With the first training, we obtained a mean squared error of the order $10^{-2}$, as detailed in Table V. To refine the network and obtain an MSE order of at least $10^{-3}$ as in the mentioned recent works found in the literature [47]-[49], the second training was performed, using the results of the first one.

● Second Training

1) A new set of input-output data were collected from a closed-loop simulation with the first trained network as the control allocation method and a PD-PID error controller.

2) 3-2-1-1 maneuvers were used as attitude commands, and the network obtained in the first step converted the virtual control commands ($\mathbf{T}_c, {\boldsymbol\tau}_c$) into angular velocities ($\omega_i$).

3) Once more, the collected angular velocities data were applied to the aerodynamic effects equations (11) and (13) aiming to collect the real mechanical conjugates ($\mathbf{T}, {\boldsymbol\tau}$) established at the aircraft.

4) The input-output data ($\mathbf{T}, {\boldsymbol\tau}, {\boldsymbol\Omega}, \mathbf{v} \rightarrow \omega_1, \omega_2, \omega_3, \omega_4$) were then used for the second neural network training, performed using the Levenberg-Marquadt algorithm. This time, we obtained the MSE of the order $10^{-4}$.

Both training steps were carried out for a maximum of 1000 epochs or until the validation error failed to decrease for 6 iterations, and we obtained Regression values of $R = 0.99999$ for the two steps. Regression values measure the correlation between outputs and targets. An R value near 1 means that the variables are strongly related, while an R value near $0$ means that they have little to no impact on each other.

The whole training process was performed offline, meaning that during the flight this allocation method requires low computational power, because the network is already trained and no update is made on-line. The classic mixing matrix was only used in the first step of the training to collect data.

While the mixing matrix control allocation uses only the thrust and torques to perform the matrix inversion and obtain the squared angular velocities, more variables were added to the neural network training process for the purpose of obtaining a better fit, according to the aerodynamic effects modeled in Section III-C. See Fig. 7, which details the input and output variables used in training. Therefore we used 10 input elements ($T_{\textit{z}}, \tau_{\phi}, \tau_{\theta}, \tau_{\psi}, \Omega_{\phi}, \Omega_{\theta}, \Omega_{\psi}, v_x, v_y, v_{\textit{z}}$), with $T_{\textit{z}}$ as the z-axis component of the trust $T$, which comes from the height controller, to obtain 4 output elements ($\omega_1, \omega_2, \omega_3, \omega_4$). During the UAV flight, these additional state variables ($\mathbf{v}$ and ${\boldsymbol\Omega}$) are obtained from sensing.

V. Results

The platform used in the flight experiments was a PARROT Mambo minidrone, illustrated in Fig. 1, which contains a downward-facing camera for positioning using optical-flow odometry, an ultrasonic sensor, a barometric sensor, and an inertial measurement unit (IMU). All these sensors are used for stabilization purposes. The minidrone dimensions are $180$ mm × $180$ mm × $40$ mm, with a weight of $70$ g. Power is provided by a $3.7$ V, $550$ mAh, 2 Wh Li-Po battery, enabling a flight time of 7–10 minutes, depending on the nature of the task to be performed and the accessories implemented. The applied control strategy for the quadcopter operation consisted of the PD-PID controller mentioned in Section III-D, again implemented in discrete-time, with the Bogacki-Shampine method and a step size of $5$ ms. In Fig. 4, the actuators (rotors) and quadrotor dynamics blocks represent the minidrone, while the other blocks correspond to the algorithm embedded in the aircraft.

We evaluated the proposed method by comparing the mechanical conjugates at the output of the controller (the virtual control commands $\mathbf{T}_c, {{\boldsymbol\tau}}_c$) to those established at the aircraft ($\mathbf{T}, {\boldsymbol\tau}$), which were calculated using the inflow and blade flapping equations. We also analyzed changes in the controller performance. The neural network was able to closely match $\mathbf{T}, {\boldsymbol\tau}$ and $\mathbf{T}_c, {{\boldsymbol\tau}}_c$, while the classic control allocation with the mixing matrix could not achieve the same performance. The error between the two mechanical conjugates (${\boldsymbol\tau} - {\boldsymbol\tau}_c$ and $\mathbf{T} - \mathbf{T}_c$), called allocation error, was also determined in both cases. We were able to enhance the performance of the PD-PID controller with the new control allocation, without altering the controller gains, that is, the tuning was the same for both the traditional mixing matrix and our proposed neural network control allocations.

Our position reference consisted of a set of waypoints for the quadrotor to follow. They formed a square with a $1.5$ m side and a transition radius of $0.3$ m to maneuver through the corners, as well as take-off and landing commands. The experiment consisted of the minidrone taking-off in the bottom-left corner of the square. Then the UAV performs a pitch forward movement until it reaches the next corner, then a roll lateral movement, a pitch backward movement, and then again a roll lateral movement. This path is repeated once more, and then the minidrone lands, resulting in a total of 9 waypoints [0 0 1; 1.5 0 1; 1.5 1.5 1; 0 1.5 1; 0 0 1; 1.5 0 1; 1.5 1.5 1; 0 1.5 1; 0 0 –1]. Here the z-axis value 1 represents the height (1 m) at which the drone navigates along the square path. Fig. 8 shows the trajectory reference and output for the complete path of the waypoint following, with take-off and landing, highlighting the aircraft positioning in black, for clarity’s sake. The yaw orientation reference of the quadrotor was fixed throughout the experiment.

Figure  8.  Trajectory reference and output for both loops of the waypoint following, with take-off and landing, highlighting the aircraft positioning, for the sake of clarity.

DownLoad: Full-Size Img PowerPoint

The test scenario was planned for a minidrone, which is a type of aircraft especially susceptible to aerodynamic effects influence, as aforementioned. The path was designed to induce mildly abrupt roll and pitch maneuvers (the corners of the square, that is, the waypoints) that would also allow an adequate observation of said effects. The behavior of the quadrotor throughout the trajectory however was not intended to be aggressive as mentioned in previous works in the literature, because we wanted to highlight the importance of the consideration of aerodynamic effects in the control allocation even without aggressive maneuvers. Note the trajectory used in this experiment differs from the maneuvers used for the neural network training explained in Section IV-A, hence the training does not need to be repeated for each specific maneuver. Furthermore, Table III shows some constants relative to the quadcopter in use, which were obtained from the documentation made available by the manufacturer. These constants are used in the aerodynamic effects modeling in Table II.

Fig. 8 shows the three-dimensional representation of the commanded and followed trajectory, with take-off and landing, for the mixing matrix and neural network control allocations. Figs. 9 and 10 show the trajectory projection in the $xy$-plane of the first loop for the same respective cases. We omit the aircraft’s second loop through the square in Figs. 9 and 10 to facilitate visualization. We set a transition radius as $0.3\;{\rm{m }}$, shown as the red circle in the figures. It means that when the drone is approaching a waypoint, this value is used as the radius around the waypoint where the drone starts moving towards the next point. If the transition radius is set as $0.3\;{\rm{m}}$, and the drone is $0.3\;{\rm{m}}$ away from the waypoint, the drone starts moving towards the next waypoint. The accuracy of reaching a specified waypoint increases with a decrease in the value of transition radius. Notice that the path following is more smoothly achieved with the neural network control allocation approach, while with the mixing matrix some stuttering is present when pursuing the commanded trajectory. This is emphasized in Fig. 11, where some parts of the trajectory are highlighted to show the difference in path following for the mixing matrix and neural network control allocations.

Figure  9.  Projection $(x,y)$ of trajectory reference and output when using the mixing matrix control allocation.

DownLoad: Full-Size Img PowerPoint

Figure  10.  Projection $(x,y)$ of trajectory reference and output when using the neural network control allocation.

DownLoad: Full-Size Img PowerPoint

Figure  11.  Augmented parts of the trajectory to show the difference in path following behavior for the mixing matrix (MMCA) and neural network (NNCA) control allocations.

DownLoad: Full-Size Img PowerPoint

The respective mean squared error (MSE) calculations for each axis are shown in Table VI, for angular and linear positions. For every case except the yaw angle and position z, the MSE is smaller when using the neural network, especially in the case of the roll angle, where we had the most improvement when using the new control allocation.

Table  VI.  Mean Squared Errors for the Angular and Linear Positions for Both Mixing Matrix and Neural Network Control Allocations

	MSE Roll	MSE Pitch	MSE Yaw
Mixing Matrix	$7.3017 \times 10^{-4}$	$4.7728 \times 10^{-4}$	$3.2070 \times 10^{-6}$
Neural Network	$1.4444 \times 10^{-4}$	$4.2851 \times 10^{-4}$	$7.2123 \times 10^{-5}$
	MSE $x$	MSE $y$	MSE z
Mixing Matrix	$0.0372$	$0.0347$	$0.0163$
Neural Network	$0.0300$	$0.0300$	$0.0176$

 | Show Table

DownLoad: CSV

The performance of the UAV during the roll angle maneuver can be seen in Fig. 12, for the mixing matrix, and in Fig. 13, for the neural network control allocation. The mixing matrix approach shows a much noisier behavior and a higher amplitude variation when compared to the neural network control allocation approach, where the adherence of the output to the commands is much more reliable, stable, and therefore preferable for the quadcopter flight. The difference between the reference signals in both cases is due to the consideration of the aerodynamic effects in the control allocation algorithm.

Figure  12.  Roll commands and output when using the mixing matrix control allocation

DownLoad: Full-Size Img PowerPoint

Figure  13.  Roll commands and output when using the neural network control allocation

DownLoad: Full-Size Img PowerPoint

Figs. 14 and 15 show the comparison between the commanded virtual roll torques, which are commanded at the output of the controller. These figures also show the achieved roll torques at the aircraft, which were calculated using (13) with experimental flight data. Fig. 14 refers to the experiment with the mixing matrix, while Fig. 15 uses the neural network control allocation. The amplitude of the roll torques when using the neural network is smaller, requiring less effort from the quadrotor. Also, a comparison between the allocation errors for both cases is shown in Fig. 16, proving that the neural network allocation considering aerodynamic effects significantly reduces the allocation error.

Figure  14.  Comparison between commanded roll torques at the control output and those achieved at the plant for the mixing matrix control allocation.

DownLoad: Full-Size Img PowerPoint

Figure  15.  Comparison between commanded roll torques at the control output and those achieved at the plant for the neural network control allocation.

DownLoad: Full-Size Img PowerPoint

Figure  16.  Comparison between the allocation errors for the achieved roll torques ($\tau_\phi$) at the plant for both the mixing matrix (in blue) and neural network (in red) control allocations.

DownLoad: Full-Size Img PowerPoint

Likewise, Figs. 17 and 18 present the corresponding comparison for commanded virtual thrusts and achieved thrusts, for mixing matrix and neural network control allocations, respectively. We note a similar behavior as the roll torque case mentioned above. There are less noisy oscillations when using the neural network, as well as smaller amplitudes. Fig. 19 shows the allocation errors for the two methods, where the error when applying the neural network control allocation is several times smaller (in the order of $10^{-1}$) than the traditional alternative with the mixing matrix.

Figure  17.  Comparison between the commanded thrust at the control output and that achieved at the plant for the mixing matrix control allocation.

DownLoad: Full-Size Img PowerPoint

Figure  18.  Comparison between the commanded thrust at the control output and that achieved at the plant for the neural network control allocation.

DownLoad: Full-Size Img PowerPoint

Figure  19.  Comparison between the allocation errors for thrusts, in the z-axis ($T_{\textit{z}}$), achieved at the plant for both the mixing matrix (in blue) and neural network (in red) control allocations.

DownLoad: Full-Size Img PowerPoint

Finally, Figs. 20 and 21 show the thrust $\mathbf{T}^B$ components on the $x$ and $y$ axes, $T_x$ and $T_y$, during the quadrotor flight when using the mixing matrix, which considers these components to be zero or negligible. Not only these components do exist, but they are not of negligible order when compared to the order of the system torques they have an influence on. For instance, Fig. 21 shows $T_y$ which has influence on the roll torque, as in (19). Note that the behavior of the roll torques seen in Fig. 14 greatly resembles that of $T_y$ in Fig. 21.

Figure  20.  Thrust component in the x-axis, $T_x$, when using the mixing matrix control allocation.

DownLoad: Full-Size Img PowerPoint

Figure  21.  Thrust component in the y-axis, $T_y$, when using the mixing matrix control allocation.

DownLoad: Full-Size Img PowerPoint

Measures were also taken during the flight experiments to minimize the error due to nonlinear uncertainties, since they were performed in a controlled indoors environment. The uncertain nonlinear terms were not directly treated in the modeling presented in this paper, as it was not the focus of this work. However, the practical results discussed here already show the advantages of the proposed method. In a future work, an approach considering these uncertainties could be included.

VI. Conclusions

In this paper, we demonstrated that the disregard of aerodynamic effects in the quadrotor control intending to simplify the process of converting the conjugates at the controller output to angular velocities to be commanded to the rotors greatly affected the flight performance, most specifically thrust and torque behaviors. We proposed that the aerodynamic effects could be taken into account by applying a more refined control allocation process that used a function fitting neural network to transform virtual control inputs into actuator commands. The control algorithm, which was a PD-PID controller, did not need to acknowledge the aerodynamic equations directly.

The network training was performed offline, requiring very low computational power when compared to other optimization-based control allocation algorithms. The neural network required a two-step training process that was independent of the classic simplified allocation method on its second step. The first training was performed in an open-loop system, via PRBS pulses. A second training was executed with the purpose of refining the network in an equivalent closed-loop system, and in this case, the allocation was implemented using the open-loop trained neural network alongside the PD-PID error controller.

The experimental flight results showed that the aircraft behavior was smoother when using the neural network control allocation approach, the allocation error was smaller and the signals were less noisy, that is, the consideration of aerodynamic effects into the system interfered directly and positively with the quadrotor performance. Finally, the performance of the controller also improved with the use of the neural network control allocation, as well as the quality of the thrust and torque signals.

The focus of this paper was the control allocation technique, not the controller project itself, which will be the focus of future works that will consider the aerodynamic effects in a nonlinear dynamic inversion controller using neural networks.

Appendix. Thrust Rotation

Assumption 3: The thrust $\mathbf{T}^B_i$ produced by each of the four propellers $i = 1, 2, 3, 4$ is equal to the vector $\alpha_T \omega_i^2 [0 \;\; 0 \;\; -1]^T$ rotated around the $y$ and $x$ axes, respectively, by the longitudinal and lateral flapping angles $a_{1,s,i}$ and $b_{1,s,i}$, as shown below:

where

Note that, with this development, we can actually rewrite (12), shown in [10], using (35), as a function of the rotation of the flapping angles, as in (32) shown in Proposition 3.

Lemma 1: The longitudinal and lateral flapping angles $a_{1,s,i}$ and $b_{1,s,i}$ are

and

where $\mathbf{e}_{\textit{z}}= [ 0\;\; 0 \;\; 1]^T$, $\mathbf{e}_y = [0 \;\; 1 \;\; 0]^T$ and $\mathbf{e}_x = [ 1 \;\; 0 \;\; 0]^T$ $\in \mathbb{R}^3$; $\mathbf{P}_{xy}$ and $\mathbf{P}_{\textit{z}}$ represent the projections on the $xy$-plane and z-axis respectively; and $\Theta = \frac{8}{3}\theta_{b0} + 2\theta_{b1}$.

Proof: Now, analyzing the other terms on Table II, let us start with the sideslip rotation matrix $\mathbf{J}_{ss,i}$. It can be interpreted as the rotation of $j_i$, defined in Table II and in Fig. 22, along the z axis, as in (38) below:

Figure  22.  $\mathbf{v}_{ra,i}$ projection on the $xy$-plane.

DownLoad: Full-Size Img PowerPoint

It is important to remember that (39) and (40) are equivalent because the inverse of a rotation matrix is equal to its transpose. Therefore, we can rewrite the longitudinal flapping $\beta_{lf,i}$, also on Table II, as

Now for the planar components $\mu_i$ and the non-dimensionalized normal inflow $l_i$, from Table II, we know that

therefore,

Also,

Hence, if $\mathbf{P}_{xy}$ and $\mathbf{P}_{\textit{z}}$ represent the projections on the $xy$-plane and z-axis respectively, we can write

and, from (43) and (44),

where $r = {d}/{2}$.

Now, taking $\cos(j_i))$ and $\sin(j_i))$ as in Fig. 22

we can rewrite the rotation matrix $R_{\textit{z}}(-j_i)$, shown in (40), as

If we consider $\mathbf{e}_{\textit{z}} = [ 0 \;\; 0 \;\; 1]^T$, $\mathbf{e}_y = [ 0 \;\; 1 \;\; 0]^T$ and $\mathbf{e}_x = [ 1 \;\; 0 \;\; 0]^T$ $\in \mathbb{R}^3$, we have

Using $\mathbf{P}_{xy}$, $\mathbf{P}_{\textit{z}}$ and (43) and (44), we have

and $a_{1,s,i}$ and $b_{1,s,i}$ can be expressed as

■

Corollary 1: Finally, from Assumption 3 and Lemma 1, it is possible to write an expression for $\mathbf{T}^B_i$ as

where $\Theta = \frac{8}{3}\theta_{b0} + 2\theta_{b1}$.

Therefore, the thrust $\mathbf{T}^B_i$ produced by the propeller $i$ is equal to the thrust $\alpha_T \omega_i^2 [0 \;\; 0 \;\; -1]^T$ rotated around the $y$ and $x$ axes, because of the longitudinal and lateral flapping angles. Such rotations not only depend on the rotational and translational states ${\boldsymbol\Omega}$ and $\mathbf{v}$ but also depend on the angular velocities of the propellers $\omega_i$. Thus, both the intensity and the rotation of the thrust $\mathbf{T}^B_i$ are functions of $\omega_i$.

Acknowledgment

The authors thank PPgEE, DEE-UFCG, CAPES and CNPq for the award of research grants and study fellowship during these investigations