Modeling and Simulation of a Quadcopter

August 2024 - December 2024

GitHub Repository

Introduction

My fourth major course at USNA was EW301 (System Modeling and Simulation). We learned the fundamentals of developing, simulating, and validating models for dynamic systems.

The final project for this course was to research a technical paper from an academic journal or conference proceeding and create a figure from the paper.

Final Report

Abstract — Flying robots such as quadcopters can perform many tasks efficiently and safely across a variety of industries such as agriculture, construction, photography, the military, and more. With current advances in sensor, actuator, and energy storage technology, quadcopters can operate on a larger scale and tackle increasingly complex challenges. As a result, it is vital to understand how to model and simulate these complex systems for development purposes. The paper studied provides a way of modeling and simulating the movement of quadcopters with PD control. The following report delves into further analysis of the underlying modeling equations and utilizes a Euler’s method scheme to compose a simulation capable of reproducing the results found in the paper studied.

Introduction

Unmanned Aerial Vehicles, or UAVs, are robots that have the ability to fly. This gives them an upper hand over land-based robots as UAVs are capable of navigating around difficult terrain and move at faster speeds. As a result, these systems have proven to be effective at aiding a variety of industries. Some examples include infrastructure inspections, exploring dangerous regions, spraying crops, and more. With the current advancements in electronic technology and control theory, UAVs have become relatively simple to manufacture at scale. Being able to model UAVs, such as quadcopters, is essential during the developmental process. With a model, simulations can evaluate the performance of different controller designs, minimizing the risk of damage to physical components.

In [1], the authors propose a nonlinear dynamic model of a quadcopter and simulate its performance via a PD controller. The equations of motion for the six degrees of freedom were provided, along with the parameters used in simulation. In the following report, a manual reconstruction of the referenced paper was conducted supplemented by an in-depth analysis on the equations of motion and the implementation of a custom Euler integration scheme. Additionally, the report offers criticism to the control design provided and explains why it lacks practical feasibility. 

Problem Statement and Objectives

Quadcopters offer a complex modelling problem as they involve four motors connected to a rigid body with some means to move about 3D space. It is a classic underactuated system with only throttle, roll, pitch, and yaw inputs to maneuver throughout six degrees of freedom. It is up to the engineer to come up with some method of converting the individual speeds of each motor into a position and orientation of the quadcopter. A controller must then be designed to manipulate the four motor speeds to reach a desired state. 

The paper studied documents a nonlinear mathematical model of a quadcopter with three translational ($x, y, z$) and three rotational ($\varphi, \theta, \psi$) degrees of freedom (DOF). Euler angles are used to model rotation as they are more intuitive. However, they are prone to “gimbal lock” which loses one degree of freedom. This issue can be resolved by using quaternions, but doing so can be confusing for the reader as they live in 4-dimensional space. Additionally, the authors propose a PD controller to achieve altitude, roll, pitch, and yaw angle setpoints. Details on the applied simulation method were excluded from the paper, but plots of the 4 setpoint responses were provided along with the associated gain values.

Modeling

A quadcopter is defined as a multirotor with four motors, typically arranged in a “+” or “x” configuration. We decided to model the “+” variant of the quadcopter rather than the “x” variant as it is more intuitive. The front and back motors, $M_2$ and $M_4$, rotate in the counterclockwise direction whereas the left and right motors, $M_1$ and $M_3$, rotate in the clockwise direction as depicted in Figure 1. This configuration ensures that the torque due to the clockwise spinning motors counteract the torque due counterclockwise spinning motors.

Figure 1 Body Frame of Quadcopter

When modeling with Euler angles, it is crucial to maintain a standard rotation order as each rotation is relative to the new reference frame due to the previous rotation. For example, to yaw 90 degrees the quadcopter rotates 90 degrees about the original z axis. This establishes a new reference frame for the next rotation. Consequently, to roll 45 degrees, the quadcopter rotates 45 degrees about the new x axis. The new x axis due to the yaw rotation should not be confused with the original x axis before the yaw rotation. This paper will assume the order: yaw, pitch, roll.

