1. Introduction

This paper is an extension of work originally presented in the 6th International Conference on Control, Automation and Robotics [1].

Multicopters as unmanned aerial vehicles (UAVs) are gaining massive attention in research because of their broad military and civilian applications [2]. Compared with other flying vehicles, multicopters have unique capabilities such as static hovering, vertical take-off and landing (VTOL) [3], and the ability to switch between any two arbitrary directions very quickly any time [4]. The simplest, most efficient, and most economical multicopter is the quadcopter since it contains the minimum number of propellers required to fully control its attitude and position. The quadcopter’s dynamic model is used as a basis in designing and adjusting controllers for best flight performance. This new work presents cheaper, alternative methods to model the quadcopter using only basic mechanical tools and electrical measuring devices. The simpler alternative modeling methods presented in this paper can be used to successfully operate the quadcopter in basic flight or even in dynamic fail-safe flight mode. In fact, the modeling methods in this paper were employed in another study [1] on fail-safe control if a motor failure were to occur. The complete methods presented in this work are able to capture the effect of air resistance and propeller torque on the quadcopter dynamics without requiring expensive equipment such as a reaction torque sensor. This work is thus very useful to small-scale researchers who have limited equipment but who nevertheless wish to perform research works on quadcopters.

2. Background Information

2.1. Quadcopter Kinematics

In geometry, a coordinate is a set of values that locates a point. A coordinate frame is a system that uses a minimum of  coordinates to uniquely determine the position of points of geometric objects in an -dimensional manifold such as the 3- dimensional Euclidean space. This is done by assigning an -tuple of values to every point on the manifold. One coordinate frame is linearly transformed into another through rotations and translations [5]. Thus, every rotation and translation of a coordinate frame is a transformation into another coordinate frame and because the operations involved in this transformation are linear, this transformation may be completely captured by matrix multiplication for the rotation and matrix addition for the translation.

In this paper, characters in boldface are vector or matrix quantities. Consider now the quadcopter kinematic diagram in Figure 1.a. There are three rotation angles , ,  and three axial coordinates , , . Four propellers are arranged symmetrically about the quadcopter’s center of mass (COM). Propellers 1 and 3 rotate opposite to the rotation of propellers 2 and 4. Three coordinate axes , ,  are drawn along the quadcopter arms and with origin at the COM. The direction of roll, pitch, and yaw angles and their corresponding body frame angular velocity components are indicated by rotation arrows along the coordinate axes.

Since the motion of the Earth is relatively very slow compared with the motion of the quadcopter in flight, a non-rotating Earth may be assumed. Moreover, since the radius of curvature of the Earth at any point on its surface is very large in comparison with the distances traveled by a quadcopter in flight, a flat Earth may be further assumed.

Figure 1: Quadcopter’s (a) Kinematic [1] and (b) Dynamic Model

In this paper, a superscript indicates the coordinate frame with which an entity is defined. A rotation matrix that rotates  to  will be denoted by . Note that for a rotation matrix [5],

The quadcopter state variables are listed in Table 1. To control a quadcopter, state space readings are required. Under basic flight conditions, a state space that completely specifies the six degrees of freedom of a quadcopter and may be used for dynamic computations is .

Table 1: Quadcopter State Variables [1]

A generally accepted quadcopter dynamic model is shown in Figure 1.b. The four motors and propellers are assigned a mass  centered at a distance  from the quadcopter’s geometric center such that all masses  are  apart from each other with respect to the quadcopter’s geometric center. The remaining body mass of the quadcopter excluding the four motors and propellers is  which is centered at the quadcopter’s geometric center. The centers of the four masses  and of mass  are assumed to be coplanar. The total mass of the quadcopter is thus

Since the quadcopter’s dimensions are negligibly small in comparison with the Earth’s radius of curvature, it is safe to assume that the center of gravity of any part of the quadcopter coincides with the object’s COM.

2.2. Moment of Inertia

Using the orthonormal , ,  axes in Figure 1.a, a symmetric inertia matrix may be defined as follows:

The diagonal components correspond to the axial moment of inertias (MOIs) of the quadcopter. The non-diagonal elements represent any inertial coupling that may be present between two axes of rotation. For a quadcopter that is symmetric about all three axes, the non-diagonal components of the inertia matrix are zero:

2.3. Drag Torque

Due to air resistance, the drag vector  acting about  may further be added into the model especially since the quadcopter will exhibit spinning about this axis after losing one or more propellers as implemented in this project. This drag vector is experimentally modeled in [6] as linearly increasing with the yaw angular velocity:

where  is the drag coefficient. However, this is only true when the angular speed is very low and there is no turbulence [7]. In this work, it was experimentally found out that the drag torque is better approximated by

where “ ” denotes the signum function. The approximately quadratic relationship provided by , especially for high angular speeds, is theoretically supported [8]. Meanwhile, if the quadcopter is tied to a bearing,  is an offset provided by the rotational kinetic friction after the bearing’s rotational static friction has been overcome. Without the bearing,  is simply set to zero.

