Mathematical Modeling of the Coaxial Quadrotor Dynamics for Its Attitude and Altitude Control

Abstract

In this paper, an easily implementable coaxial quadrotor model and its validation on data from a real unmanned aerial vehicle (UAV), are presented. The proposed mathematical model consists of two parts: description of orientation and position of the UAV in the three-dimensional space. It takes into consideration the gyroscopic effect, influence of the Coriolis force, viscous friction and a several drag-like effects (blade flapping, rotor drag, translational drag and profile drag). In contrast to multirotor models available in the literature, this one is characterized by complementarity in relation to the available control techniques. Depending on selection of these techniques, the model can be narrowed (simplified) to meet the needs without the loss of behaviour adequacy to a real UAV.

1. Introduction

Coaxial quadrotor (coaxial quadcopter, octoquad) is usually defined as a multidimensional control plant, strongly non-linear, highly dynamic (because of the use of eight independent drive units located in pairs on the common cross-frame) and in general with parameters which are non-stationary. In some works, it is called octarotor, though it may be confusing since octarotors have in general a planar mechanical structure. Literature studies in the field of modeling coaxial quadrotors allow one to extract a number of key concepts. The first of them is accompanied by practical (constructional) terms, where the representation of a mathematical model of unmanned aerial vehicle (UAV) is intended to reflect the behaviour of a real quadrotor at the maximum simplification of its physical structure [1], neglect of certain aspects of aerodynamics—for instance: viscous friction [2] or drag-like effects [3]. In many cases, even basic laws of physics are sufficient for control purposes.

On the other side of available modeling concepts appears an approach, where one wishes to obtain a complete model—with possibly complete description of all physical laws that affect the quadrotor [4]. Thus object can be described using a wide range of mathematical techniques and benefits (e.g., quaternion notation [5]). This approach is justified and practically necessary, since it permits in a major way in simulation conditions to check aspects difficult to verify on the real UAV (e.g., aggressive maneuvers [6]) and those that are expensive. It may be noted that in order to seek for a “complete model” one can fall into the endless search for the sources and effects of physical laws that may still be relevant and have not been considered so far in already existing models of multirotor UAVs and coaxial units (analysis of the aerodynamic complexities [7], ground effect issues [8], hover/low-speed bare airframe dynamics [9], blade element momentum theory selected issues [10]). The passage of time provides new aspects for modeling (hitherto unknown phenomena, analyzed incorrectly or insufficiently precise). The technology development (measurements in wind tunnels [11], the use of Finite Element Method [12], new vision systems [13], parallelization of calculations, etc.) and software development [14] open new doors and opportunities for further research. Unfortunately, very complex models generate a high demand for computing power, which greatly limits their application in terms of time corresponding to a real coaxial quadrotor flight.

There are naturally also more consensual ways in modeling than the relationship: simple model vs. complex, and developed mathematical models are characterized more by their focused use (e.g., estimation and rejection of dynamic disturbances [15,16], payload improvement [17], multi-layered drift-stabilization [18]). A distinct center of gravity is placed for example in models developed for the purposes of control algorithms design (e.g., cascaded PD control [19], gain scheduling-based PID control [20], adaptive [21] and hierarchical integrating backstepping control [22], fault tolerant control [23], self-tuning fuzzy PID control [24], LQR control [25], multi-model control using fuzzy techniques [26]). Their complexity is lower than of the models used for simulation analysis only. Here the aim is to provide the highest control quality in a coaxial quadrotor, while modeling is limited only to the aspects needed for a potentially applied control algorithm [27]. Position and attitude control approaches from the last 10 years of UAV research may be found in [28].

Model simplification entails certain consequences, which however can be minimized through the use of estimation algorithms [29] to estimate neglected or unmeasurable values and state variables. Therefore, in accordance with the authors’ proposal, it is advantageous to adopt reasoning based on the search for balance between modeling, estimation and the used type of control. To make it easier, a new and relatively complex model is proposed. It is predefined in the designing of control algorithms and can be simplified or expanded with a more comprehensive analysis of the dynamics and electrical aspects [30].

To articulate contribution of the paper, to the authors’ best knowledge, this is the first journal paper about coaxial quadrotor UAVs that summarizes all necessary details of model (with the proposed important physical aspects) for control purposes in one place. The work presents a step-by-step method for modeling a quadrotor equipped with rotors whose number and location are nonstandard (rotors are in coaxial configuration and displaced in relation to the UAV’s center of gravity). From the idea, through the physical effects to mathematical conventions and formulas, the authors explain the model, tune it and validate it on real data from a coaxial quadrotor. This dynamical model takes into consideration the gyroscopic effect, influence of the Coriolis force, viscous friction and a several drag-like effects (blade flapping, rotor drag, translational drag and profile drag). In contrast to multirotor models available in the literature, this one is characterized by complementarity in relation to the available control techniques. Depending on selection of these techniques, the model can be narrowed (simplified) to meet the needs without the loss of behaviour adequacy to a real UAV.

The article is organized as follows: Section 2 presents principles of UAV control and the ways of defining quadrotor position and orientation in the 3D space. Section 3 contains the proposed coaxial quadrotor model. A brief description of the physical phenomena affecting the UAV, is given. Section 4 presents validation tests with the use of a proposed model of the Falcon V5 robot. Conclusions are provided in Section 5.

2. Coaxial Quadrotor Control Basics

2.1. Physical Construction of Coaxial Quadrotor