Control of the quadcopter involves varying the speed of the motors. Rotating the motor, and thereby the propeller, generates two main forces: thrust (1) and torque (2).

$F_i = b\omega_i^2$

(1)

$T_i=d\omega_i^2$

(2)

Where $F$ is the thrust produced, $T$ is the torque, $i$ references the motor number, $l$ is the length between the motor and center of gravity, and $\omega$ is the angular velocity of the motor. The constants $b$ and $d$ are a simplified representation of the parameters affecting the thrust and torque produced by the propeller. This value can be estimated empirically and is influenced by factors such as air density and size of the propeller. Thrust will produce a lifting force as denoted in Figure 1 and torque will produce a moment in the direction opposite of the propeller rotation. Both of these phenomenon’s are derived from Newton’s third law which states, that every action has an equal and opposite reaction. The thrust generated by the propellers will push air down, causing the quadcopter to lift in the opposite direction. Additionally, the act of rotating a propeller clockwise will in turn apply a force to rotate the body counterclockwise.

In order to maneuver the quadcopter, there are four independent control inputs: throttle, roll, pitch, and yaw. To change the throttle, $U_1$, all four motors are increased or decreased equivalently (3).

$U_1=F_{tot}=F_1+F_2+F_3+F_4=b(\omega_1^2+\omega_2^2+\omega_3^2+\omega_4^2)$

(3)

To make the quadcopter roll, $U_2$, in the positive $\varphi$ direction, keep $F_1$ and $F_3$ the same, but increase $F_4$ and decrease $F_2$ (4). To roll the opposite direction, decrease $F_4$ and increase $F_2$.

$U_2=T_\varphi=l(F_4-F_2)=lb(\omega_4^2-\omega_2^2)$

(4)

To make the quadcopter pitch, $U_3$, in the positive $\theta$ direction, keep $F_2$ and $F_4$ the same, but increase $F_1$ and decrease $F_3$ (5). To pitch in the opposite direction, decrease $F_1$ and increase $F_3$.

$U_3=T_\theta=l(F_1-F_3)=lb(\omega_1^2-\omega_3^2)$

(5)

To make the quadcopter yaw, $U_4$, in the positive $\psi$ direction, increase both $F_2$ and $F_4$, but decrease both $F_1$ and $F_3$ (6). To yaw in the opposite direction, decrease both $F_2$ and $F_4$, but increase both $F_1$ and $F_3$.

$U_4=T_\psi=-T_1+T_2-T_3+T_4=d(-\omega_1^2+\omega_2^2-\omega_3^2+\omega_4^2)$

(6)

Equation (7) represents (3), (4), (5), and (6) in matrix form.

$\begin{bmatrix}U_1\\U_2\\U_3\\U_4\end{bmatrix}=\begin{bmatrix}b&b&b&b\\0&-lb&0&lb\\lb&0&-lb&0\\-d&d&-d&d\end{bmatrix}\begin{bmatrix}\omega_1^2\\\omega_2^2\\\omega_3^2\\\omega_4^2\end{bmatrix}$

(7)

From a practical control perspective, it is more beneficial to calculate the resulting motor speeds based on the four control inputs (8).

$\begin{bmatrix}\omega_1^2\\\omega_2^2\\\omega_3^2\\\omega_4^2\end{bmatrix}=\begin{bmatrix}\frac{1}{4b}&0&\frac{1}{2bl}&\frac{1}{4d}\\\\\frac{1}{4b}&-\frac{1}{2bl}&0&\frac{1}{4d}\\\\-\frac{1}{4b}&0&-\frac{1}{2bl}&-\frac{1}{4d}\\\\\frac{1}{4b}&\frac{1}{2bl}&0&\frac{1}{4d}\end{bmatrix}\begin{bmatrix}U_1\\U_2\\U_3\\U_4\end{bmatrix}$

(8)