2.4. Rotational Dynamics

The complete rotational dynamic model for the quadcopter is [5]

3. Experimental Methods

3.1. Hardware Design

The quadcopter’s skeletal body is the DJI Flame Wheel 450 with the commercial DJI E600 propulsion system [11]. The electronic speed controllers (ESCs) of this propulsion system work by current control and, as with most ESCs in the industry [12], vary the speed of the propellers depending on the on-time of the sent commands and not on the duty cycle. Thus, twice the duty cycle at twice the frequency would result in the same effective on-time and thus, the same propeller speed. This sets a limit on the maximum allowable frequency for the ESCs as the period of the signal sent must be greater than the minimum on-time required to turn on the propellers. The DJI E600 ESCs can receive motor signals between  and . Note that this limit on the actuation frequency also sets a limit on the effectivity of the flight controller since even if the controller loop frequency is increased, a constant delay in the actuation will persist. In this project, the target loop frequency of the controller was set at .

Figure 2: Quadcopter Platform Designed

A top view of the quadcopter for this project is shown in Figure 2. The letters N, E, W, and S indicate the designated North, East, West, and South directions respectively along the quadcopter’s body frame axes. Since cross configuration is used in this paper,  is assigned to point along the North;  is assigned to point towards the reader; and  completes a right-handed coordinate frame. Note that the quadcopter is not symmetric about the  plane nor is it symmetric about the  plane. However, the quadcopter is symmetric about the  and  where  is  rotated  about  and  is  rotated  about . An axis formed by any plane, the  plane for example, and any plane orthogonal to  that bisects the body symmetrically is an axis of symmetry while such a bisecting plane is a plane of symmetry. If more than one axis of symmetry exists, then the intersection of any two axes of symmetry will be the body’s origin. Assuming that the body’s mass is evenly distributed and placing  such that  passes through the body’s COM, the existence of one plane of symmetry means that rotation of the body about an axis created by any other plane orthogonal to  and passing through the origin will be as if the body is rotating about an axis of symmetry. This means that for this quadcopter, even though  and  are not axes of symmetry, rotating about  and  will be as if they are axes of symmetry and so Equation (4) will still hold using .

3.2. Quadcopter Modeling

Methods to obtain the propeller parameters and subsequently complete the quadcopter’s dynamic model are discussed in this section.

3.2.1. Propeller Parameters

Typically, a force sensor and a reaction torque sensor are necessary to come up with the propeller thrust-torque and thrust–pulse width modulation (PWM) relationships [13, 14]. However, these sophisticated sensors are highly expensive and not readily available. It is also extremely difficult, if not impossible, to find any force sensor or torque sensor product that will readily couple with a commercial UAV propeller to establish a working setup with enough airflow space for propeller modeling purposes. With these considerations in mind, an alternative setup was designed in this work to obtain such relationships for the model. The setup is shown in Figure 3.

The microcontroller sends PWM command to the propeller. As the propeller spins, a weighing scale measures the thrust while a tachometer measures the propeller angular speed. The reaction torque is computed from the voltage and current used to power the spinning propeller.

Figure 3: Alternative Propeller Modeling Method

The succeeding discussion on reaction torque appears in [1]. For every rotating propeller , the produced thrust magnitude , reaction torque magnitude , and angular speed  are related by [15]

where , ,  are positive constants. At steady state, electrical torque  equals mechanical torque  [16]. The propeller exerts a torque in a direction opposite its sense of rotation [17]—a part of this is the reaction torque  exerted by the spinning propeller onto the body frame while the remaining is dissipated by aerodynamic drag . Thus, at steady state,

where  is also a positive scalar. Since  and  behave similarly with respect to , it is impossible to distinguish between the two. We can set  so that

is lumped together with the drag torque  on the body frame via experimental modeling. This effectively eliminates whatever excess torque allotment that went into .

3.2.2. Axial Moments of Inertia

The axial MOIs of the quadcopter were obtained by suspending the quadcopter on three different axes as shown in Figure 4. The MOI of the quadcopter about the axis located at the pivot point is

where  is the measured period of oscillation and  is the distance between the COM and the pivot point; Equation (17) is derived in [19]. Note that  is equal to the quadcopter arm length plus possibly a small additional length from the rope that tied the quadcopter to the pivot point. Using parallel axis theorem, the MOI of the quadcopter about the body frame axis passing through the quadcopter’s COM is

To accurately measure , the quadcopter is filmed while swinging freely at 240 fps. The measured time intervals between 16 complete swings since the start of filming are averaged to yield an estimate of . This is repeated three times. The average of the three trials is the assumed value for . This averaging is important to protect the data from measurement errors. With this value, the quadcopter axial MOIs are then obtained using Equation (17) and Equation (18).

Figure 4: Measuring the Quadcopter MOI About (a) , (b) , and (c)

Meanwhile, to obtain the MOI of the motors and propeller blades about their axis of rotation, they were first weighed. Assuming that the propeller blade is a circular disk with an evenly distributed mass “m” _”b” and radius “r” _”b” , its axial MOI can be shown to equal

