Modeling and Implementation of Quadcopter Autonomous Flight Based on Alternative Methods to Determine Propeller Parameters

To properly simulate and implement a quadcopter flight control for intended load and flight conditions, the quadcopter model must have parameters on various relationships including propeller thrust-torque, thrust-PWM, and thrust--angular speed to a certain level of accuracy. Thrust-torque modeling requires an expensive reaction torque measurement sensor. In the absence of sophisticated equipment, the study comes up with alternative methods to complete the quadcopter model. The study also presents a method of modeling the rotational aerodynamic drag on the quadcopter. Although the resulting model of the reaction torque generated by the quadcopter's propellers and the model of the drag torque acting on the quadcopter body that are derived using the methods in this study may not yield the true values of these quantities, the experimental modeling techniques presented in this work ensure that the derived dynamic model for the quadcopter will nevertheless behave identically with the true model for the quadcopter. The derived dynamic model is validated by basic flight controller simulation and actual flight implementation. The model is used as basis for a quadcopter design, which eventually is used for test purposes of basic flight control. This study serves as a baseline for fail-safe control of a quadcopter experiencing an unexpected motor failure.


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 failsafe 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.

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 n coordinates to uniquely determine the position of points of geometric objects in an n-dimensional manifold such as the 3dimensional Euclidean space. This is done by assigning an n-tuple of values to every point on the manifold. One coordinate frame is ASTESJ ISSN: 2415-6698 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 x, y, z. 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 x, y, z 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. In this paper, a superscript indicates the coordinate frame with which an entity is defined. A rotation matrix that rotates C a to C b will be denoted by R a b . 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 (ϕ, θ, ψ, p, q, r, x, y, z, x, y, z).  [1] x inertial position along A generally accepted quadcopter dynamic model is shown in Figure 1.b. The four motors and propellers are assigned a mass m centered at a distance l from the quadcopter's geometric center such that all masses m are 90° 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 M B which is centered at the quadcopter's geometric center. The centers of the four masses m and of mass M B 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.

Moment of Inertia
Using the orthonormal x, y, z 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:

Drag Torque
Due to air resistance, the drag vector τ d ≜τ d b =(0, 0, τ d ) acting about ẑb 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 " sgn " denotes the signum function. The approximately quadratic relationship provided by γ 1 , especially for high angular speeds, is theoretically supported [8]. Meanwhile, if the quadcopter is tied to a bearing, γ 2 is an offset provided by the rotational kinetic friction after the bearing's rotational static friction has been overcome. Without the bearing, γ 2 is simply set to zero.

Rotational Dynamics
The complete rotational dynamic model for the quadcopter is [5] which in simplified matrix form is where the algebraic sum of the propeller angular speeds is denoted by [9,10] Ω=

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 30 Hz and 450 Hz. 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 450 Hz.  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, x b is assigned to point along the North; z b is assigned to point towards the reader; and y b completes a right-handed coordinate frame. Note that the quadcopter is not symmetric about the x b z b plane nor is it symmetric about the y b z b plane. However, the quadcopter is symmetric about the An axis formed by any plane, the x b y b plane for example, and any plane orthogonal to x b y b 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 x b y b such that x b y b 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 x b y b 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 x b and y b are not axes of symmetry, rotating about x b and y b will be as if they are axes of symmetry and so Equation (4) will still hold using C b .

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

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. The succeeding discussion on reaction torque appears in [1]. For every rotating propeller i, the produced thrust magnitude f i , reaction torque magnitude τ i , and angular speed ω i are related by [15] f i =k f ω i 2 (10) where k f , k τ , k are positive constants. At steady state, electrical torque τ e equals mechanical torque τ m [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 τ d prop . Thus, at steady state, From [15], where k τ is a positive scalar. It was also shown in [18] that where k d is also a positive scalar. Since τ ψ and τ d prop behave similarly with respect to ω, it is impossible to distinguish between the two. We can set τ d prop =0 so that τ ψ =τ e (16) τ d prop is lumped together with the drag torque τ d on the body frame via experimental modeling. This effectively eliminates whatever excess torque allotment that went into τ ψ .

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 T is the measured period of oscillation and r is the distance between the COM and the pivot point; Equation (17) is derived in [19]. Note that r 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 T , 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 T. This is repeated three times. The average of the three trials is the assumed value for T. 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). 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.

Drag Coefficients
The behavior of the drag torque due to air resistance may be experimentally determined by allowing the quadcopter body to rotate about ẑb and applying different torques along the clockwise or counterclockwise direction. Increasing the applied torque will accelerate the rotation of the quadcopter about ẑb until the drag torque due to the quadcopter's angular speed equals the applied torque, in which case the quadcopter's rotation about ẑb reaches its steady state angular speed. Thus, different applied torques will correspond to different steady state angular speeds about ẑ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]. 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 ẑ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.

