MATHEMATICAL MODELLING OF UNMANNED AERIAL VEHICLES WITH FOUR ROTORS

Mathematical model of an unmanned aerial vehicle with four propulsors (quadcopter) is indispensable in quadcopter movement simulation and later modelling of the control algorithm. Mathematical model is, at the same time, the first step in comprehending the mathematical principles and physical laws which are applied to the quadcopter system. The objective is to define the mathematical model which will describe the quadcopter behavior with satisfactory accuracy and which can be, with certain modifications, applicable for the similar configurations of multirotor aerial vehicles. At the beginning of mathematical model derivation, coordinate systems are defined and explained. By using those coordinate systems, relations between parameters defined in the earth coordinate system and in the body coordinate system are defined. Further, the quadcopter kinematic is described which enables setting those relations. Also, quadcopter dynamics is used to introduce forces and torques to the model through usage of Newton-Euler method. Final derived equation is Newton’s second law in the matrix notation. For the sake of model simplification, hybrid coordinate system is defined, and quadcopter dynamic equations derived with the respect to it. Those equations are implemented in the simulation. Results of behavior of quadcopter mathematical model are graphically shown for four cases. For each of the cases the propellers revolutions per minute (RPM) are set in a way that results in the occurrence of the controllable variables which causes one of four basic quadcopter movements in space.


INTRODUCTION
Quadcopter (or quadrotor) is an unmanned aerial vehicle with four propulsors that are enabling vertical takeoff and landing.It has six degrees of freedom -6 DOF and four controllable variables which makes it an underactuated and dynamically unstable system.Number of controllable variables is equal to the number of propulsors which affect position and attitude of quadcopter in space.Quadcopter cannot move translationally without the finite rotation around one of the axes, i.e. without the inclination of quadcopter.To achieve the inclination it is necessary to change the torque with respect to the one of the axes.In order to do so, it is necessary to increase or decrease the thrust on one or two propulsors.If the change of thrust happens on just one propulsors, that could cause instability in torque around the Z B axis of rotation.To achieve stable flight, it is necessary to combine several high accuracy sensors with fast and robust control algorithm.
Only moving parts on quadcopter are propellers on propulsors and which are fixed in a propulsors axis.Quadcopter can have cross (+) configuration (X B and Y B axes are oriented in the directions of propulsors), or it can have X configuration (X B and Y B axes are oriented in the directions between the propulsors).For further mathematical modelling, + configuration is assumed.Quadcopter frame is a symmetric, light and thin construction that mechanically connects propulsors.Propulsor motor and propeller are directly connected, with all the propulsors axes being fixed and parallel.Propeller rotation causes airflow in the negative direction of the Z B axis which results in thrust in the positive direction of the Z B axis.Quadcopter frame is assumed to be rigid.The only thing that has direct influence on the quadcopter movement are each motor's RPM.

MATHEMATICAL MODELLING OF QUADCOPTER
Mathematical model describes quadcopter movement and behavior with the respect to the input values of the model and external influences on quadcopter.Mathematical model can be observed as a function that is mapping inputs on outputs.By using mathematical model, it is possible to predict position and attitude of quadcopter by knowing the four angular velocities of propellers, i.e. it enables computer simulation of quadcopter behavior in different conditions.Computer simulation is relatively simple, cheap and harmless method for control algorithm verification.More detailed mathematical model describes quadcopter behavior more accurately, but it also requires more computer resources which leads to longer simulation time or even to inability to successfully complete the simulation.Depending on the requirements, it is necessary to find compromise between model complexity and accuracy.For the purpose of mathematical model, quadcopter movement is directly dependent on angle velocities of the propellers.Quadcopter frame with propulsors are the only elements that will come into the consideration in this mathematical model.Quadcopter structure along with the body frame and corresponding angular velocities are shown in Figure 1.