With (7) and (8), understanding the relationship between motor speed and change in translation or rotation is straight forward. It is important to notice that the four control inputs $U_1$, $U_2$, $U_3$, and $U_4$ affect movements relative to the body. These changes however, are not equal to the movements which are observed from a fixed point of view. This problem is exacerbated when attempting to integrate sensors such as a barometer and GPS. For example, when using these sensors to estimate altitude, a fixed reference frame is being used. Altitude is defined as the location along the z axis relative to Earth. However, if the quadcopter is pitched or rolled, the z axis relative to the body no longer aligns with the z axis relative to Earth. This can cause problems when trying to maintain altitude as maintaining position on the z axis relative to the body, will not always equate to maintaining altitude or position on the z axis relative to Earth. As a result, it is important to define two coordinate systems, one body frame and one Earth frame. The body frame is represented in Figure 1 where the reference frame is always fixed to the body of the quadcopter. No matter the orientation, the body frame’s z axis will always be orthogonal to the quadcopter’s body. Furthermore, the Earth frame is fixed, with the z axis orthogonal to the surface of the Earth.

Moving between the two coordinate systems can be achieved using a rotation matrix accounting for the three successive Euler angle rotations. Before the first rotation, the body frame and Earth frame are aligned. Sticking with the yaw, pitch, roll order convention, first rotate about the z axis (yaw) through an angle $\psi$ (9). The frame $e$ refers to the original Earth frame and $e_1$ refers to the rotated reference frame.

$R_e^{e_1}(\psi)=\begin{bmatrix}\cos(\psi)&\sin(\psi)&0\\-\sin(\psi)&\cos(\psi)&0\\0&0&1\end{bmatrix}$

(9)

Next, pitch about the new y axis by the angle $\theta$ (10).

$R_{e_1}^{e_2}(\theta)=\begin{bmatrix}\cos(\theta)&0&-\sin(\theta)\\0&1&0\\\sin(\theta)&0&\cos(\theta)\end{bmatrix}$

(10)

Finally, roll about the new x axis by the angle $\varphi$ (11). The frame $b$ refers to the final body frame.

$R_{e_2}^b(\varphi)=\begin{bmatrix}1&0&0\\0&\cos(\varphi)&\sin(\varphi)\\0&-\sin(\varphi)&\cos(\varphi)\end{bmatrix}$

(11)

Equation (12) represents the derivation of the rotation matrix of the original Earth frame, $e$, in the fully transformed body frame, $b$.

$R_e^b(\varphi,\theta,\psi)=R_{e_2}^b(\varphi)R_{e_1}^{e_2}(\theta)R_e^{e_1}(\psi)$

(12)

Equation (13) is the expanded form of $R_e^b$ where $\text{s}$ represents $\sin$ and $\text{c}$ represents $\cos$.

$R_e^b=\begin{bmatrix}\text{c}(\theta)\text{c}(\psi)&\text{c}(\theta)\text{s}(\psi)&-\text{s}(\theta)\\\text{s}(\theta)\text{s}(\varphi)\text{c}(\psi)-\text{s}(\psi)\text{c}(\varphi)&\text{s}(\psi)\text{s}(\theta)\text{s}(\varphi)+\text{c}(\psi)\text{c}(\varphi)&\text{s}(\varphi)\text{c}(\theta)\\\text{s}(\theta)\text{c}(\varphi)\text{c}(\psi)+\text{s}(\psi)\text{s}(\varphi)&\text{s}(\psi)\text{s}(\theta)\text{c}(\varphi)-\text{c}(\psi)\text{s}(\varphi)&\text{c}(\varphi)\text{c}(\theta)\end{bmatrix}$

(13)

With a method of converting between both reference frames, deriving the translational equations of motions is straight forward. Newton’s 2nd law of motion states that the acceleration of an object is directly proportional to the net force acting on it and inversely proportional to its mass (14).

$F=ma$

(14)

For simplification purposes, the model assumes no drag and claims that the gyroscopic effects due to the spinning motors are negligible. Therefore, the only forces acting on the body is the total thrust, $U_1$, and gravity, $mg$ (15). Notice how the accelerations $\ddot{x}$, $\ddot{y}$, $\ddot{z}$ and $g$ are in relation to the earth frame but $U_1$ is relative to the body frame. The transformation of $U_1$ to the earth frame is achieved by multiplying the rotation matrix $R_e^b$.