Each coaxial quadrotor is a multirotor UAV equipped with eight propellers. They are placed coaxially on the robot’s frame [31] (see Figure 1). This results in an increased payload at the similar mechanical structure compared to classical quadrotor. Thrust and torque are generated by the rotational movement and appropriate geometrical construction of propellers mounted directly on the Brushless Direct Current (BLDC) motor.

The propeller rotates in one direction, which is due to a constant pitch of the blade. Changes of voltages transmitted to the BLDC motor, i.e., Pulse-Width Modulation (PWM) duty cycles from the Electronic Speed Controller (ESC), result in corresponding changes of the rotational speed of the propeller.

Avionics of the UAV is based first and foremost on measurements from the available inertial measurement unit (IMU) containing three gyroscopes and triaxial accelerometer and magnetometer, which are used for the stabilization of the UAV’s position and orientation. More precisely, gyroscopes are used to determine the orientation relative to the base, and accelerometers with magnetometers are used to reset the gyro drift. Gyroscopes provide measurements of the angular velocity in the local coordinate system of the quadrotor, while accelerometers and magnetometers provide measurements of linear acceleration and magnetic induction, with respect to the global coordinate system ($\mathcal{EF}$—Earth Frame), projected onto the axes of the local coordinate system ($\mathcal{BF}$—Body Frame).

2.2. Change of the UAV’s Position and Orientation

The coaxial quadrotor can be considered as a flying robot with six degrees of freedom (DOFs). Three degrees are used to describe position, and the next three, orientation (via Euler or so-called Tait–Bryan angles). These angles allow the mathematical description to be converted from $\mathcal{BF}$ to the $\mathcal{EF}$ coordinate system. In the modeling process presented in the next chapters, the North–East–Down (NED) configuration of the reference system is used (Figure 2).

Every change of the UAV’s position $x,y,z$ expressed in relation to the beginning of the global coordinate system $X,Y,Z$, and every change of its orientation (rotations about Roll, Pitch and Yaw angles) measured in the local coordinate system and implemented in relation to the axis of the global coordinate system from Figure 2, is possible only by changing the rotational speed of appropriate propulsion units (and thus by the change of thrust and torque). Figure 3 shows how this is accomplished. From the perspective of control, the preferred solution is the architecture from Figure 4 which is based on a cascade control loop realization of movement and orientation of the UAV.

Thanks to this solution (cascade control structure), the problem of control of a plant with six DOFs becomes simplified to the problem of X, Y, Z position control with automatic conversion of proper orientation and adjustment of rotational speed that will ensure this. Six controllers are used for this purpose.

Note, that an interesting approach to the problem of UAV attitude control in cascade control structure, when there are discrepancies in the characteristics of the 4 actuators (electrical motors and propellers) with the use of cascade control strategy with a PD controller in the inner loop (to achieve stabilization), and PI controller in the outer loop (to ensure disturbance rejection), one may find in [32]. In [33], a model of constrained control scheme to steer an UAV to the desired position while ensuring constraints satisfaction at all times (based on the Explicit Reference Governor framework), is presented.

In order to facilitate the tuning of real UAV controllers, mathematical modeling based on physical phenomena affecting the coaxial quadrotor and physical modeling of the robot itself is used. In following sections key aspects necessary to build the models of coaxial quadrotor are shown.

3. The Coaxial Quadrotor Model

The mathematical model of coaxial quadrotor proposed by the authors describes the orientation and position of a robot in the 3D space. The basic formulas are as follows:

$$\begin{array}{c}\hfill \dot{\omega }={\mathbf{I}}^{-1}\left(\tau -\omega \times \left(\mathbf{I}\omega \right)-{\tau }_{\mathrm{gyro}}-{\tau }_d-{\tau }_t\right),\end{array}$$

(1)

$$\begin{array}{c}\hfill \dot{\mathbf{v}}=\frac{1}{m}\left(\mathbf{f}-{\mathbf{f}}_g-{\mathbf{f}}_{\mathrm{co}}-{\mathbf{f}}_d-{\mathbf{f}}_t\right),\end{array}$$

(2)

where left sides of (1) and (2) are angular and linear accelerations respectively. The components of (1) are input torque (from propellers), gyroscope torque due to the robot’s frame and propulsion units system, as well as aerodynamical effects and friction (${\tau }_d$ and ${\tau }_t$). In case of Equation (2), the elements are similar to the torques (those with the same subscripts), extended with projected gravitational force and Coriolis force (${\mathbf{f}}_g$ and ${\mathbf{f}}_{co}$, respectively). The presented model is explained in the further part of the section. It is divided into two basic parts: the position model (relates to positions, linear velocities and accelerations in reference and local Cartesian coordinate systems) and the orientation model (describes angular velocities and angles in X, Y and Z axis).

3.1. Euler Angles

The Formulas (1) and (2) describe the relation between the quadrotor state and input forces/torques. These equations are expressed in the local coordinate system. It means that all state variables are directly linked with the UAV. However, in case of control tasks, inertial coordinate systems are used. Thus, one needs to transfer the quadrotor state to the reference system chosen based upon the purposes. Basically, two methods are used: quaternions and Euler angles (or Tait–Bryan angles). A rotation matrix transferring an arbitrary vector from $\mathcal{BF}$ to $\mathcal{EF}$ is needed [34]. In this work, the authors have chosen the flight notation (“3–2–1”), which assumes that the final rotation matrix ${\mathbf{R}}_e^b$ is a result of multiplication of three basic rotation matrices, in the following order:

$$\begin{array}{c}\hfill {\mathbf{R}}_e^b={\mathbf{R}}_x\left(\varphi \right){\mathbf{R}}_y\left(\theta \right){\mathbf{R}}_z\left(\psi \right).\end{array}$$

(3)

Note that the basic rotation matrix in a given axis transfers a vector from an old coordinate system $\mathcal{SF}$ to a new one, ${\mathcal{SF}}^{\prime }$ (see Remark no. 1 in Appendix A):

$$\begin{array}{c}\hfill {\mathbf{v}}^{\left\{S{F}^{\prime }\right\}}={\mathbf{R}}_{x\left|y\right|z}{\mathbf{v}}^{\left\{SF\right\}}.\end{array}$$

(4)

This rule enforces the following rotations matrices forms:

$$\begin{array}{c}\hfill {\mathbf{R}}_x\left(\varphi \right)=\left[\begin{array}{ccc}1& 0& 0\\ 0& cos\varphi & sin\varphi \\ 0& -sin\varphi & cos\varphi \end{array}\right],\end{array}$$

(5)

$$\begin{array}{c}\hfill {\mathbf{R}}_y\left(\theta \right)=\left[\begin{array}{ccc}cos\theta & 0& -sin\theta \\ 0& 1& 0\\ sin\theta & 0& cos\theta \end{array}\right],\end{array}$$

(6)

$$\begin{array}{c}\hfill {\mathbf{R}}_z\left(\psi \right)=\left[\begin{array}{ccc}cos\psi & sin\psi & 0\\ -sin\psi & cos\psi & 0\\ 0& 0& 1\end{array}\right].\end{array}$$

(7)

The “3–2–1” assumes that during the rotations, the following steps are performed (Figure 5):

• Rotation Yaw about angle $\psi$ in the z axis of the $\mathcal{EF}$ (global coordinate system),
• Rotation Pitch about angle $\theta$ in the y axis of the ${\mathcal{EF}}^{\prime }$ (new coordinate system),
• Rotation Roll about angle $\varphi$ in the x axis of the ${\mathcal{EF}}^{″}$ (new local coordinate system), yielding $\mathcal{BF}$ as a product.

Under this assumption, the transformation matrix given in (3) may be written as:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & {\mathbf{R}}_e^b={\mathbf{R}}_x\left(\varphi \right){\mathbf{R}}_y\left(\theta \right){\mathbf{R}}_z\left(\psi \right)=\hfill \\ \hfill & =\left[\begin{array}{ccc}c\theta c\psi & c\theta s\psi & -s\theta \\ s\varphi s\theta c\psi -c\varphi s\psi & s\varphi s\theta s\psi +c\varphi c\psi & s\varphi c\theta \\ c\varphi s\theta c\psi +s\varphi s\psi & c\varphi s\theta s\psi -s\varphi c\psi & c\varphi c\theta \end{array}\right],\hfill \end{array}\end{array}$$

(8)

where s and c stands for $sin$ and $cos$, respectively.

In case of control tasks, the reverse transformation is needed. The required matrix ${\mathbf{R}}_b^e$ transforms any given vector ${\mathbf{v}}^b$ from $\mathcal{BF}$ to $\mathcal{EF}$, yielding ${\mathbf{v}}^e\in \mathcal{EF}$. Since the ${\mathbf{R}}_e^b$ is an orthonormal matrix, it can be written that:

$$\begin{array}{c}\hfill {\mathbf{R}}_e^b={\left({\mathbf{R}}_b^e\right)}^T.\end{array}$$

(9)

In the notation, Euler (or Tait–Bryan) angles vector $\theta$ is as follows:

$$\begin{array}{c}\hfill \theta ={\left[\begin{array}{ccc}\varphi & \theta & \psi \end{array}\right]}^T.\end{array}$$

(10)

Euler angle vector $\theta$ from Equation (16) is used in the paper (in angular velocities derivation) for the definition of the transformation matrices.

3.2. Angular Velocities

Euler’s notation presented in the previous subsections is the base for the mathematical formulas for angular velocities of a coaxial quadrotor. The vector $\dot{\theta }$ (Euler rates) expresses the projections of the UAV’s angular velocities (written in the $\mathcal{BF}$) onto the axes of the $\mathcal{EF}$. It is an alternative notation to $\omega \in \mathcal{BF}$ (from gyroscopes one may directly measure the angular velocity $\omega$), but to obtain the Euler rates, one must project $\omega$ on axes of the $\mathcal{EF}$.

The algorithm for the Euler rates calculation has been described in [35] (for more information see also Remark no. 2 in Appendix A). The first step is the calculation of the angular velocity $\omega$ in a local coordinate system $\mathcal{BF}$:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \omega ={\left[\begin{array}{ccc}{\omega }_x& {\omega }_y& {\omega }_z\end{array}\right]}^T=\hfill \\ \hfill & =\left[\begin{array}{c}\dot{\varphi }\\ 0\\ 0\end{array}\right]+{\mathbf{R}}_x\left(\varphi \right)\left[\begin{array}{c}0\\ \dot{\theta }\\ 0\end{array}\right]+{\mathbf{R}}_x\left(\varphi \right){\mathbf{R}}_y\left(\theta \right)\left[\begin{array}{c}0\\ 0\\ \dot{\psi }\end{array}\right].\hfill \end{array}\end{array}$$

(11)

A convenient representation of (11) may be given in a matrix form:

$$\begin{array}{c}\hfill \omega ={\mathbf{P}}_e^b\dot{\theta },\end{array}$$

