✅作者简介：热爱科研的Matlab仿真开发者，擅长毕业设计辅导、数学建模、数据处理、程序设计科研仿真。

🍎完整代码获取 定制创新 论文复现点击：Matlab科研工作室

👇 关注我领取海量matlab电子书和数学建模资料

🍊个人信条：做科研，博学之、审问之、慎思之、明辨之、笃行之，是为：博学慎思，明辨笃行。

🔥 内容介绍

一、引言

四旋翼飞行器因其结构简单、机动性强，在航拍、物流、监测等领域得到广泛应用。深入理解其六自由度动力学特性，并设计有效的 PID 控制系统，对于实现四旋翼飞行器的稳定、精准飞行至关重要。通过动力学仿真，我们能够模拟飞行器在各种情况下的运动状态，为 PID 控制器的设计与优化提供依据。

二、六自由度四旋翼动力学模型

（一）基本假设与坐标系定义

1. 基本假设

   ：假设四旋翼飞行器为刚体，忽略其弹性变形；旋翼产生的升力和扭矩与旋翼转速呈线性关系；不考虑空气压缩性等复杂空气动力学效应。

⛳️ 运行结果

📣 部分代码

function dydt = Motion_ABC(t, y, Omega_hover, I, m, g, l, k, b)

% Get phi, theta, psi, p, q, and r

phi = y(7); theta = y(8); psi = y(9);

p = y(10); q = y(11); r = y(12);

% Create omega for rotors

Omega = Omega_hover * ones(4,1);

% If statements to find Omega

if t < 1

Omega = Omega_hover + 70 * sin(2*pi*t/4);

Omega = ones(4,1) * Omega;

elseif t < 2

Omega = Omega_hover - 77 * sin(2*pi*t/4);

Omega = ones(4,1) * Omega;

elseif t < 3

Omega(2) = sqrt(Omega_hover^2 - 70^2 * sin(2*pi*(t-2)/4));

Omega(4) = sqrt(Omega_hover^2 + 70^2 * sin(2*pi*(t-2)/4));

elseif t < 4

Omega(2) = sqrt(Omega_hover^2 + 70^2 * sin(2*pi*(t-2)/4));

Omega(4) = sqrt(Omega_hover^2 - 70^2 * sin(2*pi*(t-2)/4));

elseif t < 5

Omega(1) = sqrt(Omega_hover^2 - 70^2 * sin(2*pi*(t-4)/4));

Omega(3) = sqrt(Omega_hover^2 + 70^2 * sin(2*pi*(t-4)/4));

else

Omega(1) = sqrt(Omega_hover^2 + 70^2 * sin(2*pi*(t-4)/4));

Omega(3) = sqrt(Omega_hover^2 - 70^2 * sin(2*pi*(t-4)/4));

end

% Get T and Tau

T = k * sum(Omega.^2);

L = k * l * (-Omega(2)^2 + Omega(4)^2);

M = k * l * (-Omega(1)^2 + Omega(3)^2);

N = b * (Omega(1)^2 - Omega(2)^2 + Omega(3)^2 - Omega(4)^2);

Tau = [L; M; N];

% Rotation matrix (inertial to body)

C_IB = [ cos(theta)*cos(psi), cos(theta)*sin(psi), -sin(theta);

sin(phi)*sin(theta)*cos(psi)-cos(phi)*sin(psi), sin(phi)*sin(theta)*sin(psi)+cos(phi)*cos(psi), sin(phi)*cos(theta);

cos(phi)*sin(theta)*cos(psi)+sin(phi)*sin(psi), cos(phi)*sin(theta)*sin(psi)-sin(phi)*cos(psi), cos(phi)*cos(theta)];

% Body to inertial rotation

C_BI = C_IB';

Thrust = C_BI(:,3);

% Get r double dot

r_dd = [0;0;-g] + (T/m) * Thrust;

% Get w and Omega_dot

w = [p; q; r];

Omega_dot = I \ ( - cross(w, I*w) + Tau);

% Euler angle rates

A = [1, sin(phi)*tan(theta), cos(phi)*tan(theta);

0, cos(phi), -sin(phi);

0, sin(phi)/cos(theta), cos(phi)/cos(theta)];

Euler = A \ w;

if t > 4

Euler(1) = 0;

end

% dydt

dydt = zeros(12,1);

dydt(1:3) = y(4:6);

dydt(4:6) = r_dd;

dydt(7:9) = Euler;

dydt(10:12) = Omega_dot;

end

🔗 参考文献

🍅更多免费数学建模和仿真教程关注领取