Assuming that the motor is a solid cylinder with an evenly distributed mass “m” _”m” and radius “r” _”m” , its axial MOI can be shown to equal

Since the propeller blade and motor have the same axis of rotation, their MOIs add up so that the propeller MOI about its axis of rotation is simply

These were the same methods used in [6] to obtain “J” _”zz” ^”P” . Since there are four propellers, four “J” _”zz” ^”P” were computed. Since the propellers are manufactured to be identical, their “J” _”zz” ^”P” values were averaged and the average value was taken as the “J” _”zz” ^”P” of any of the propellers. Again, this averaging was done since the motors and propellers are manufactured to be identical; this averaging would protect from experimental errors and remove biases.

3.2.3. Drag Coefficients

The behavior of the drag torque due to air resistance may be experimentally determined by allowing the quadcopter body to rotate about “z” ̂^”b” and applying different torques along the clockwise or counterclockwise direction. Increasing the applied torque will accelerate the rotation of the quadcopter about “z” ̂^”b” until the drag torque due to the quadcopter’s angular speed equals the applied torque, in which case the quadcopter’s rotation about “z” ̂^”b” reaches its steady state angular speed. Thus, different applied torques will correspond to different steady state angular speeds about “z” ̂^”b” and the relationship of the applied torques (which become equal to the drag torques at steady state) with respect to the measured steady state angular speeds will determine the behavior of the drag torque due to air resistance. This behavior will be described by a quadratic equation consisting of drag coefficient(s). The same curve-fitting method based on experimental data plot, albeit using a simpler linear assumed model equation, was employed in [6].
Since air’s viscosity increases with greater pressure and lower temperature, the air resistance will decrease with higher temperature and higher elevation (which lowers air pressure). It is possible to model the variation of the air resistance with respect to temperature and pressure if both the pressure and temperature of the environment can be controlled and maintained. However, without a closed room facility and equipment that maintains temperature and pressure, it is impossible to keep the temperature and pressure constant throughout the modeling of the air resistance versus angular velocity relationship. Moreover, a changing set of drag coefficients would imply a changing plant—this means that the periodic equilibrium solutions and controller decisions should also be dependent on temperature and pressure. Due to the lack of the aforementioned facility and equipment, and to simplify the problem, the air resistance versus the angular speed relationships were assumed constant and invariant to changes in pressure and temperature; this was also assumed in [6]. In a similar way, temperature and pressure could also affect the PWM to thrust relationship of the propellers. However, due again to the lack of the required facility and equipment to control ambient temperature and pressure and to simplify the problem, the PWM to thrust relationship was assumed constant and invariant to changes in pressure and temperature; this was also assumed in [13].

Figure 5: Rope and Bearing Setup for Air Resistance Characterization [1]

The drag torque on the quadcopter’s body due to the air resistance was modeled by suspending the quadcopter along a bearing with a single axis of rotation as shown Figure 5. For the air resistance characterization setup, U-shaped metal rods were welded onto a circular bearing, resulting in an axially rotating pivot point. The quadcopter is made to rotate about “z” ̂^”b” by commanding propellers 1 and 3 to exert a certain torque of the same value while disabling propellers 2 and 4. Each torque value will correspond to a steady state yaw angular velocity: the higher the torque value, the greater the magnitude of the steady state angular velocity of the quadcopter will be. The parametrization of the plot of the combined torque versus the steady state angular velocities will give the coefficients for the drag torque.
Because of metal welding, the metal junction areas expanded, making it hard for the bearing to rotate; this problem was solved using a rubbing compound to scrape off the expanded junction areas. A spherical ball bearing would allow for an omnidirectional rotation and would appear to be more preferable but under the same setup in Figure 5, the rotational friction it yielded was too high that it failed to rotate because of the weight of the quadcopter. In the end, the circular bearing was retained. Ideally, the rotational friction of the bearing would be unaffected by the mass that it carries. However, in reality, the heavier the mass, the greater the rotational kinetic and static friction. Dry lubricant was applied onto the bearing’s interior to reduce the rotational friction.

3.3. Simulation of Basic Flight Controller with the Model Derived

The proportional–integral–derivative (PID) controller output is

where ■(“u” )(“t” ) is the control output; ■(“y” )(“t” ) is the state variable to be controlled; ■(“y” _”d” )(“t” ) is the desired reference value; “k” _”p” , “k” _”d” , and “k” _”i” are non-negative constants to be tuned before running the controller. Implementing Equation (22) in a microcontroller or computer will result in a controller output that updates in discrete instances according to the controller loop frequency wherein the integral term becomes a running sum. Examining Equation (22), a sharp impulse called derivative kick can saturate the actuator and push away the system from the linear zone whenever a task requires ■(“y” _”d” )(“t” ) to change from one constant value to another [20]. Because of this, most of PID architectures use the derivative of ■(“y” )(“t” ) only. Thus,

