MATLAB for Quadcopter Modeling — A Complete Tutorial

From MATLAB fundamentals to a validated 6-DOF simulation you can tune a real flight controller against.

Table of Contents

Part I — MATLAB Fundamentals

1. Environment and Orientation
2. Vectors, Matrices, and Indexing
3. Scripts, Functions, and Code Organization
4. Plotting and Visualization
5. Solving ODEs — the Engine of Simulation

Part II — Quadcopter Theory 6. Frames, Notation, and Conventions 7. Rotation Representations 8. Rigid-Body Dynamics 9. Rotor Aerodynamics and the Mixer 10. Motor and ESC Dynamics

Part III — Building the Simulation 11. Parameter File 12. The Dynamics Function 13. Integrating and Animating 14. Trimming and Linearization

Part IV — Control 15. Cascaded PID Architecture 16. Rate Loop 17. Attitude Loop 18. Position and Velocity Loops 19. Motor Saturation and Anti-Windup 20. LQR as an Alternative

Part V — Realism and Validation 21. Sensor Models and Noise 22. Complementary and Extended Kalman Filters 23. Discretization and Fixed-Point Effects 24. Monte Carlo and Robustness Testing 25. Exporting to Embedded Code

Appendices — Simulink track, symbolic derivations, parameter identification, common errors

Part I — MATLAB Fundamentals

1. Environment and Orientation

1.1 The Interface

MATLAB’s window is divided into panels that matter in this order:

• Command Window — the REPL. Type an expression, get an answer. Good for experiments, bad for anything you want to repeat.
• Editor — where scripts (.m files) live. All real work goes here.
• Workspace — every variable currently in memory, with size and type. Double-click any variable to open it in the Variable Editor, which is essentially a spreadsheet view. This is invaluable when a matrix is the wrong shape and you can’t tell why.
• Current Folder — MATLAB only sees files on its path. If a function “doesn’t exist,” this is usually why.

1.2 Toolboxes You Need

Check what you have:

ver

For this tutorial:

The core tutorial runs on base MATLAB alone. Everything else is optional enrichment. GNU Octave will run most of Parts I–III with minor syntax adjustments.

1.3 First Commands

clear; clc; close all;

Put this at the top of every script. clear wipes variables, clc clears the console, close all shuts figure windows. Skipping this is the single most common source of “why did my results change?” — a stale variable from a previous run silently persists.

Semicolons suppress output. Omit one and MATLAB prints the result, which is a fast debugging tool but a disaster if the variable is a 10000×12 matrix.

a = 5          % prints: a = 5
b = 5;         % silent

1.4 Getting Help

help ode45          % terse text summary
doc ode45           % full documentation with examples
which ode45         % where the file lives — catches name shadowing

which matters more than it looks. If you name a script plot.m, you have just broken plotting for your entire session.

2. Vectors, Matrices, and Indexing

MATLAB is built around matrices. Everything is a matrix; a scalar is 1×1.

2.1 Construction

r = [1 2 3];           % row vector, 1×3
c = [1; 2; 3];         % column vector, 3×1
A = [1 2; 3 4];        % 2×2 matrix

z = zeros(3,1);        % column of zeros
o = ones(2,4);
I = eye(3);            % identity

Column vectors are the convention in dynamics. State vectors, positions, velocities — all columns. Mixing rows and columns produces either a dimension error (good, you find it) or an unintended outer product via implicit expansion (bad, silent garbage).

[1 2 3] + [1; 2; 3]    % NOT an error — gives a 3×3 matrix!

This is implicit expansion, and it has silently corrupted more simulations than any other MATLAB feature. Be deliberate about orientation.

2.2 Ranges

t = 0:0.01:10;              % 0 to 10, step 0.01 → 1×1001
t = linspace(0, 10, 1001);  % same thing, specify count not step

linspace is safer when the endpoint matters — floating-point accumulation in the colon operator can drop the final element.

2.3 Indexing

MATLAB indexes from 1, not 0.

v = [10 20 30 40 50];