$m\begin{bmatrix}\ddot{x}\\\ddot{y}\\\ddot{z}\end{bmatrix}=\begin{bmatrix}0\\0\\-mg\end{bmatrix}+R_e^b\begin{bmatrix}0\\0\\U_1\end{bmatrix}$

(15)

Expanding (15) and isolating the linear accelerations $\ddot{x}$, $\ddot{y}$, and $\ddot{z}$ is illustrated in (16).

$\begin{bmatrix}\ddot{x}\\\ddot{y}\\\ddot{z}\end{bmatrix}=\begin{bmatrix}(\sin(\theta)\cos(\varphi)\cos(\psi)+\sin(\psi)\sin(\varphi))\frac{U_1}{m}\\\\(\sin(\psi)\sin(\theta)\cos(\varphi)-\cos(\psi)\sin(\varphi))\frac{U_1}{m}\\\\(\cos(\varphi)\cos(\theta))\frac{U_1}{m}-g\end{bmatrix}$

(16)

Deriving the rotational equations of motion is a little more difficult due to the fact that the time rate of change of Euler angles does not directly equate to angular velocity of the body frame. This makes sense because the yaw angle, applied first, is not rotated about the final body yaw axis. Instead, pitch and roll rotations have to be incorporated before calculating the time rate change for yaw relative to the final body frame. Equation (17) illustrates this concept in mathematical form.

$\begin{bmatrix}p\\q\\r\end{bmatrix}=R_{e_2}^b(\varphi)R_{e_1}^{e_2}(\theta)\begin{bmatrix}0\\0\\\dot{\psi}\end{bmatrix}+R_e^{e_1}(\varphi)\begin{bmatrix}0\\\dot{\theta}\\0\end{bmatrix}+\begin{bmatrix}\dot{\varphi}\\0\\0\end{bmatrix}$

(17)

The vector, $[p\;q\;r]^T$, denotes the angular velocity of the body frame and $\dot{\varphi}$, $\dot{\theta}$, and $\dot{\psi}$ denote the instantaneous change of Euler angles. Notice how the time rate change vector for yaw, $[0\;0\;\dot{\psi}]^T$, is multiplied by the rotation matrices $R_{e_2}^b(\varphi)$ and $R_{e_1}^{e_2}(\theta)$ to account for the roll and pitch rotations.

Equation (18) is the expanded form of (17).

$\begin{bmatrix}p\\q\\r\end{bmatrix}=\begin{bmatrix}1&0&-\sin(\theta)\\0&\cos(\varphi)&\sin(\varphi)\cos(\theta)\\0&-\sin(\varphi)&\cos(\varphi)\cos(\theta)\end{bmatrix}\begin{bmatrix}\dot{\varphi}\\\dot{\theta}\\\dot{\psi}\end{bmatrix}$

(18)

It is more beneficial to calculate the time rate change of Euler angles from the angular velocity of the body frame (19). Numerically integrating $\dot{\varphi}$, $\dot{\theta}$, $\dot{\psi}$ will get the Euler rotation angles $\varphi$, $\theta$, $\psi$ to be used in (16).

$\begin{bmatrix}\dot{\varphi}\\\dot{\theta}\\\dot{\psi}\end{bmatrix}=\begin{bmatrix}1&\sin(\varphi)\tan(\theta)&\cos(\varphi)\tan(\theta)\\0&\cos(\varphi)&-\sin(\varphi)\\0&\sin(\varphi)\sec(\theta)&\cos(\varphi)\sec(\theta)\end{bmatrix}\begin{bmatrix}p\\q\\r\end{bmatrix}$

(19)

A method of calculating the angular velocity can be derived using Euler’s rotation equations. Equation (20) denotes the general vector form.

$M=I\dot{\omega}+\omega\times(I\omega)$

(20)

$M$ refers to the vector of the applied torques (21). $T_\varphi$, $T_\theta$, and $T_\psi$ denote the torques about the x, y, and z axis. Notice how these torques also equate to $U_2$, $U_3$, and $U_4$ as defined in (4), (5), and (6).