where Equation (24) is the discretized control law version of Equation (23) that can be implemented in microcontrollers.
The simulations were done at a controller frequency of “450 Hz” ; the higher the controller frequency, the closer Equation (24) will be to Equation (23). This frequency is chosen since it is the maximum signal frequency that can be sent to the motor ESCs.
The differential reference values used in tuning the basic flight controller in the simulations were “±2 m” for the altitude, “±5°” for the roll and pitch, and “±30°” for the yaw. Furthermore, the settling period was set to “3 s” ; the desired number of oscillations for the altitude was set to “1” ; the desired number of oscillations for the yaw was set to “2” ; the desired number of oscillations for the roll and pitch was set to “11” . The greater the desired number of oscillations for a variable, the stronger its proportional controller, the higher the variable overshoots, and the faster the system response will be for the variable.
In tuning the PID gains, “k” _”i” is initially set to zero. Afterwards, “k” _”p” is manually increased until the desired number of oscillations is reached. Then “k” _”d” is manually increased until no oscillation remains. The allowable error was set to “1 cm” for the altitude and “0.1°” for the angles. Note that in this PID tuning method, a single initial overshoot in ■(“y” )(“t” ) is allowed. “k” _”i” may be introduced in the end if steady state errors are to be removed. This tuning procedure was employed for the PID controllers of “ϕ” , “θ” , “ψ” , and “z” with controller outputs “u” _”ϕ” , “u” _”θ” , “u” _”ψ” , and “u” _”z” . For the quadcopter,

where “B” is the base weight which is set to “1” /”4″ “Mg” or one-fourth the total weight of the quadcopter; this is to ensure equal sharing of load among the four motors at hover. “u” _”z” , before being used in Equation (25), is first rescaled by dividing it by “z” ̂^”b” ” ∙” “z” ̂^”i” “=cos” (“ϕ” )”cos” (“θ” ) or the directional cosine between the quadcopter’s body frame z-axis and the direction pointing opposite to gravity in order to take into account the quadcopter tilt which reduces the lift force [1].

4. Results and Discussion

4.2. Quadcopter Model and Simulation

The obtained parameters of the quadcopter used in this project are listed in Table 2. The standard acceleration due to gravity value is assumed for “g” . The rest of the parameters were obtained using methods in Section 3.
Table 2: Parameters of the Quadcopter

 

For the propeller characteristics, the following models were experimentally obtained:

where “f” _”i” is the force value in “N” , “P” _”i” is a command value from “0” to “255” with “0” corresponding to “0%” duty cycle and “255” corresponding to “100%” duty cycle, “τ” _”i” is the torque value in “N m” , and “ω” _”i” is the propeller speed in “rad” ⁄”s” of propeller “i” . Note that Equation (27) and Equation (28) are based on Equation (10) and Equation (12).
Following the alternative method designed, as described in Section 3, the propeller thrust to PWM, propeller thrust to angular speed, and propeller thrust to torque relationships were obtained. This final set of parameters was also used in computing the total torque exerted by the propellers for the drag coefficients. Figure 6, Figure 7, and Figure 8 show plots of these combinations. The plots in yellow correspond to the final model whose parameters are the average of those of the obtained models for the four propellers.

Figure 6: Propeller Thrust as a Function of PWM Command Value

 

Figure 7: Propeller Thrust as a Function of Propeller Angular Speed

 

Figure 8: Reaction Torque as a Function of Propeller Thrust

The values of “h” _”1″ , “h” _”2″ , “c” _”1″ , “g” _”1″ , and “g” _”2″ were estimated from experimental data using the setup in Section 3.2.1. In particular, “h” _”1″ and “h” _”2″ were used to convert thrust commands into motor command values. The numbers from “0” to “255” were scaled to “40000” which was the limit of the microcontroller used. It was eventually found out after tuning the fail-safe flight controller in a different work [1] that this improves quadcopter performance. This reduced the discretization of thrust commands by two orders of magnitude and resulted in a highly smoothened or less shaky flight.

Figure 9: Steady State Yaw Angular Velocities of the Quadcopter Body as a Function of Propeller Torque

Figure 9 shows a plot of steady state angular velocities of the quadcopter body about “z” ̂^”b” achieved at particular torque values exerted by propeller 1 and propeller 3 as well as a plot of the torque values predicted by the drag coefficients “γ” _”1″ and “γ” _”2″ in Table 2. Notice that the model for the drag is more accurate at higher angular velocities. This property proved useful in fail-safe mode [1] where the quadcopter’s yaw rate would exceed normal operation values.

4.2. Basic Flight Controller

Figure 10 shows the breakdown of the basic flight controller to various sub-algorithms with their execution times. In the Main Loop, the Timing part updates the time parameters for the PID controller and for the flight instructions. The SPI Angles is the part where angles-related information is requested and received in a buffer. An Adjustable Delay is inserted to control the loop frequency. The Decode Angles part converts the information stored in the buffer into angular information that can be accessed in the code. The Basic Flight Control part is where the PID control law is computed. The Flight Instructions part is where certain flight decisions are made and where the PWM commands are sent to the motors.