QUADCOPTER COORDINATE SYSTEMS
To get the mathematical model, it is necessary to define two coordinate systems:  Earth fixed frame (E-frame, ℱ E )  Body fixed frame (B-frame, ℱ B ) Some quadcopter physical properties are measured in ℱ E (roll, pitch and yaw angles, angular velocities), while some properties are measured in ℱ B (linear accelerations).
ℱ E is the inertial right-handed coordinate system where positive direction of  E axis is in the direction from the earth.Quadcopter position  and attitude  are defined in ℱ E .ℱ B is fixed on quadcopter body.Positive direction of the  B axis goes through the propulsor 1 which is located on the front side of quadcopter.Positive direction of the  B axis goes through the propulsor 4 which is located on the left side of quadcopter.The  B axis is perpendicular to  B and  B axes and its positive direction is in the direction of propulsors thrust forces.Assumption is that the origin of ℱ B coincides with the center of gravity of the quadcopter.Linear velocities  B , angular velocities  B , forces  B and torques  B are defined in ℱ B .Quadcopter position is defined with vector  between origins of ℱ E and ℱ B (Fig. 2) Quadcopter attitude  is defined with the orientation of ℱ B with the respect to the ℱ E .The orientation is defined with three consecutive rotations around the ℱ E coordinate axes.Roll-pitch-yaw order is applied Motion equations are more suitable to formulate with the respect to the ℱ B for several reasons: system inertia matrix is time-invariant, equations simplification because of quadcopter frame symmetry, sensors measurements are easily converted to ℱ B and control variables equations simplification.

QUADCOPTER KINEMATICS
Kinematics of a rigid body with 6 DOF is given with ̇=  , Because of the need to transform measured values from one coordinate system to another, the rotation matrix is introduced, which by matrix multiplication, transfers linear velocity vector from one coordinate system to another.Matrix  is the orthogonal matrix.
Angles and angular velocities are measured in ℱ E .Matrix  is the transformation matrix that transfers angular velocities from To transfer angular velocities from ℱ E to ℱ B , angular velocity vector in ℱ E has to be matrix multiplied by inverse of the transformation matrix.

QUADCOPTER DYNAMICS: NEWTON-EULER METHOD
Quadcopter dynamics is described by differential equations that were derived by using the Newton-Euler method.Dynamics of a rigid 6 DOF body takes into consideration the mass m and the inertia of the body I.By applying the assumption that the quadcopter frame has symmetrical structure where four propulsors carriers are in line with X B and Y B axes, i.e. the principal inertia axes coincides with the ℱ B coordinate axes, inertia matrix becomes the diagonal matrix where Quadcopter dynamics is described by where  3×3 is the identity matrix of size 3  3 (different from I), ̇B is the linear acceleration vector, ̇B is the angular acceleration vector,  B is the force vector acting on quadcopter,  B is the torque vector acting on quadcopter, all with the respect to the ℱ B .

Generalized force vector 𝛌
Now, equation ( 10) can be written in the form of  B ̇+  B ()  = , (12) where ̇ is the generalized acceleration vector,  B is the system inertia matrix,  B () is the Coriolis-centripetal matrix, all with the respect to the ℱ B .System inertia matrix is diagonal matrix comprised from quadrotor mass and moments of inertia with the respect to the coordinate axes of ℱ B .Because the quadcopter dynamics is described in ℱ B , it is important to introduce the Coriolis-centripetal matrix to the model which will upgrade the model with centripetal force that is acting on the rotational body and also with the Coriolis effect which affects moving body in rotational frame.Equation ( 12) generally applies to all rigid bodies to which the made assumptions can be applied (symmetry and same origin for ℱ B and COG).
Force vector  can be divided into three components: gravitational vector, gyroscopic torque vector and movement vector.

Gravitational vector
The first component of  is gravitational vector which introduces the gravitational force to the model.It only affects the linear components of the model, and not the angular.In ℱ E , the gravitational force always has the direction in the negative direction of  E , if the ℱ  coordinate axes are set in a way previously described.Considering that dynamics is described in ℱ B , it is necessary to transfer gravitational vector from ℱ E to ℱ B by matrix multiplication with .Due to quadcopter movement, Z B is not always parallel with  E axis.In that case, only a certain component of gravity vector will act in the direction of Z B , depending of the roll and pitch angles ( and ) where  G B is the gravitational vector with respect to ℱ B , and  G E is the gravitational vector with respect to ℱ E .

Gyroscopic torque vector
Gyroscopic torque vector introduces the gyroscopic effect to the model.It manifests itself in a form of a rotation around unwanted axis, which is perpendicular to the axis of the propulsor and the axis around which the wanted rotation is being achieved