$M=\begin{bmatrix}T_\varphi\\T_\theta\\T_\psi\end{bmatrix}=\begin{bmatrix}U_2\\U_3\\U_4\end{bmatrix}$

(21)

$I$ refers to the principle inertia matrix where $I_x$, $I_y$, and $I_z$ denotes the moment of inertia about the x-axis, y-axis, and z-axis respectively (22).

$I=\begin{bmatrix}I_x&0&0\\0&I_y&0\\0&0&I_z\end{bmatrix}$

(22)

$\omega$ denotes the angular velocity of the body frame, $[p\;q\;r]^T$. $\dot{\omega}$ denotes the angular acceleration of the body frame $[\dot{p}\;\dot{q}\;\dot{r}]^T$. Equation (23) represents (20) rearranged with the vectors plugged in.

$I\begin{bmatrix}\dot{p}\\\dot{q}\\\dot{r}\end{bmatrix}=\begin{bmatrix}U_2\\U_3\\U_4\end{bmatrix}-\begin{bmatrix}p\\q\\r\end{bmatrix}\times I\begin{bmatrix}p\\q\\r\end{bmatrix}$

(23)

Equation (24) is the result of simplifying and isolating $\dot{\omega}$, the angular velocity vector.

$\begin{bmatrix}\dot{p}\\\dot{q}\\\dot{r}\end{bmatrix}=\begin{bmatrix}\frac{U_2}{I_x}+qr\frac{I_y-I_z}{I_x}\\\\\frac{U_3}{I_y}+rp\frac{I_z-I_x}{I_y}\\\\\frac{U_4}{I_z}+pq\frac{I_x-I_y}{I_z}\end{bmatrix}$

(24)

All equations required to model the quadcopter have been derived. Table 1 refers to the parameters used.

Table 1 Quadcopter Model Parameters

$\textbf{Parameter}$	$\textbf{Description}$	$\textbf{Value}$
$m$	$\text{Quadcopter Mass}$	$\text{0.65 kg}$
$l$	$\text{Quadcopter Arm Length}$	$\text{0.23 m}$
$b$	$\text{Thrust Factor}$	$3.13\times 10^{-5}\text{ Ns}^2$
$d$	$\text{Torque Factor}$	$7.5\times 10^{-7}\text{ Ns}^2$
$I_x$	$\text{Moment of Inertia about x-axis}$	$7.5\times 10^{-3}\text{ Ns}^2$
$I_y$	$\text{Moment of Inertia about y-axis}$	$7.5\times 10^{-3}\text{ Ns}^2$
$I_z$	$\text{Moment of Inertia about z-axis}$	$7.5\times 10^{-3}\text{ Ns}^2$
$g$	$\text{Acceleration of Gravity}$	$9.81\text{ m/s}^2$

A PD controller was utilized for altitude, roll, pitch, and yaw angle control. A PD controller for altitude is plausible because the force of gravity on the quadcopter is known (25). $K_p$ and $K_d$ denote the proportional and derivative gains of the controller. $e_z$ refers to the altitude error.

$u_1(t)=\frac{e_zK_p+mg-K_d\frac{d}{dt}z}{\cos(\varphi)\cos(\theta)}$

(25)

Equations (26), (27), and (28) refer to the PD controllers for the roll, pitch, and yaw angles respectively. $e_\varphi$, $e_\theta$, and $e_\psi$ denote the associated roll, pitch, and yaw errors.

$u_2(t)=e_\varphi K_p-K_d\frac{d}{dt}\varphi$

(26)

$u_3(t)=e_\theta K_p-K_d\frac{d}{dt}\theta$

(27)

$u_4(t)=e_\psi K_p-K_d\frac{d}{dt}\psi$

(28)

Table 2 refers to the gains used.

Table 2 PD Gains

$\textbf{Gain}$	$\textbf{Altitude}$	$\textbf{Roll}$	$\textbf{Pitch}$	$\textbf{Yaw}$
$K_p$	$\text{1.8}$	$\text{0.4}$	$\text{0.6}$	$\text{0.3}$
$K_d$	$\text{2}$	$\text{0.2}$	$\text{0.2}$	$\text{0.2}$