Figure 10: Basic Flight Control

The Main Loop is interrupted every “1” /”10″ of the time by the Request Attitude part where a function triggers the ultrasonic sensor to begin altitude measurement. Consequently, a Height Control part interrupts every “1” /”10″ of the time to calculate the

requested altitude reading and corresponding altitude PID output. The Main Loop’s total execution time was  corresponding to a  frequency while the execution time of the altitude interrupts was  which occurred at a frequency of around .

The formula for computing the average execution time of the algorithms will now be discussed. Consider algorithm 1 with execution time  interrupted by algorithm 2 with execution time  every  as illustrated in Figure 11.

Figure 11: Algorithm 1 Interrupted by Algorithm 2 Every

Then the processes can happen as demonstrated in Figure 12.

Figure 12: An Example of an Execution of the Processes in Figure 11

The average loop execution time is thus

Using Equation (29) with , , and , the average loop period is  so that the average loop frequency of the basic flight controller will be . The observed  to  main loop frequency in the implementation was due to intermittent drops in the clock frequency of the Raspberry Pi’s processor.

4.3. Basic Flight Simulation

This section documents the simulation of basic flight for the modeled quadcopter. As shown in Table 3, the PID outputs were capped based on the characteristic quantity  discussed in Section 3.3.

Table 3: Limits of Output Variables for Basic Flight

The basic flight simulations were done through MATLAB R2017a’s ode45() solver using the parameters modeled and the experimentally determined drag coefficients in Table 3. In the absence of the bearing system that was later used in the fail-safe implementation [1], the rotational kinetic friction contribution of the bearing  was set to zero. These parameters corresponded to the quadcopter implemented.

The state variables of interest are , , , and . The tuned PID gains of the quadcopter in the simulation are listed in Table 4. The second symbol in the subscript indicates the variable to control. Since for this project’s basic flight requirement, small steady state errors are of no concern, the integral gain was zeroed out for all variables to avoid unwanted system oscillations. The altitude was first tuned with the remaining variables’ PID controllers disabled. Then with only the tuned altitude controller enabled, the roll and pitch controllers were separately tuned. Finally, the yaw controller was tuned with the rest of the tuned controllers enabled.

Table 4: Tuned PID Gains for Basic Flight Simulation

Results of the simulation with the tuned PID gains are shown in Figure 13, Figure 14, Figure 15, and Figure 16. Figure 13 demonstrates the quadcopter’s stability and altitude response to ascent and landing commands. Figure 14 demonstrates the roll, pitch, yaw, and altitude stability of the modeled quadcopter as it accomplishes a positive roll command. The model was also tested for negative roll, and behaved similarly stable. Figure 15 demonstrates the quadcopter’s stable response to a negative pitch command. The model was also tested for positive pitch and behaved similarly stable. The quadcopter remains stable and achieves the desired values as expected. Finally, Figure 16 demonstrates both positive and negative yaw responses. As expected, the positive and negative responses in all variables are symmetric, and the attitude and altitude controls are capable of maintaining stability upon a command change in a single flight parameter.

Figure 13: Simulated Quadcopter Altitude Response. In (a), the quadcopter is commanded to ascend from  to . In (b), the quadcopter is commanded to perform controlled landing from  to . The quadcopter’s attitude was unaffected during ascent or landing.

Figure 14: Simulated Quadcopter Response when Banking from  to  Along the Roll Axis. In (a), the quadcopter’s roll angle, initially , smoothly achieves the  command. In (b), the variation in the quadcopter’s pitch angle with time while banking is negligibly small. In (c), the variation in the quadcopter’s yaw angle with time while banking is similarly negligible. In (d), the gradual drop in the quadcopter’s altitude while banking is also negligible.

Figure 15: Simulated Quadcopter Response when Banking from  to  Along the Pitch Axis. In (a), the variation in the quadcopter’s roll angle with time while banking is negligibly small. In (b), the quadcopter’s pitch angle, initially , smoothly achieves the  command. In (c), the variation in the quadcopter’s yaw angle with time while banking is similarly negligible. In (d), the gradual drop in the quadcopter’s altitude while banking is also negligible.

Figure 16: Simulated Quadcopter Response when Rotating  Along the Yaw Axis. In (a), the quadcopter’s yaw angle, initially , smoothly achieves the  command. In (b), the quadcopter’s yaw angle, initially , smoothly achieves the  command. The quadcopter’s roll, pitch, and altitude were unaffected during rotation.

4.4. Basic Flight Implementation

This section documents the simulation of basic flight for the modeled quadcopter.

4.4.1. Control and Actuation Frequency

As simulated, the average loop frequency of the basic flight controller was around  with the altitude controller interrupting every  of the time. This was found to work well with the ultrasonic proximity sensor having a working frequency of  in its datasheet [21]. The frequency used in the motor commands for the basic flight was  since the ESC command frequency limit was .

4.4.2. Sensor Fusion Technique

The roll and pitch of the quadcopter can be extracted from the , , and  accelerometer axial readings using [22]