v(1)          % 10
v(end)        % 50
v(2:4)        % [20 30 40]
v([1 5])      % [10 50]
v(v > 25)     % [30 40 50] — logical indexing

For matrices, (row, column):

A = magic(4);
A(2,3)        % single element
A(2,:)        % entire row 2
A(:,3)        % entire column 3
A(1:2, 3:4)   % submatrix
A(:)          % flatten to column vector

The colon means “everything along this dimension.” You will use A(:,k) constantly to pull the k-th time sample from a state history matrix.

2.4 Element-wise vs Matrix Operations

This distinction is central and trips up nearly everyone.

A * B      % matrix multiplication — inner dimensions must agree
A .* B     % element-wise — dimensions must match exactly

A^2        % A*A (matrix power)
A.^2       % each element squared

A / B      % right division: A*inv(B)
A ./ B     % element-wise division

A'         % conjugate transpose
A.'        % plain transpose (use this for real matrices)

The dot means “element-wise.” In dynamics code, rotation matrix applications use *; scaling a vector of motor speeds uses .*.

2.5 Preallocation

MATLAB arrays that grow inside a loop force a full memory reallocation each iteration — O(n²) behavior.

% Bad — grows each iteration
for k = 1:10000
    result(k) = k^2;
end

% Good — preallocate
result = zeros(1, 10000);
for k = 1:10000
    result(k) = k^2;
end

For a simulation logging 12 states over 100,000 steps, this is the difference between seconds and minutes.

2.6 Vectorization

Loops in MATLAB are slow; array operations are fast because they drop into optimized libraries.

% Loop
for k = 1:length(t)
    y(k) = sin(t(k)) * exp(-0.1*t(k));
end

% Vectorized — same result, far faster
y = sin(t) .* exp(-0.1*t);

Write dynamics functions to accept vectors where practical. That said: correctness first. A clear loop that works beats a clever vectorization that doesn’t.

2.7 Structs

Structs group related data under one name. This is how you’ll manage parameters.

p.m     = 0.68;      % mass, kg
p.g     = 9.81;
p.Ixx   = 6.9e-3;
p.arm   = 0.17;

disp(p.m)
fieldnames(p)

Passing one struct into a function beats passing fourteen loose arguments. It also means adding a parameter later doesn’t require changing every function signature.

3. Scripts, Functions, and Code Organization

3.1 Scripts

A script is a .m file of commands executed top to bottom, sharing the base workspace.

% simulate_hover.m
clear; clc; close all;
params = quad_params();
x0 = zeros(12,1);
[t, x] = ode45(@(t,x) quad_dynamics(t,x,params), [0 10], x0);
plot(t, x(:,3));

3.2 Functions

Functions have their own scope — variables inside don’t leak out. This is what makes them safe.

function xdot = quad_dynamics(t, x, params)
%QUAD_DYNAMICS  6-DOF rigid body derivative for a quadcopter.
%   XDOT = QUAD_DYNAMICS(T, X, PARAMS) returns the 12×1 state derivative.
%
%   Inputs:
%     t      - time (s), unused in autonomous form but required by ode45
%     x      - 12×1 state vector [pos; vel; euler; omega]
%     params - struct from quad_params()
%
%   Output:
%     xdot   - 12×1 derivative

    % ... body ...
end

The comment block immediately after the signature is what help quad_dynamics prints. Write it. Six months from now you will not remember whether element 7 was roll or pitch.

3.3 Local and Anonymous Functions

Local functions live at the bottom of a file and are visible only within it — good for helpers you don’t want cluttering the namespace.

function main()
    y = helper(3);
end

function out = helper(x)
    out = x^2;
end

Anonymous functions capture variables at creation time:

f = @(x) x.^2 + 3;
f(2)                          % 7

k = 5;
g = @(x) k*x;
k = 100;
g(2)                          % still 10 — k was captured as 5

That capture behavior is exactly why @(t,x) quad_dynamics(t,x,params) works: params is frozen into the handle, letting you pass extra arguments through ode45, which only ever calls f(t,x).

3.4 Project Layout