Notice the use of the derivative on measurement versus taking the derivative on error. This is implemented to eliminate derivative kick by assuming the setpoint is constant (29)

$\frac{d\text{Error}}{dt}=\frac{d\text{Setpoint}}{dt}-\frac{d\text{Input}}{dt}=-\frac{d\text{Input}}{dt},\space\frac{d\text{Setpoint}}{dt}=0$

(29)

Simulations

Simulation of the quadcopter plant system can be accomplished using the equations of motion (16) and (24), with (19) serving as an equation to transform the angular velocities into the time rate change of the Euler angles. Numerical integration of variables $\ddot{x}$, $\dot{x}$, $\ddot{y}$, $\dot{y}$, $\ddot{z}$, $\dot{z}$, $\dot{p}$, $\dot{q}$, $\dot{r}$, $\dot{\varphi}$, $\dot{\theta}$, and $\dot{\psi}$ is achieved using the Euler method (30). Where $y$ denotes the integrated variable and $h$ the variable to integrate. The step size, $h$, was set to 0.001 seconds.

$y(t+h)\approx y(t)+h\dot{y}(t)$

(30)

The feedback control system is simulated using the PD control equations (25), (26), (27), and (28). The derivative term is calculated by taking the backward difference (31). Where $\dot{y}$ denotes the derivative of the variable $y$. The sampling rate, $h$, is equal to the step size of 0.001 seconds.

$\dot{y}(t)=\frac{y(t)-y(t-h)}{h}$

(31)

Results

Figure 2 Simulated Altitude PD Response

Figure 3 Paper Altitude PD Response [1]

Figure 2 illustrates the altitude step response of the simulated PD controller with a settling time of 2.510 seconds and percent overshoot of 2.9683%. Figure 3 illustrates the altitude plot provided from the paper that was emulated.

Figure 4 Simulated Roll Angle PD Response

Figure 5 Paper Roll Angle PD Response [1]

Figure 4 illustrates the roll angle step response with a settling time of 1.933 seconds and negligible percent overshoot. Figure 5 illustrates the roll angle plot provided from the paper that was emulated.

Figure 6 Simulated Pitch Angle PD Response

Figure 7 Paper Pitch Angle PD Response [1]

Figure 6 illustrates the pitch angle step response with a settling time of 1.166 seconds and percent overshoot of 0.0248%. Figure 7 illustrates the pitch angle plot provided from the paper that was emulated.

Figure 8 Simulated Yaw Angle PD Response

Figure 9 Paper Yaw Angle PD Response [1]

Figure 8 illustrates the yaw angle step response with a settling time of 2.576 seconds and negligible percent overshoot. Figure 9 illustrates the yaw angle plot provided from the paper that was emulated.

The paper did not provide the settling time or overshoot of the responses, but by visual inspection, both associated plots for altitude, roll, pitch, and yaw are identical. Additionally, while the article did not specify the simulation method used, the Euler method is a suitable approach for effectively modeling the quadcopter control system, as evidenced by similar results.

The use of a PD controller for altitude served as a means of testing the quadcopter model rather than a reasonable control approach. It is unrealistic to program a controller that will only work if the mass of the system is always known. A more practical solution would be a PID controller for $U_1$. This allows $mg$ in $u_1(t)$ to be replaced by an integral term, which will gradually converge to the gravitational force over time.

Conclusion

The reproduction of the proposed quadcopter model in [1] was conducted with further analysis on the equations of motion and the implementation of a custom Euler integration scheme. The similarity of the setpoint responses for altitude, roll, pitch, and yaw angles indicate that the simulation was properly replicated. Furthermore, this report includes an analysis of both the transient and steady-state responses, which was not covered in the referenced paper. Future refinements could include the use of quaternions and the modeling of wind resistance effects. Nevertheless, the model produced in this report serves as a stepping stone for further developments, such as implementing a trajectory controller.

References

[1] Karahan, Mehmet, and Cosku Kasnakoglu. “Modeling and Simulation of Quadrotor UAV Using PID Controller.” Conference-processing. ECAI 2019 - International Conference – 11th Edition. IEEE, 2019.