while the yaw can be extracted from the , , and  magnetometer axial readings together with the estimates for the roll and the yaw using [23, 24]

Although these are absolute information about the attitude, these are very noisy estimates because of the noisy accelerometer and magnetometer readings. The gyroscope readings which are used to compute , , and  are less noisy but they can only give differential estimates for the angles. To get accurate and clean enough attitude estimates, complementary filter may be used to fuse the absolute but noisy estimates with the less noisy but differential estimates according to the weighting equation:

with  as the parameter.

The result is a low-pass filtering of the accelerometer and magnetometer’s absolute estimates to remove the noise and a high-pass filtering of the gyroscope’s differential estimates to remove the drift which worsens with each iteration. The time constant  is for both filters for a loop frequency given by  [25]. The tuning of  was done based on the value of  which was adjusted in powers of  of the average loop period  [1].

4.4.3. EMA Filter for Gyroscope

Using the smoothing factor, , an exponential moving average (EMA) filter was applied on the gyroscope readings , , and  in order to remove noise according to the formula

The time constant of an exponential moving average is the amount of time for the smoothened unit step response to reach . The relationship between this time constant, , and  is given by

where  is the iteration frequency [26].

4.4.4. Sensor Tuning

The exponential smoothing of , ,  were done first before performing the sensor fusion technique [1]. The  for each gyroscope reading was tuned while the quadcopter was being held mid-air. The tuning of  was done based on the value of , in powers of  of the average loop period  [1].  for each of the variables was doubled until the amplitude of oscillations in  and  were within . The same  was adopted for . Since the  readings of the sensor was inherently less noisy, the amplitude oscillations in  using the same  in  and  were within .

After tuning the value of , the quadcopter was made to sit still on a flat surface for  minute. The gyroscope readings along each rotation axis for the minute were averaged to get the gyroscopic axial bias which must be subtracted from the future , ,  readings of the gyroscope. The corresponding angular values from the accelerometer readings for the minute were also averaged to get the biases on the estimates of  and ; these must be subtracted from the future accelerometer estimates of  and  prior to sensor fusion. Meanwhile, the  value for the complementary filters in the roll, pitch, and yaw were increased until the amplitude of the oscillations were within . Table 5 lists the filter values.

Table 5: EMA and Complementary Filter Values

4.4.5. PID Tuning

A lot of non-idealities that appeared in the actual implementation of basic flight were not captured in the simulation. These include irregular and intermittent delays in the processor, ESCs, and actuation of motors in addition to sensor noise and possibly slight modeling discrepancies. Thus, although implementing the simulated controller gains worked for the quadcopter’s basic flight, its performance was improved by altering the simulated tuned controller gain values to better adapt to these non-idealities. Using the simulated gains in Table 4 as the starting point, the basic flight controller PID gains were retuned to the values listed in Table 6. Note that small integral gains for the roll and pitch angles, which were zeroed out in the simulation, became necessary in the implementation to ensure that the quadcopter did not drift from the reference command values over time, since drifting along the roll and pitch axes would quickly accelerate the quadcopter to dangerously high translational speeds.

Table 6: Tuned PID Gains for Basic Flight Implementation

To avoid external airflow disturbances such as wind gusts during the implementation of basic flight, the quadcopter was flown indoors. The hovering test produced the altitude, roll, pitch, and yaw response in Figure 17. In the following plots, the actual response is shown in blue whereas the reference command is shown in red. As expected, because of the aforementioned non-idealities, more jitter is present in the actual quadcopter response than in the ideal simulation results, but the settling of each parameter emphasizes the correct controller design which in turn verifies the quadcopter model derived.

Compared with the obtained altitude response using the flight controller in this work, the altitude control performance of typical modern UAVs [27] are more than twice as jittery. Moreover, the altitude tracking error of modern quadcopters using a similar ultrasonic sensor [28] can be more than 20 times worse than that in Figure 17. Although quadcopters using more sophisticated sensors like laser radar [29] can reduce altitude tracking error to almost one order of magnitude better than the altitude response presented in this paper, the controller in this work nevertheless performs more smoothly and with an almost equal fluctuation in the altitude.

Moreover, the graph responses for positive and negative banking along the roll and pitch directions are shown in Figure 18 and Figure 19. The quadcopter was set to hover at the beginning of both tests, and upon achieving stability, commanded to bank before returning to hover mode. In both plots, the quadcopter was able to find and settle around the command values in each change of state. The yaw responses to positive and negative rotation commands as shown in Figure 20 also demonstrate quick stabilization in either direction.

In terms of attitude control, the roll, pitch, and yaw responses of the quadcopter system developed in this work are typically more jittery compared with the responses of quadcopters using popular firmware systems like ArduPilot [30] and Pixhawk [31, 32] but the developed system can nonetheless achieve a lower average tracking error in the roll and pitch angles.

5. Summary and Conclusion