quadsim/
├── params/
│   └── quad_params.m
├── dynamics/
│   ├── quad_dynamics.m
│   ├── rotation_matrix.m
│   └── mixer.m
├── control/
│   ├── pid_controller.m
│   └── attitude_controller.m
├── analysis/
│   ├── trim_quad.m
│   └── linearize_quad.m
├── plotting/
│   ├── plot_states.m
│   └── animate_quad.m
└── run_simulation.m

Add subfolders to the path at the top of your main script:

addpath(genpath(fileparts(mfilename('fullpath'))));

mfilename('fullpath') gives the location of the currently running file, so this works regardless of where MATLAB was launched from.

4. Plotting and Visualization

4.1 Basics

t = linspace(0, 10, 500);
y = exp(-0.3*t) .* sin(2*pi*t);

figure;
plot(t, y, 'LineWidth', 1.5);
grid on;
xlabel('Time (s)');
ylabel('Amplitude');
title('Damped Oscillation');

Always label axes with units. A plot without units is not evidence.

4.2 Multiple Series

figure; hold on;
plot(t, x(:,1), 'r-',  'DisplayName', 'x');
plot(t, x(:,2), 'g--', 'DisplayName', 'y');
plot(t, x(:,3), 'b-.', 'DisplayName', 'z');
hold off;
legend('show', 'Location', 'best');
grid on;

hold on prevents each new plot from wiping the axes. DisplayName plus legend('show') is more maintainable than passing a list of strings to legend — reorder your plot calls and the labels follow automatically.

4.3 Subplots

figure('Position', [100 100 900 700]);

subplot(3,2,1);
plot(t, x(:,1)); ylabel('x (m)'); grid on;

subplot(3,2,2);
plot(t, rad2deg(x(:,7))); ylabel('\phi (deg)'); grid on;
% ... etc

subplot(rows, cols, index) — index runs left-to-right, top-to-bottom. LaTeX-style markup in labels (\phi, \omega) renders automatically.

4.4 3D Trajectory

figure;
plot3(x(:,1), x(:,2), -x(:,3), 'LineWidth', 1.5);
grid on; axis equal;
xlabel('North (m)'); ylabel('East (m)'); zlabel('Up (m)');
view(45, 25);

Note the -x(:,3). In NED coordinates z points down, so negate it for an intuitive “up” axis. axis equal prevents the visual distortion that makes a gentle spiral look like a violent corkscrew.

4.5 Saving Figures

exportgraphics(gcf, 'response.png', 'Resolution', 300);

exportgraphics (R2020a+) produces tight, correctly-sized output. The older saveas adds enormous whitespace margins.

5. Solving ODEs — the Engine of Simulation

Everything in Part III rests on this.

5.1 The Standard Form

MATLAB solvers require first-order form:

$\dot{\mathbf{x}} = f(t, \mathbf{x})$

A second-order system becomes first-order by introducing velocity as a state. For a mass-spring-damper $m\ddot{y} + c\dot{y} + ky = 0$, let $x_1 = y$, $x_2 = \dot{y}$:

$\dot{x}_1 = x_2, \qquad \dot{x}_2 = -\frac{c}{m}x_2 - \frac{k}{m}x_1$

function dx = msd(t, x)
    m = 1; c = 0.5; k = 4;
    dx = [ x(2);
          -(c/m)*x(2) - (k/m)*x(1) ];
end

[t, x] = ode45(@msd, [0 20], [1; 0]);
plot(t, x(:,1)); grid on;

Your 12-state quadcopter is the same idea, just larger.

5.2 Choosing a Solver

Quadcopter dynamics with fast motor time constants (~20 ms) alongside slow position dynamics (~2 s) are mildly stiff. ode45 usually copes but may take very small steps. If your simulation crawls, try ode15s.

5.3 Options

opts = odeset('RelTol', 1e-8, ...
              'AbsTol', 1e-10, ...
              'MaxStep', 0.01, ...
              'Stats',  'on');

[t, x] = ode45(@(t,x) quad_dynamics(t,x,p), [0 10], x0, opts);