Movement vector
Movement vector introduces thrust force of propulsors  1 and torques around coordinate axes ℱ B :  2 ,  3 and  4 .Number of controllable variables is equal to the number of propulsors Four basic movements enable the quadcopter to reach defined altitude and attitude.This is the only vector that can be directly influenced, which together with previously named vectors influence to quadcopter accelerations.
Movement vector is influenced through propellers angular velocity   , which is, naturally, influenced by propulsors RPM.
where   is angular speed of the propeller , and  is the thrust coefficient (measurement unit Ns 2 ) which can be calculated according the equation where  T is thrust factor,  is air density,  is the area of propeller disk and  is the propeller radius.
Controllable variable  2 is the roll torque  ϕ around X B axis.By increasing the 4 th propeller angular velocity and decreasing 2 nd propeller angular velocity, positive quadcopter rotation with respect to the X B axis is achieved, which causes the linear quadcopter movement in the negative direction of the Y B axis.Differences in the angular velocities are determined in a way which does not change vertical thrust, i.e. during the quadcopter rotation around X B axis, there cannot be any unwanted changes of the quadcopter altitude.Controllable variable  2 equation is where  is the distance between the quadcopters Center Of Gravity (COG) and propulsors axis.
Controllable variable  3 is the pitch torque  θ around Y B axis.By increasing the 3 rd propeller angular velocity and decreasing 1 st propeller angular velocity, positive quadcopter rotation with the respect to the Y B axis is achieved, which causes the linear quadcopter movement in the positive direction of the X B axis.Similar to the  2 controllable variable, the differences in the angular velocities are determined in a way which does not change vertical thrust, i.e. during the quadcopter rotation around Y B axis, there cannot be any unwanted changes of the quadcopter altitude.Controllable variable  3 equation is Controllable variable  4 is the yaw torque  ψ around Z B axis.By increasing the angular velocity of a pair of propellers (on 2 nd and 4 th propulsor) and by decreasing the angular velocity of a pair of propellers (on 1 st and 3 rd propulsor), positive rotation with respect to the Z B axis is achieved thanks to the fact that one pair of propellers (on 2 nd and 4 th propulsor) is rotating clockwise, while other pair of propellers (on 1 st and 3 rd propulsor) is rotating counterclockwise.Created resulting moments causes unbalance in the sum of moments with respect to Z B , which is the cause of rotation.Differences in the angular velocities are determined in a way which does not change vertical thrust.Controllable variable  4 equation is where  is the drag coefficient which can be calculated according the equation  =  P    3 , (21) where  P is the power factor.
Considering aerodynamic effects, it follows that forces and moments are proportional to the squared angular velocities of the propellers.Movement vector  B () is represented by the product of movement matrix  B and the vector of the squared angular velocities of the propellers The system inertia matrix inverse is multiplied by the sum of all the forces and torques, relevant to the mathematical model, acting on the quadcopter.

NEWTON-EULER EQUATIONS IN HYBRID COORDINATE SYSTEM (𝓕 𝐇 )
Quadcopter dynamics equation ( 24) is written with respect to ℱ B .It can be useful to rewrite the dynamics equations with respect to ℱ H which is comprised from the linear accelerations equations with respect to ℱ E and angular accelerations equations with respect to ℱ B .In such coordinate system, it is relatively easy to show the dynamics of the linear accelerations in the combination with regulation.It is especially suitable for showing the vertical position in ℱ E .