The study came up with alternative methods to determine thrust vs. PWM command, thrust vs. angular speed, and reaction torque vs. thrust relationships of the quadcopter’s propellers as well as the coefficients of the aerodynamic drag acting on the quadcopter, thereby completing a dynamic model for the quadcopter. With the specific model, the study extended into designing and simulating a basic flight controller. The simulation led to determining PID values and limits which served as starting points for the actual controller implemented. Plots of the simulated responses to thorough actual tests demonstrated the expected stable flight of the quadcopter. These results further verified the completeness and correctness of the model obtained. This method was repeated for multiple versions of quadcopter with a completely different set of motors, ESCs, and propellers. Using the same methodology, corresponding values were obtained and similar behavior was demonstrated.

The quadcopter implemented in this paper was used in another work [1] which is the implementation of a fail-safe controller for a quadcopter that experiences a motor failure. A working basic flight controller is a stepping stone to solving the motor failure problem.

Conflict of Interest

The authors declare no conflict of interest.

Acknowledgment

This project was supported in part by grants from the UP Engineering Research & Development Foundation, Inc. Publication of this paper was made possible by the Engineering Research and Development for Technology program of the Department of Science and Technology in the Philippines.

Figure 17: Quadcopter Hover Response. In (a), the quadcopter’s altitude remains stably close to the  reference value with time. In (b), the quadcopter’s roll angle remains stably close to the  reference value with time. In (c), the quadcopter’s pitch angle also remains stably close to the  reference value with time. In (d), the quadcopter’s yaw angle remains similarly stable and close to the  reference value with time.

Figure 18: Quadcopter Roll Response to Positive and Negative Banking Commands. In (a), the quadcopter’s roll angle, initially , achieves a  command for  before returning to . In (b), the quadcopter’s roll angle, initially , achieves a  command for  before returning to . The quadcopter’s pitch, yaw, and altitude were not significantly affected during banking.

Figure 19: Quadcopter Pitch Response to Positive and Negative Banking Commands. In (a), the quadcopter’s pitch angle, initially , achieves a  command for  before returning to . In (b), the quadcopter’s pitch angle, initially , achieves a  command for  before returning to . The quadcopter’s roll, yaw, and altitude were not significantly affected during banking.

Figure 20: Quadcopter Yaw Response to Positive and Negative Rotation Commands. In (a), the quadcopter’s yaw angle, initially , achieves the  command. In (b), the quadcopter’s yaw angle, initially , achieves the  command. The quadcopter’s roll, pitch, and altitude were not significantly affected during rotation.