(12)

where:

$$\begin{array}{c}\hfill {\mathbf{P}}_e^b=\left[\begin{array}{ccc}1& 0& -s\theta \\ 0& c\varphi & s\varphi c\theta \\ 0& -s\varphi & c\varphi c\theta \end{array}\right].\end{array}$$

(13)

Since an inverse transformation ($\mathcal{BF}\to \mathcal{EF}$) is needed, the ${\mathbf{P}}_b^e$ matrix must be derived:

$$\begin{array}{c}\hfill {\mathbf{P}}_b^e={\left({\mathbf{P}}_e^b\right)}^{-1}=\left[\begin{array}{ccc}1& s\varphi \frac{s\theta }{c\theta }& c\varphi \frac{s\theta }{c\theta }\\ 0& c\varphi & -s\varphi \\ 0& \frac{s\varphi }{c\theta }& \frac{c\varphi }{c\theta }\end{array}\right].\end{array}$$

(14)

The matrix ${\mathbf{P}}_b^e$ is used directly to determine the Euler rates:

$$\begin{array}{c}\hfill \dot{\theta }={\mathbf{P}}_b^e\omega .\end{array}$$

(15)

Note that in case of multirotor sensory system, the value of $\omega$ is provided from measurements.

3.3. Angular Acceleration

In the derivation it is assumed that the UAV is considered as a perfectly stiff object in the 3D space. That allows the description of the angular acceleration of the UAV to use the basic Euler equation for the motion of a rigid body ([34], p. 168), and in general form one can write that:

$$\begin{array}{c}\hfill \frac{d\mathbf{L}}{dt}+\omega \times \mathbf{L}=\tau ,\end{array}$$

(16)

where angular momentum $\mathbf{L}$ for a rigid body is

$$\begin{array}{c}\hfill \mathbf{L}=\mathbf{I}\omega .\end{array}$$

(17)

Based on the constancy of the inertia moments and by substituting (17) to (16), one obtain the torque equation:

$$\begin{array}{c}\hfill \tau =\mathbf{I}\dot{\omega }+\omega \times \left(\mathbf{I}\omega \right).\end{array}$$

(18)

Using (18), the angular acceleration is written as:

$$\begin{array}{c}\hfill \dot{\omega }={\mathbf{I}}^{-1}(\tau -\omega \times \left(\mathbf{I}\omega \right)).\end{array}$$

(19)

Equation (29), however, due to the gyroscopic effect, viscous friction and several drag-like effects, must be extended (described in the following subsections) to the complete form presented in subsection Equation (1).

3.3.1. Input Torques

The forces generated by the propellers can be written as

$$\begin{array}{c}\hfill {\mathbf{f}}_i={\left[\begin{array}{ccc}0& 0& b_i{\omega }_i^2\end{array}\right]}^T,\end{array}$$

(20)

where ${\mathbf{f}}_i\in \mathcal{BF}$ and $b_i$ are the force and thrust factor of the i-th propulsion unit, respectively.

The input torque vector $\tau$ is a sum of torques (${\tau }_i$) due to subsequent propellers. A particular i-th torque defined as a cross-product of the i-th propeller position ($p_i$) in respect to the $\mathcal{BF}$ (22) and the i-th force (20) can be calculated as

$$\begin{array}{c}\hfill {\tau }_i={\left[\begin{array}{ccc}{\tau }_{x_i}& {\tau }_{y_i}& {\tau }_{z_i}\end{array}\right]}^T={\mathbf{p}}_i\times {\mathbf{f}}_i,\end{array}$$

(21)

$${\mathbf{p}}_{1/2}={\left[\begin{array}{ccc}{l}_w& l_w& \mp h\end{array}\right]}^T,$$

(22)

$${\mathbf{p}}_{3/4}={\left[\begin{array}{ccc}-l_w& l_w& \mp h\end{array}\right]}^T$$

$${\mathbf{p}}_{5/6}={\left[\begin{array}{ccc}-l_w& -l_w& \mp h\end{array}\right]}^T,$$

$${\mathbf{p}}_{7/8}={\left[\begin{array}{ccc}{l}_w& -l_w& \mp h\end{array}\right]}^T,$$

where $l_w$=$\frac{\sqrt{2}}{2}l$ and l-arm’s length in the UAV cross-frame, h-distance between the arm and propeller’s plane in respect to the Z axis of $\mathcal{BF}$.

Equation (21) covers the torques generated in the X and Y axes in $\mathcal{BF}$, while the Z axis will also suffer from the reaction torques, given in (23):

$${\tau }_z=d_1({\omega }_2^2+{\omega }_3^2+{\omega }_6^2+{\omega }_7^2)-d_2({\omega }_1^2+{\omega }_4^2+{\omega }_5^2+{\omega }_8^2).$$

(23)

In Equation (23), the parameters $d_i$ refer to the reaction torque gains. Two of them are introduced, which stems from the opposite type of propellers system. Propellers 1, 4, 5, 8 rotate clockwise (CW) and propellers 2, 3, 6, 7 rotate counterclockwise (CCW) (Figure 6). Note that CW rotation contributes to counterwise reaction torque.

$$\begin{array}{c}\hfill \dot{\omega }={\mathbf{I}}^{-1}(\tau -\omega \times \left(\mathbf{I}\omega \right)).\end{array}$$

(24)

3.3.2. Gyroscopic Effect

The gyroscopic effect is described by the basic equation of the gyroscope:

$$\begin{array}{c}\hfill {\tau }_{gyro}={\omega }_p\times {\mathbf{L}}_r,\end{array}$$

(25)

where: ${\mathbf{L}}_r$—rotors–propellers angular momentum vector, ${\omega }_p$—angular velocity precession vector.

Each object with a spinning mass is subjected to a gyroscopic precession motion effect, so it can be written that

$$\begin{array}{c}\hfill {\omega }_p=\omega .\end{array}$$

(26)

Similarly to the angular momentum of the robot frame, given in Equation (17) the angular momentum of the propulsion unit can be expressed as follows:

$$\begin{array}{c}\hfill {\mathbf{L}}_r=I_r{\omega }_r.\end{array}$$

(27)

In Equation (27), $I_r$ scalar denotes simplified motor–propeller inertia, while ${\omega }_r$ refers to the sum of the rotors’ angular velocities (in the Z-axis). The angular momentum vector (27) is obtained by multiplying the moment of inertia ($I_r$) by the rotor’s angular velocity about the axis where the rotation occurs (${\omega }_r$):

$${\omega }_r=\left[\begin{array}{c}0\\ 0\\ {\omega }_1-{\omega }_2-{\omega }_3+{\omega }_4+{\omega }_5-{\omega }_6-{\omega }_7+{\omega }_8\end{array}\right].$$

(28)

The rotor angular velocity vector (${\omega }_r$) contains zeros for the components in the axes X, Y and the sum of the rotors angular velocities ${\omega }_i$ (where i varies in the range of $[1,8]$) in the Z axis. Knowing ${\mathbf{L}}_r$ and ${\omega }_r$, Equation (25) can be rewritten as follows:

$$\begin{array}{c}\hfill {\tau }_{gyro}=\left[\begin{array}{c}{\omega }_x\\ {\omega }_y\\ {\omega }_z\end{array}\right]\times \left[\begin{array}{c}0\\ 0\\ I_r{{\omega }_r}_z\end{array}\right]=\left[\begin{array}{c}{\omega }_y{I}_r{{\omega }_r}_z\\ -{\omega }_x{I}_r{{\omega }_r}_z\\ 0\end{array}\right].\end{array}$$

(29)

3.3.3. Viscous Friction

According to [36], friction model for air flows has to include viscous friction. In a standard form, with respect to the UAV’s angular velocity, viscous friction can be expressed as (30):

$$\begin{array}{c}\hfill {\tau }_t=b_w\left(\omega \right)\left(\omega \right)\omega =\left[\begin{array}{c}{b}_w\left({\omega }_x\right){\omega }_x^2\\ b_w\left({\omega }_y\right){\omega }_y^2\\ b_w\left({\omega }_z\right){\omega }_z^2\end{array}\right].\end{array}$$

(30)

3.3.4. Drag-Like Effects

As it may be found in [37], there are a number of drag-like effects: blade flapping, rotor drag, translational drag, profile drag and parasitic drag (which may be ignored for speeds up to 10 $\mathrm{m}/\mathrm{s}$ [37]).

The blade flapping effect is especially noticeable on multirotor UAVs, when its rotors undergoes translational motion. This effect is seen on the rotor’s tip path plane which deviates by an angle called as flapping angle ${\beta }_{flap}$. More information about the flapping effect can be found, among others in [38], and further in the article, the final form of the equation for the i-th rotor is introduced:

$$\begin{array}{c}\hfill {\delta }_i={\mathbf{f}}_i\left({\mathbf{A}}_{flap}\frac{{\mathbf{v}}_p}{{\omega }_i}+{\mathbf{B}}_{flap}\frac{\omega }{{\omega }_i}\right),\end{array}$$

(31)

where ${\mathbf{v}}_p={\left[\begin{array}{ccc}{v}_x& v_y& 0\end{array}\right]}^T$ is horizontal velocity in $\mathcal{BF}$, and ${\mathbf{A}}_{flap}$, ${\mathbf{B}}_{flap}$ are lumped parameter matrices defined as:

$$\begin{array}{c}\hfill {\mathbf{A}}_{flap}=\frac{1}{r}\left[\begin{array}{ccc}-A_{1c}& A_{1s}& 0\\ -A_{1s}& -A_{1c}& 0\\ 0& 0& 0\end{array}\right],\end{array}$$

(32)

$$\begin{array}{c}\hfill {\mathbf{B}}_{flap}=\left[\begin{array}{ccc}\hfill -B_2& \hfill B_1& \hfill 0\\ \hfill B_1& \hfill -B_2& \hfill 0\\ \hfill 0& \hfill 0& \hfill 0\end{array}\right],\end{array}$$

(33)

where r—rotor radius, $A_{1c}$ and $A_{1s}$ are positive scalar constants that depend on the geometry of the rotor blades (see Equation (4.45) and Equation (4.47) in [39]). In practice agile coaxial quadrotors require stiff rotors in order that rotation in the airframe is directly translated into rotation in the thrust vector. That is to reduce the magnitude of ${\mathbf{B}}_{flap}$ [40]. Note that parameters: $A_{1c}$, $A_{1s}$, $B_1$, $B_2$ needs to be identified from flight tests [38]. In the paper [29], the authors provided example of the identification of matrix parameters using the Particle Swarm Optimization method for a typical quadrotor attached on a special testbed that allows to record measurement data for the used propulsion units with propellers, on which drag-like effects were analyzed.

The second important drag-like effect is so-called rotor drag, which usually occurs in small multirotor UAV, where its rigid or semirigid rotor blades causes an inbalance in flapping forces. Mathematically it is modeled as:

$$\begin{array}{c}\hfill {\mathbf{D}}_I=K_I{\mathbf{v}}_p,\end{array}$$

(34)

where $K_I$ is a scalar parameter and $v_p$-planar velocity.

Another effect depends on speed in the air and is called translational drag. It is observed when the air is redirected from any translational movement to the downward of $\mathcal{BF}$. During slow movements (according with [37]), this effect can be modelled via following equation:

$$\begin{array}{c}\hfill {\mathbf{D}}_T=K_{T_1}{\mathbf{v}}_p.\end{array}$$

(35)

Considering fast maneuvers, one needs to add the air velocity $v_i$ induced by the i-th rotor:

$$\begin{array}{c}\hfill {\mathbf{D}}_T=K_{T_2}{(-v_z+v_i)}^4{\mathbf{v}}_p,\end{array}$$

(36)

where $K_{T_1}$, $K_{T_2}$ are scalar parameters, and $v_z$ is translational the speed in z axis.

There are two more important drag-like effects: parasitic drag and profile drag (37). The first one comes from nonlifting surfaces of the multirotor UAV only for high speeds (>10 $\mathrm{m}/\mathrm{s}$). Therefore one can ignore it in the case considered in this paper [37] and applied in to the model just the profile drag. This drag effect is gained with the transverse velocity of the rotor blades in the air and can be formalized as:

$$\begin{array}{c}\hfill {\mathbf{D}}_P=K_P{\mathbf{v}}_p,\end{array}$$

(37)

where $K_P$ is a lumped parameter.

After considering all of the above-mentioned drag effects, the total drag effect force acting on the i-th rotor is:

$$\begin{array}{c}\hfill {\mathbf{D}}_i={\delta }_i+{\mathbf{D}}_I+{\mathbf{D}}_T+{\mathbf{D}}_P,\end{array}$$

(38)

and the torques corresponding to them, for the i-th rotor are:

$$\begin{array}{c}\hfill {\tau }_{D_i}={\mathbf{p}}_i\times {\mathbf{D}}_i.\end{array}$$

(39)

The net torque ${\tau }_D$ is the sum of all ${\tau }_{{\mathbf{D}}_i}$, where $i=1,\dots ,8$ and $i\in \mathbb{Z}$.

3.3.5. Orientation Model—Angular Acceleration

The final form of upgraded Euler’s equation given as (1) is:

$$\begin{array}{cc}\hfill & \tau =I\dot{\omega }+\omega \times \left(\mathbf{I}\omega \right)+{\tau }_{gyro}+{\tau }_d+{\tau }_t=I\dot{\omega }+\omega \times \left(\mathbf{I}\omega \right)+\hfill \\ \hfill & +\omega \times I_r{\omega }_r+\sum _{i=1}^8{\tau }_{d_i}+b_w\left(\omega \right)\left(\omega \right)\omega .\hfill \end{array}$$

(40)

In (40), the resulting quantity is the angular acceleration, so after simple algebra it may be written that:

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \dot{\omega }=I^{-1}(\tau -\omega \times \left(\mathbf{I}\omega \right)-\omega \times I_r{\omega }_r+\hfill \\ \hfill & -\sum _{i=1}^8{\tau }_{d_i}-b_w\left(\omega \right)\left(\omega \right)\omega ).\hfill \end{array}\end{array}$$

(41)

Analysis of previous and present sections sugests that using measurements of angular accelerations from gyroscopes (to retrieve the UAV’s orientation in $\mathcal{EF}$), one needs to integrate angular acceleration and after that, transform obtained result into the Euler rates (15). Note that Euler (or Tait–Bryan) angles are given by the second integration of the Euler rates.

3.4. Linear Acceleration

The basic equation for linear acceleration includes the sum of forces acting on a body with a certain mass, which stems from Newton’s law. Since the UAV can rotate and translate in 3D, the dynamics of both of these types of motion has to be captured. There are several types of forces explained below. All of them act directly on the UAV in $\mathcal{BF}$.

3.4.1. Thrust Force

Thrust force or simply input force is the sum of all forces generated by the propulsion units. If the force of the i-th rotor is given as ${\mathbf{f}}_i$:

$$\begin{array}{c}\hfill {\mathbf{f}}_i={\left[\begin{array}{ccc}0& 0& b_i{\omega }_i^2\end{array}\right]}^T,\end{array}$$

(42)

then the total thrust force vector $\mathbf{f}$ is as follows:

$$\begin{array}{c}\hfill \mathbf{f}=\sum _{i=1}^8{\mathbf{f}}_i.\end{array}$$

(43)

Note that $b_i$ is the thrust coefficient of the i-th propeller. These parameters may change according to the spatial location of the propellers (especially in coaxial multirotors), and what is more, they can differentiate the type of flight, for example the ground effect can be reflected in a proper change of $b_i$.

3.4.2. Gravitational Force

The force of gravity is normally provided in $\mathcal{EF}$. Since the NED is chosen as the configuration of $\mathcal{EF}$, the gravity force is as follows:

$$\begin{array}{c}\hfill {\mathbf{f}}_g^e={\left[\begin{array}{ccc}0& 0& mg\end{array}\right]}^T.\end{array}$$

(44)

In $\mathcal{BF}$, a projection of ${\mathbf{f}}_g^e$ is used:

$$\begin{array}{c}\hfill {\mathbf{f}}_g={\mathbf{f}}_g^b={\mathbf{R}}_e^b{\mathbf{f}}_g^e=\left[\begin{array}{c}-mgsin\theta \\ mgsin\varphi cos\theta \\ mgcos\varphi cos\theta \end{array}\right].\end{array}$$

(45)

3.4.3. Coriolis Force

The Coriolis force is a pseudo force, which occurs only in the local coordinate system:

$$\begin{array}{c}\hfill {\mathbf{f}}_{co}=-2m\omega \times \mathbf{v}.\end{array}$$

(46)

3.4.4. Drag Forces

Drag forces are due to aerodynamic effects, explained in the previous subsection. The ${\mathbf{f}}_d$ is a sum of all drag-like effects of all propulsion units:

$$\begin{array}{c}\hfill {\mathbf{f}}_d=\sum _{i=1}^8{\mathbf{D}}_i.\end{array}$$

(47)

3.4.5. Friction in Linear Motion

The friction force ${\mathbf{f}}_t$ is results from the viscous drag:

$$\begin{array}{c}\hfill {\mathbf{f}}_t=b_v\left(\mathbf{v}\right)\left(\mathbf{v}\right)\mathbf{v}=\left[\begin{array}{c}{b}_v\left(v_x\right)v_x^2\\ b_v\left(v_y\right)v_y^2\\ b_v\left(v_z\right)v_z^2\end{array}\right].\end{array}$$

(48)

3.4.6. Linear Motion Model—Linear Acceleration

The overall model of linear motion is given by (49). Note that the left side of the equation is a linear acceleration in the local coordinate system, i.e., $\dot{\mathbf{v}}\in \mathcal{BF}$.

$$\begin{array}{c}\hfill \begin{array}{cc}\hfill & \dot{\mathbf{v}}=\frac{1}{m}(\sum _{i=1}^8{\mathbf{f}}_i-{\mathbf{R}}_e^b{\mathbf{f}}_g^e+2m\omega \times \mathbf{v}+\hfill \\ \hfill & -\sum _{i=1}^8{\mathbf{D}}_i-b_v\left(\mathbf{v}\right)\left(\mathbf{v}\right)\mathbf{v}).\hfill \end{array}\end{array}$$

(49)

It is important to mention that the presented mathematical model is still just a simplification, which means that in general, the main sources of potential uncertainties and disturbances need to be considered, analyzed and remembered in its application for control purposes. The main source of disturbances is the flapping effect, while the key role in the model uncertainty is played by thrust coefficients which change in time, and dynamic states due to air stream coupling between coaxial propellers. The impact and scale of uncertainty, as well as unmeasured (stochastic) disturbances have been assessed in experiments on data from a real robot (Section 5, Experiment no. 3).

4. Model Validation

4.1. Introduction

In the following part, the paper presents test results obtained for the Falcon V5 flying robot (see Figure 7), which is the octoquad based on a X-shape skeleton construction made of Carbon Fiber Reinforced Polymer (CFRP). It has eight propulsion units protected by additional covers. Based on conducted research results [41], it was decided that the robot should be driven by BLDC engines from RC Tiger Motors (model: MN10–15) in conjunction with the company’s electronic speed controllers (model: Opto–Pro) and 10 × 3.3 carbon fiber propellers. The main computing unit was based on the ARM microcontroller (Cortex M4F family). The sensing system uses the ADIS16488 sensor and the additional BMP085 pressure sensor (for more accurate altitude and vertical speed estimation of the robot).

The proposed model of the Falcon V5 was implemented in Matlab 2015a in order to verify the adequacy of the model behaviour to simulate dynamic properties of a real coaxial quadrotor (i.e., correlation between the current orientation and position of model and real UAVs). For this purpose, a series of experiments was conducted.

In the first stage, it was necessary to obtain the values of parameters in the mathematical equations. By having access to the full technical documentation of the Falcon V5 robot, the authors could use its 3D model developed in the Inventor software (Figure 6a) to determine all the necessary geometric dimensions, masses and moments of inertia. The transparency of the model is guaranteed by the knowledge of:

• the total mass of the flying robot (m),
• the mass of the individual propulsion unit ($m_w$),
• the arm’s length in the robot cross frame (l),
• the distance between propeller’s plane and the arm in respect to the Z axis of the $\mathcal{BF}$ (h),
• moments of inertia in X, Y, Z axes, respectively: $I_{xx}$, $I_{yy}$, $I_{zz}$, $I_{xy}$, $I_{yz}$, $I_{xz}$, $I_w$,
• coefficients ($b_1$, $b_2$, $d_1$, $d_2$, ${\mathbf{A}}_{flap}$, ${\mathbf{B}}_{flap}$) used in input torques, viscous friction and drag effects models.

For the purpose of further tests, appropriate measurements were gathered and numerical values of coefficients were identified (see [31]). For control of position and orientation in the 3D space, one can use six PID controllers. This solution is still the most popular in coaxial quadrotors because of its universality, tuning simplicity and effectiveness. In the research the controller parameters are tuned via a heuristic algorithm, i.e., the Particle Swarm Optimization (PSO) method. All necessary information about the PSO method is found in [29]. The coaxial quadrotor’s PID controller tuning is outside the scope of this paper, but to give a brief insight into the main configuration parameters of the PSO algorithm and the PID tuning strategy, some of these are shown.

During the optimization, 18 different gain values were searched for i.e., the gains $K_P$, $K_I$ and $K_D$ of six controllers. Due to the complexity of the search space, the optimization was performed as follows:

• Simultaneous search for controller gains for orientation in the axes: X, Y and Z, where the minimized cost function is defined as a sum of squared errors between the set and obtained angle (note that $\varphi \in [-\pi /2,\pi /2]$),
• Search for the Z-axis position controller gains. The minimized cost function is a sum of squared errors between the set and measured altitude.
• Search for the position controller gains in the X and Y axes, the minimized cost function is a sum of squared errors between the set and measured position.

The PSO basic velocity equation was prepared according to the original one [42], with modified weights:

$$\begin{array}{cc}\hfill & v_k^i=\xi \left(\omega v_{k-1}^i+c_1\ast rand(pbes{t}^i-presen{t}_{k-1}^i)\right)\hfill \\ \hfill & +c_2\ast rand(gbes{t}^i-presen{t}_{k-1}^i),\hfill \end{array}$$

(50)

where $c_1=c_2=1.8$, $\xi =0.73$ and $\omega =1.2$.

4.2. Experiment No. 1

The behaviour of the coaxial quadrotor model within the time horizon equal to 50 $\mathrm{s}$ was analyzed. At the moment of beginning, the local reference system coincided with the system of the observer. The aim of the test was to maintain the desired altitude of 3 $\mathrm{m}$ at simultaneously linear change of position in the other reference system axis. The behaviour of the inner control loop (angles) and changes in rotational speed of particular propulsion units of the robot were analyzed. It was also assumed that the coaxial quadrotor cannot rotate about the z axis. The rotational speeds of particular propulsion units were limited to 1200 $\mathrm{rad}/\mathrm{s}$ (equivalent of $11,459.16$$\mathrm{RPM}$).

Figure 8 shows the results obtained in the first 12 $\mathrm{s}$ of flight. The applied control system allowed the quadrotor (after about 3 $\mathrm{s}$) to rise to the fixed altitude with simultaneous adjustment of $\varphi$ and $\theta$ angles.

4.3. Experiment No. 2

Experiment no. 2 extended the previous test, ie. in the 21 $\mathrm{s}$ of flight (at the altitude of 3 $\mathrm{m}$) a step change of the set position in the Y axis was forced in the opposite direction to the current flight. Furthermore, the rotation of the coaxial quadrotor about the $\psi$ angle was introduced during the movement of the robot relative to the z axis. The experiment results are provided in Figure 9 and Figure 10.

The stated goals of these two tests allowed us to evaluate both the behaviour of a real quadrotor through the implementation of its mathematical model, as well as the effectiveness of control of continuous rotation of the coaxial quadrotor during the change of its position in 3D (without a cascade control system, manual maneuverability by the change of quadrotor’s angles and thrust is a challenge in itself). Experiment results highlight that the rapid changes of rotational speed of particular UAV’s units allow us to effectively compensate the tilt of robot after the change of set position. Naturally, the effectiveness of the control system can be increased (for example, through the selection of a different strategy for controllers tuning or using a different type of controller than the one used for tests).

4.4. Experiment No. 3 (with Real-World Data from the Falcon V5 Flying Robot)

Tests of the proposed coaxial quadrotor model would not be complete without the validation with real, measured data recorded from the Falcon V5 flying robot. The experiment of length of 40 $\mathrm{s}$ was performed (see Supplementary Materials). The semi-autonomously controlled Falcon V5 moved to a fixed altitude (2 $\mathrm{m}$). Then he flew forward about 4–5 m and return to the base point on the set altitude. In the next stage, the UAV was tasked to fly about 4–5 m to the left and return to the base point, ending the attempt with landing near the starting point. In straight flight, the altitude was changed simultaneously, and in flight from side to side—it remained unchanged.

In Figure 11, angular accelerations (in x,y,z ∈$\mathcal{EF}$ coordinates) calculated for the Falcon V5 UAV model and real, measured data (i.e., the derivative of measured angular velocity), are shown. Figure 12 presents, respectively, forces (in $\mathit{x},\mathit{y},\mathit{z}$∈$\mathcal{EF}$ coordinates) calculated from the model and measured by the accelerometer (precisely: accelerometer measurements after filtration and subtraction of the force of gravity projected on its axies). It should be remembered that during the hovering (when the force of gravity is equal to the thrust), the accelerometer does not show the state of equilibrium (zero on axes). It returns the value of the force of gravity.

The obtained results are promising, but it should be noted that the UAV model requires continuous modification of rotors thrusts gains [31]. Otherwise one will not be able to obtain precise estimates of model parameters. The experiments conducted by the authors did not cover this estimation issue.

5. Summary

It is said that modeling is the key to a better understanding of the world around us. The more difficult and more complex behavioral aspect is analyzed, the more complex mathematical-physical formulas accompan its description.

From the perspective of efficient control algorithms design, it is necessary to build models that are flexible enough to be compared to one another and to be simplified according to predefined requirements. The model of coaxial quadrotor proposed in this paper fulfills this task as it can be used effectively both to analyze the orientation of the robot in the 3D space and to verify the quality of the estimation algorithms [31], as well as for control design.

The obtained results allowed the authors to plan further research—specifically, on increase of the universality of coaxial quadrotor model. These studies include modeling of the dynamics of drive units mounted in pairs (so far approximated in tests by simple, linear, inertial mathematical models) and research on the effects of aerodynamic behaviour of UAVs during the flight by use of such a solution. Initial trials and results are presented in [41].

The second equally interesting, priority direction of research already begun by the authors is a quantitative investigation the role of small aerodynamic effects (usually neglected are known to be small) on controllers (especially relying only on basic inertial measurements).