Generalized velocity vector with respect to
Quadrotor dynamics with respect to ℱ H in matrix form    ̇+   ()  =   () +   ()  +     , (26 where  ̇ is the generalized acceleration vector with respect to ℱ H . System inertia matrix is defined in a similar way to the one defined with respect to ℱ B .Coriolis-centripetal matrix defined with respect to ℱ H excludes members that affects linear accelerations (because they are defined with respect to ℱ E , and not in the rotational ℱ B ). Gravitational vector is defined with respect to ℱ E , so the gravitational force is influencing only the component in the direction of the  E axis.Gyroscopic effects that occurs as the result of propellers rotation are the same as in the mathematical model with respect to ℱ B .Movement vector defined with respect to ℱ H is different than the one defined in ℱ B .Controllable variable  1 influences all three linear accelerations equation through the rotation matrix .
By modifying the equation (26), generalized acceleration vector with respect to ℱ H can be calculated There are many effects that were not included in the derivation of this mathematical model because of the rise of the model complexity and because they have little significance to the model accuracy.

SIMULATION RESULTS
The behavior of the outputs of the mathematical model is dependent on the input values.
Even though the mathematical model deals with angular velocities as one of the variables, the input values are chosen to be the RPM of the propulsors propellers ( = 30  π ⁄ ) because that is the more natural way of setting the propeller rotation.The output values are chosen to be the position coordinates and quadrotor attitude.
Based on the given and calculated values which are needed for calculating the output values, input values can be given for which the quadcopter will behave in the expected way.In all four cases the state of constant altitude is wanted to be achieved.In order to achieve that state, controllable variable  1 has to be the same value as the gravitational force.In cases 2, 3 and 4 some change in quadcopter altitude occurs because of the change of thrust force component value in the direction of the  E axis due to the change in either roll or pitch angle.Also, the model does not have implemented control algorithm which would compensate for the quadcopter attitude changes.In the first case, only  1 is different from 0. In the rest of the cases,  1 stays approximately constant, while other variables change their values.In the first case the hovering mode is simulated (Fig. 5), for which all four propellers had the same RPM.It can be seen that  1 is equal to the gravitational force.The rest of the controllable variables are equal to 0.
Also, it can be seen that the quadcopter position (Fig. 6 is showing  E coordinate) and attitude is without change, i.e. the quadcopter is in the hovering mode.
Controllable variable  2 value is very small, but still it causes the change in the quadcopter attitude. 3 and  4 are equal to 0.
Next, the second case checks quadcopter behavior due to the change of  2 value (Fig. 7).
Due to the  2 , the change in roll angle  value occurs, as well as the change in  E coordinate.Change of the  E coordinate is previously explained.
In the third case, behavior due to the change of  3 value (Fig. 10) is tested.
Just like  2 , controllable variable  3 value is also very small, but still it causes the change in the quadcopter attitude.In this case,  2 and  4 are equal to 0.
Due to the  3 , the change in pitch angle  value occurs, as well as the change in  E coordinate.
Again, controllable variable  4 value, in this case, is very small, but still it causes the change in the quadcopter attitude.In this case,  2 and  3 are equal to 0.
Due to the  4 , the change in yaw angle  value occurs, i.e. the quadcopter is rotating with respect to Z B .

CONCLUSIONS
In this article, the complexity of the quadcopter mathematical model derivation is shown.The influences of various parameters on quadcopter behavior are modelled from the aspects of physics and mathematics.This mathematical model does not include influential parameters that would complicate model without the significant improvement of the model accuracy.
Quadcopter behavior for different RPMs is shown by simulation.Simulation results show the roll, pitch and yaw angles as a function of time.As it can be seen, the function is of the 2 nd order, which is expected considering that controllable variables are a function of squared angular velocities.Also, position in ℱ B is shown as a function of time.
Further work will include mathematical model improvements by modelling the brushless DC electric motor (BLDC) and the propeller, including the effects of air as a medium, modelling the ground effect and possible disturbances.

Figure 3 .
Figure 3. Linear and angular accelerations with the respect to the ℱ B .Controllable variable  1 , i.e. propulsors thrust force   , is directed in the positive direction of Z B axis.By increasing or decreasing propellers angular velocities by the same amount, thrust force   is also increasing or decreasing, enabling the quadcopter movement along the Z B axis or hovering.If the quadcopter is on horizontal position (roll and pitch angles are equal to 0), Z axes of ℱ B and ℱ E coincide.If the quadcopter is not in horizontal position, the thrust force creates vertical and horizontal linear accelerations in ℱ E .Controllable variable  1 equation is  1 =  Z = ∑   4 =1 =  ∑   2 4 =1 =  ( 1 2 +  2 2 +  3 2 +  4 2 ),(16)
where ̇ is generalized velocity vector in ℱ B ,  is generalized velocity vector in ℱ B , and  is generalized rotation and transformation matrix. consists of quadcopter position  and

Table 1 .
Mathematical model input parameters.