Simulation of Basic Flight Controller with the Model Derived
The proportional-integral-derivative (PID) controller output is u(t)=k p :y d (t)-y(t); +k d d dt :y d (t)-y(t); +k i ∫ :y d (t)-y(t); dt (22) 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 ẑb •ẑ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].

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. 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 A 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.  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.  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.   The formula for computing the average execution time of the algorithms will now be discussed. Consider algorithm 1 with execution time T 1 interrupted by algorithm 2 with execution time T 2 every t 0 as illustrated in Figure 11.  The average loop execution time is thus

Basic Flight Controller
Using Equation (29) with t 0 = 1 45 Hz , T 1 =2068 μs, and T 2 =26 μs, the average loop period is 2070.23899911 μs so that the average loop frequency of the basic flight controller will be 483 Hz. The observed 450 Hz to 470 Hz main loop frequency in the implementation was due to intermittent drops in the clock frequency of the Raspberry Pi's processor.

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 B discussed in Section 3.3.
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 γ 2 was set to zero. These parameters corresponded to the quadcopter implemented.
The state variables of interest are ϕ, θ, ψ, and z. 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. 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 0 m to 2 m. In (b), the quadcopter is commanded to perform controlled landing from 2 m to 0 m. The quadcopter's attitude was unaffected during ascent or landing. Figure 14: Simulated Quadcopter Response when Banking from 0° to 5° Along the Roll Axis. In (a), the quadcopter's roll angle, initially 0°, smoothly achieves the 5° 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 0° to 5° 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 0°, smoothly achieves the -5° 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 ±30° Along the Yaw Axis. In (a), the quadcopter's yaw angle, initially 0°, smoothly achieves the +30° command. In (b), the quadcopter's yaw angle, initially 0°, smoothly achieves the -30° command. The quadcopter's roll, pitch, and altitude were unaffected during rotation.

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

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

Sensor Fusion Technique
The roll and pitch of the quadcopter can be extracted from the x, y, and z accelerometer axial readings using [22] ϕ=arctan : while the yaw can be extracted from the x, y, and z magnetometer axial readings together with the estimates for the roll and the yaw using [23,24] ψ=arctan D magnet z sin(ϕ)-magnet y cos(ϕ) magnet x cos(θ)+magnet y sin(θ)sin(ϕ)+magnet z sin(θ)cos(ϕ)

E (32)
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 highpass filtering of the gyroscope's differential estimates to remove the drift which worsens with each iteration. The time constant τ= αdt 1-α is for both filters for a loop frequency given by dt [25]. The tuning of α was done based on the value of τ which was adjusted in powers of 2 of the average loop period dt [1].

EMA Filter for Gyroscope
Using the smoothing factor, α EMA , an exponential moving average (EMA) filter was applied on the gyroscope readings p, q, and r 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 1-1 e A ≈63.2%. The relationship between this time constant, τ EMA , and α EMA is given by where dt is the iteration frequency [26].

Sensor Tuning
The exponential smoothing of p, q, r were done first before performing the sensor fusion technique [1]. The α EMA for each gyroscope reading was tuned while the quadcopter was being held mid-air. The tuning of α EMA was done based on the value of τ EMA , in powers of 2 of the average loop period dt [1]. τ EMA for each of the variables was doubled until the amplitude of oscillations in p and q were within 0.01 rad s . The same α EMA was adopted for r. Since the r readings of the sensor was inherently less noisy, the amplitude oscillations in r using the same α EMA in p and q were within 0.006 rad s . After tuning the value of α EMA , the quadcopter was made to sit still on a flat surface for 1 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 p, q, r 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 ±0.1°. Table 5 lists the filter values.

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. 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 nonidealities, 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.

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. Figure 17: Quadcopter Hover Response. In (a), the quadcopter's altitude remains stably close to the 1.5 m reference value with time. In (b), the quadcopter's roll angle remains stably close to the 0° reference value with time. In (c), the quadcopter's pitch angle also remains stably close to the 0° reference value with time. In (d), the quadcopter's yaw angle remains similarly stable and close to the 0° reference value with time. Figure 18: Quadcopter Roll Response to Positive and Negative Banking Commands. In (a), the quadcopter's roll angle, initially 0°, achieves a +5° command for 3 s before returning to 0°. In (b), the quadcopter's roll angle, initially 0°, achieves a -5° command for 5 s before returning to 0°. 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 0°, achieves a +5° command for 3 s before returning to 0°. In (b), the quadcopter's pitch angle, initially 0°, achieves a -5° command for 5 s before returning to 0°. 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 0°, achieves the +30° command. In (b), the quadcopter's yaw angle, initially 0°, achieves the -30° command. The quadcopter's roll, pitch, and altitude were not significantly affected during rotation.