1. G. P. S. Rible, N. A. A. Arriola, M. C. Ramos, “Fail-Safe Controller Architectures for Quadcopter with Motor Failures” in 2020 6th International Conference on Control, Automation and Robotics (ICCAR), Singapore, 2020. https://doi.org/10.1109/ICCAR49639.2020.9108038
2. B. Fan, J. Sun, Y. Yu, “A LQR controller for a quadrotor: Design and experiment” in 2016 31st Youth Academic Annual Conference of Chinese Association of Automation (YAC), Wuhan, Hubei, China, 2016. https://doi.org/10.1109/YAC.2016.7804869
3. J. Qiao, Z. Liu, Y. Zhang, “Gain scheduling PID control of the quad-rotor helicopter” in 2017 IEEE International Conference on Unmanned Systems (ICUS), Beijing, China, 2017. https://doi.org/10.1109/ICUS.2017.8278414
4. Y. Sen, W. Zhongsheng, “Quad-Rotor UAV Control Method Based on PID Control Law” in 2017 International Conference on Computer Network, Electronic and Automation (ICCNEA), Xi’an, Shaanxi, China, 2017. https://doi.org/10.1109/ICCNEA.2017.24
5. R. W. Beard, “Quadrotor Dynamics and Control,” Tech. Rep., Brigham Young University, 2008.
6. M. W. Mueller, R. D’Andrea, “Stability and control of a quadrocopter despite the complete loss of one, two, or three propellers” in 2014 IEEE International Conference on Robotics and Automation (ICRA), Hong Kong, China, 2014. https://doi.org/10.1109/ICRA.2014.6906588
7. G. K. Batchelor, An Introduction to Fluid Dynamics, Cambridge University Press, 2000.
8. Y. Nakayama, R. Boucher, “9 – Drag and lift” in Introduction to Fluid Mechanics, Butterworth-Heinemann, 1998.
9. A. Nagaty, S. Saeedi, C. Thibault, M. Seto, H. Li, “Control and Navigation Framework for Quadrotor Helicopters” J. Intell. Robot. Syst., 70(1-4), 1–12, 2013. https://doi.org/10.1007/s10846-012-9789-z
10. H. t. M. N. ElkHoly, “Dynamic Modeling and Control of a Quadrotor Using Linear and Nonlinear Approaches,” MS Thesis, The American University in Cairo, 2014.
11. DJI E600 Multirotor Propulsion System, User Manual V1.04, DJI Technology Co., Ltd., 2014.
12. M.-G. Yoon, “On Driving Signal of Electronic Speed Controller for Small Multi-Rotor Helicopter” Int. J. Eng. Res. Technol., 4(11), 2278–0181, 2015. https://www.ijert.org/on-driving-signal-of-electronic-speed-controller-for-small-multi-rotor-helicopter
13. D. M. Filatov, A. V. Devyatkin, A. I. Friedrich, “Quadrotor parameters identification and control system design” in 2017 IEEE Conference of Russian Young Researchers in Electrical and Electronic Engineering (EIConRus), St. Petersburg, Russia, 2017. https://doi.org/10.1109/EIConRus.2017.7910684
14. D. M. Filatov, A. I. Friedrich, A. V. Devyatkin, “Parameters identification of thrust generation subsystem for small unmanned aerial vehicles” in 2017 IEEE II International Conference on Control in Technical Systems (CTS), St. Petersburg, Russia, 2017. https://doi.org/10.1109/CTSYS.2017.8109574
15. T. Bresciani, “Modelling, Identification and Control of a Quadrotor Helicopter,” MS Thesis, Lund University, 2008.
16. A. Hughes, Electric Motors and Drives: Fundamentals, Types and Applications, 3rd ed., Elsevier Science/Newnes, 2005.
17. H. Young, R. Freedman, University Physics with Modern Physics, 14th (Global) ed., Pearson Education, 2016.
18. Y. Nakayama, R. Boucher, Introduction to Fluid Mechanics, Butterworth-Heinemann, 1998.
19. J. Shigley, J. Uicker, Theory of Machines and Mechanisms, International ed., McGraw-Hill, 1980.
20. A. H. Ahmed, A. N. Ouda, A. M. Kamel, Y. Z. Elhalwagy, “Attitude stabilization and altitude control of quadrotor” in 2016 12th International Computer Engineering Conference (ICENCO), Cairo, Egypt, 2016. https://doi.org/10.1109/ICENCO.2016.7856456
21. Ultrasonic Ranging Module HC – SR04, User’s Manual, ElecFreaks, n. d.
22. M. Pedley, “Tilt Sensing Using a Three-Axis Accelerometer,” Freescale Semiconductor Application Note, 1, 2012–2013, NXP Semiconductors N.V., 2013.
23. H. Vathsangam, “Magnetometer – My IMU estimation experience,” 2010. https://sites.google.com/site/myimuestimationexperience/sensors/magnetometer
24. “3-Axis Compass Sensor Set HMC1055” in Sensing Earth’s magnetic field: Compassing, Magnetometry and Dead Reckoning Solutions, Magnetic Sensors Product Catalog, Honeywell Solid State Electronics Center, 2001. https://stevenengineering.com/Tech_Support/PDFs/31MAGSEN.pdf
25. S. Colton, “The Balance Filter,” White Paper, 2007. https://d1.amobbs.com/bbs_upload782111/files_44/ourdev_665531S2JZG6.pdf
26. G. Stanley, “Exponential Filter,” 2010-2020. https://gregstanleyandassociates.com/whitepapers/FaultDiagnosis/Filtering/Exponential-Filter/exponential-filter.htm
27. W. Liu, C. Yu, X. Wang, Y. Zhang, Y. Yu, “The Altitude Hold Algorithm of UAV Based on Millimeter Wave Radar Sensors” in 2017 9th International Conference on Intelligent Human-Machine Systems and Cybernetics (IHMSC), Hangzhou, Zhejiang, China, 2017. https://doi.org/10.1109/IHMSC.2017.106
28. D. J. Esteves, A. Moutinho, J. R. Azinheira, “Stabilization and Altitude Control of an Indoor Low-Cost Quadrotor: Design and Experimental Results” in 2015 IEEE International Conference on Autonomous Robot Systems and Competitions, Vila Real, Vila Real District, Portugal, 2015. https://doi.org/10.1109/ICARSC.2015.30
29. J. Zhao, Y. Li, D. Hu, Z. Pei, “Design on altitude control system of quad rotor based on laser radar” in 2016 IEEE International Conference on Aircraft Utility Systems (AUS), Beijing, China, 2016. https://doi.org/10.1109/AUS.2016.7748029
30. Ö. Elbir, A. U. Batmaz, C. Kasnakoğlu, “Improving quadrotor 3-axes stabilization results using empirical results and system identification” in 2013 9th Asian Control Conference (ASCC), Istanbul, Turkey, 2013. https://doi.org/10.1109/ASCC.2013.6606281
31. E. A. Niit, W. J. Smit, “Integration of model reference adaptive control (MRAC) with PX4 firmware for quadcopters” in 2017 24th International Conference on Mechatronics and Machine Vision in Practice (M2VIP), Auckland, New Zealand, 2017. https://doi.org/10.1109/M2VIP.2017.8211479
32. S. Sakulthong, S. Tantrairatn, W. Saengphet, “Frequency Response System Identification and Flight Controller Tuning for Quadcopter UAV” in 2018 Third International Conference on Engineering Science and Innovative Technology (ESIT), North Bangkok, Thailand, 2018. https://doi.org/10.1109/ESIT.2018.8665114

(Click Here)

(Click Here)