• RelTol/AbsTol — tighten these when results change with tolerance. If they do, you haven’t converged.
• MaxStep — force a ceiling on step size. Critical when your input changes discontinuously (a step command); otherwise the solver may step straight over the transition.
• Stats — reports steps taken and failed attempts. A high failure count signals stiffness.

5.4 Fixed-Step Integration

Variable-step solvers are wrong for controller design, because a real flight controller runs at a fixed rate. To match reality, integrate at fixed step with the control law running at its true frequency:

dt   = 1/400;                  % 400 Hz control loop
T    = 10;
N    = round(T/dt);
x    = zeros(12, N+1);
x(:,1) = x0;
t    = (0:N)*dt;

for k = 1:N
    u = controller(x(:,k), setpoint, p, dt);   % runs once per step
    x(:,k+1) = rk4_step(@quad_dynamics, t(k), x(:,k), u, dt, p);
end

With RK4:

function xn = rk4_step(f, t, x, u, dt, p)
    k1 = f(t,        x,           u, p);
    k2 = f(t+dt/2,   x+dt/2*k1,   u, p);
    k3 = f(t+dt/2,   x+dt/2*k2,   u, p);
    k4 = f(t+dt,     x+dt*k3,     u, p);
    xn = x + dt/6*(k1 + 2*k2 + 2*k3 + k4);
end

Use ode45 for plant validation. Use fixed-step RK4 for controller development. The variable-step solver will happily hide discretization problems that will bite you on real hardware.

5.5 Event Detection

Stop integration on a condition — ground contact, for instance:

function [value, isterminal, direction] = ground_event(t, x)
    value      = x(3);    % NED down-position; zero at ground
    isterminal = 1;       % stop integrating
    direction  = 1;       % only when increasing (descending)
end

opts = odeset('Events', @ground_event);
[t, x, te, xe, ie] = ode45(@(t,x) quad_dynamics(t,x,p), [0 20], x0, opts);

Part II — Quadcopter Theory

6. Frames, Notation, and Conventions

Conventions are where most quadcopter models go wrong. Choose, document, and never deviate.

6.1 Inertial Frame — NED

North-East-Down, fixed to the earth:

• $x_I$ → North
• $y_I$ → East
• $z_I$ → Down

Down-positive is standard in aerospace because it makes the frame right-handed with x forward, and gravity becomes simply $+g$ along $z$. The cost is that altitude is negative $z$. Get used to it, or use ENU consistently instead — but do not mix.

6.2 Body Frame — FRD

Front-Right-Down, fixed to the airframe, origin at the center of mass:

• $x_B$ → forward (nose)
• $y_B$ → right
• $z_B$ → down through the belly

Thrust acts along $-z_B$ (upward out of the top).

6.3 State Vector

Twelve states:

$\mathbf{x} = \begin{bmatrix} p_n & p_e & p_d & u & v & w & \phi & \theta & \psi & p & q & r \end{bmatrix}^T$

Velocity in the body frame is the aerodynamicist’s convention and simplifies the force equations, since thrust and drag are naturally body-fixed. Some texts use inertial velocity instead; both work, but the equations differ. This tutorial uses body-frame velocity throughout.

6.4 Rotor Numbering — X Configuration

Viewed from above, x forward:

        x_B (front)
           ↑
   M4 ◯         ◯ M1
       ＼       ／
         ＼   ／
           ✕           → y_B (right)
         ／   ＼
       ／       ＼
   M3 ◯         ◯ M2

Diagonally opposite rotors spin the same direction, so that in hover the net reaction torque cancels. Differential thrust between the CW pair and the CCW pair produces yaw.

Different autopilots number rotors differently — Betaflight, PX4, and ArduPilot all disagree. Pick one, write it down at the top of your parameter file, and verify against your actual hardware before you trust any mixer output.

7. Rotation Representations

7.1 Rotation Matrix from Euler Angles

Using the aerospace Z-Y-X (yaw-pitch-roll) sequence, the matrix rotating a body vector into the inertial frame:

$R_{IB} = R_z(\psi)R_y(\theta)R_x(\phi)$