1. 四旋翼无人机建模基础
四旋翼无人机作为一种典型的垂直起降飞行器，其建模过程需要从物理结构到动力学特性进行全面考虑。与固定翼无人机不同，四旋翼通过调节四个旋翼的转速来实现姿态和位置控制，这种独特的控制方式带来了建模上的特殊挑战。
1.1 物理结构参数定义
在MATLAB中建立四旋翼模型时，首先需要定义其物理参数。这些参数直接影响无人机的动力学特性：
matlab复制% 无人机物理参数定义
params.mass = 1.2; % 质量(kg)
params.arm_length = 0.25; % 机臂长度(m)
params.Ixx = 0.034; % X轴转动惯量(kg·m²)
params.Iyy = 0.034; % Y轴转动惯量(kg·m²)
params.Izz = 0.06; % Z轴转动惯量(kg·m²)
params.k = 1.5e-5; % 升力系数(N·s²/rad²)
params.b = 2.5e-7; % 阻力系数(N·m·s²/rad²)
params.g = 9.81; % 重力加速度(m/s²)
这些参数中，转动惯量(Ixx, Iyy, Izz)的准确测量尤为重要。实际工程中，可以通过悬摆实验或CAD软件计算得到。升力系数k和阻力系数b需要通过风洞实验或电机-螺旋桨台架测试确定。
1.2 坐标系系统建立
四旋翼建模涉及两个主要坐标系：
- 机体坐标系(Body Frame)：原点在无人机质心，X轴向前，Y轴向左，Z轴向上
- 惯性坐标系(Inertial Frame)：固定于地面的参考坐标系
两个坐标系间的转换通过欧拉角(φ,θ,ψ)描述：
- 滚转角φ：绕X轴旋转
- 俯仰角θ：绕Y轴旋转
- 偏航角ψ：绕Z轴旋转
旋转矩阵R从机体坐标系到惯性坐标系的转换关系为：
matlab复制R = [cosθ*cosψ, sinφ*sinθ*cosψ-cosφ*sinψ, cosφ*sinθ*cosψ+sinφ*sinψ;
cosθ*sinψ, sinφ*sinθ*sinψ+cosφ*cosψ, cosφ*sinθ*sinψ-sinφ*cosψ;
-sinθ, sinφ*cosθ, cosφ*cosθ];
1.3 动力学方程推导
四旋翼的六自由度动力学模型包括平移运动和旋转运动：
平移动力学：
math复制m\ddot{\mathbf{p}} = m\mathbf{g} + R\mathbf{F}_b
其中p为位置向量，F_b为机体坐标系下的总升力。
旋转动力学：
math复制I\dot{\boldsymbol{\omega}} + \boldsymbol{\omega}\times I\boldsymbol{\omega} = \boldsymbol{\tau}
其中ω为角速度向量，τ为总力矩。
在MATLAB中，这些方程通常转化为状态空间形式进行求解。典型的状态向量包含12个元素：
matlab复制x = [p_x p_y p_z v_x v_y v_z φ θ ψ ω_x ω_y ω_z]'
实际建模时需要注意：上述模型假设机体为刚体且螺旋桨升力瞬时响应，忽略了空气动力学效应和电机动态特性。对于高精度仿真，还需要考虑这些因素的影响。
2. 独立旋转机翼的实现
传统四旋翼的四个旋翼通常两两反向旋转以抵消反扭矩。而实现机翼独立旋转需要修改这一约束条件，为控制系统设计带来新的可能性和挑战。
2.1 旋翼配置方案
独立旋转机翼系统允许每个旋翼的转速和旋转方向完全独立控制。在MATLAB中，我们可以定义旋翼的旋转方向矩阵：
matlab复制% 传统四旋翼旋转方向 (1:顺时针，-1:逆时针)
rotor_direction = [1 -1 1 -1];
% 独立旋转机翼方向 (可动态变化)
independent_direction = [1 1 -1 -1]; % 示例配置
这种灵活性带来了几个优势：
- 故障容错能力：单个旋翼故障时可通过调整其他旋翼补偿
- 机动性提升：可实现更复杂的力矩组合
- 能量效率优化：根据飞行状态调整旋翼方向
2.2 控制分配矩阵
独立旋转机翼系统的控制分配矩阵与传统四旋翼不同。控制效率矩阵E将旋翼转速(ω₁²,ω₂²,ω₃²,ω₄²)映射到总升力F和力矩τ：
matlab复制% 传统X型布局控制分配矩阵
E_classic = [k k k k;
0 -k*L 0 k*L;
-k*L 0 k*L 0;
b -b b -b];
% 独立旋转机翼分配矩阵 (L为机臂长度)
E_independent = [k k k k;
-k*sin(pi/4) k*sin(pi/4) k*sin(pi/4) -k*sin(pi/4);
-k*cos(pi/4) -k*cos(pi/4) k*cos(pi/4) k*cos(pi/4);
b b -b -b];
在MATLAB中，可以通过伪逆计算所需的旋翼转速：
matlab复制omega_squared = pinv(E_independent) * [F; tau_x; tau_y; tau_z];
2.3 动态方向调整算法
实现机翼独立旋转的核心是动态调整算法。以下是一个简单的MATLAB实现示例：
matlab复制function [direction, omega] = adjust_rotor_direction(...
desired_F, desired_tau, current_attitude)
% 根据当前状态和期望力/力矩计算最优旋翼方向和转速
% 初始化参数
k = 1.5e-5; b = 2.5e-7; L = 0.25;
options = optimoptions('fmincon', 'Display', 'off');
% 优化目标：最小化转速平方和（节省能量）
fun = @(x) sum(x(1:4).^2);
% 决策变量：[ω1², ω2², ω3², ω4², d1, d2, d3, d4]
% di ∈ [-1,1]表示旋转方向
% 约束条件：满足力/力矩需求
Aeq = [k k k k;
-k*L k*L k*L -k*L;
-k*L -k*L k*L k*L;
b b -b -b];
beq = [desired_F; desired_tau];
% 优化求解
x_opt = fmincon(fun, zeros(8,1), [], [], Aeq, beq, ...
[0;0;0;0;-1;-1;-1;-1], ...
[inf;inf;inf;inf;1;1;1;1], [], options);
omega = sqrt(x_opt(1:4));
direction = sign(x_opt(5:8));
end
实际应用中，还需要考虑电机响应延迟、转速变化率限制等实际问题。建议先用简化模型验证算法可行性，再逐步增加复杂度。
3. MATLAB仿真实现
有了数学模型后，我们需要在MATLAB中实现完整的仿真系统。这一部分将介绍如何使用Simulink搭建四旋翼无人机仿真模型。
3.1 Simulink模型架构
一个完整的四旋翼仿真模型通常包含以下子系统：
- 环境子系统：实现重力、空气阻力等环境效应
- 控制器子系统：实现飞行控制算法
- 无人机动力学子系统：实现六自由度运动方程
- 传感器模型子系统：模拟IMU、GPS等传感器数据
- 可视化子系统：提供三维动画显示
在Simulink中创建模型的基本框架：
matlab复制% 新建Simulink模型
model = 'quadcopter_simulation';
new_system(model);
open_system(model);
% 添加主要子系统
add_block('simulink/Commonly Used Blocks/Subsystem', [model '/Environment']);
add_block('simulink/Commonly Used Blocks/Subsystem', [model '/Controller']);
add_block('simulink/Commonly Used Blocks/Subsystem', [model '/Dynamics']);
add_block('simulink/Commonly Used Blocks/Subsystem', [model '/Sensors']);
add_block('simulink/Commonly Used Blocks/Subsystem', [model '/Visualization']);
% 连接子系统
add_line(model, 'Controller/1', 'Dynamics/1');
add_line(model, 'Dynamics/1', 'Sensors/1');
add_line(model, 'Sensors/1', 'Controller/2');
add_line(model, 'Dynamics/2', 'Visualization/1');
3.2 动力学子系统实现
动力学子系统是模型的核心，实现六自由度运动方程。以下是关键部分的MATLAB Function Block实现：
matlab复制function dx = dynamics(x, F, tau, params)
% 状态向量: [x y z vx vy vz phi theta psi omega_x omega_y omega_z]
% 输入: 总升力F, 力矩tau=[tau_x;tau_y;tau_z], 参数params
% 位置导数
v = x(4:6);
p_dot = v;
% 速度导数 (惯性系)
R = rotation_matrix(x(7:9));
a = [0; 0; -params.g] + R*[0; 0; F]/params.mass;
% 欧拉角导数
phi = x(7); theta = x(8);
H = [1, sin(phi)*tan(theta), cos(phi)*tan(theta);
0, cos(phi), -sin(phi);
0, sin(phi)/cos(theta), cos(phi)/cos(theta)];
euler_dot = H * x(10:12);
% 角速度导数 (机体系)
omega = x(10:12);
I = [params.Ixx, 0, 0;
0, params.Iyy, 0;
0, 0, params.Izz];
omega_dot = I \ (tau - cross(omega, I*omega));
% 组合状态导数
dx = [p_dot; a; euler_dot; omega_dot];
end
function R = rotation_matrix(euler)
phi = euler(1); theta = euler(2); psi = euler(3);
R = [cos(theta)*cos(psi), sin(phi)*sin(theta)*cos(psi)-cos(phi)*sin(psi), cos(phi)*sin(theta)*cos(psi)+sin(phi)*sin(psi);
cos(theta)*sin(psi), sin(phi)*sin(theta)*sin(psi)+cos(phi)*cos(psi), cos(phi)*sin(theta)*sin(psi)-sin(phi)*cos(psi);
-sin(theta), sin(phi)*cos(theta), cos(phi)*cos(theta)];
end
3.3 独立旋转控制实现
在控制器子系统中，我们需要实现独立旋转机翼的控制逻辑。典型的级联控制结构如下：
- 外环位置控制：PID控制生成期望姿态
- 内环姿态控制：生成所需力矩
- 控制分配：将力/力矩分配到四个独立旋翼
matlab复制function [omega, direction] = independent_rotor_control(...
pos_des, pos_actual, att_des, att_actual, params)
% 位置控制器 (外环)
Kp_pos = [1; 1; 1]; Ki_pos = [0.1; 0.1; 0.1]; Kd_pos = [0.5; 0.5; 0.5];
persistent pos_error_int;
if isempty(pos_error_int)
pos_error_int = zeros(3,1);
end
pos_error = pos_des - pos_actual;
pos_error_int = pos_error_int + pos_error;
F_des = params.mass * params.g + ...
Kp_pos.*pos_error + Ki_pos.*pos_error_int + Kd_pos.*(pos_error - [0;0;0]);
% 姿态控制器 (内环)
Kp_att = [10; 10; 5]; Kd_att = [2; 2; 1];
att_error = att_des - att_actual;
tau_des = Kp_att.*att_error + Kd_att.*(att_error - [0;0;0]);
% 控制分配
[direction, omega] = adjust_rotor_direction(...
F_des(3), tau_des, att_actual);
end
3.4 可视化实现
使用MATLAB的Aerospace工具箱或Simulink 3D Animation工具箱可以实现三维可视化：
matlab复制function visualize_quadcopter(pos, att, rotor_speeds)
persistent fig quad_handle rotor_handles;
if isempty(fig) || ~isvalid(fig)
fig = figure('Name', 'Quadcopter Animation');
ax = axes('Parent', fig);
axis equal; grid on; hold on;
xlabel('X (m)'); ylabel('Y (m)'); zlabel('Z (m)');
view(3);
% 绘制无人机主体
quad_handle = plot3(0,0,0,'ro','MarkerSize',10,'MarkerFaceColor','r');
% 绘制四个旋翼
rotor_pos = [0.25 0 0; 0 0.25 0; -0.25 0 0; 0 -0.25 0];
rotor_handles = gobjects(4,1);
for i = 1:4
rotor_handles(i) = plot3(0,0,0,'bo','MarkerSize',8+rotor_speeds(i)*5,...
'MarkerFaceColor','b');
end
end
% 更新位置和姿态
R = rotation_matrix(att);
set(quad_handle, 'XData', pos(1), 'YData', pos(2), 'ZData', pos(3));
% 更新旋翼位置和大小
rotor_pos = [0.25 0 0; 0 0.25 0; -0.25 0 0; 0 -0.25 0];
for i = 1:4
rotor_pos_rotated = pos' + (R*rotor_pos(i,:)')';
set(rotor_handles(i), 'XData', rotor_pos_rotated(1),...
'YData', rotor_pos_rotated(2),...
'ZData', rotor_pos_rotated(3),...
'MarkerSize', 8+rotor_speeds(i)*5);
end
drawnow;
end
可视化时注意调整更新频率，避免影响仿真速度。对于复杂场景，可以考虑使用Simulink 3D Animation工具箱的VR Sink块连接外部VRML文件。
4. 仿真分析与案例研究
完成模型搭建后，我们需要进行系统仿真并分析独立旋转机翼对飞行性能的影响。
4.1 基本飞行测试
首先测试基本的悬停和轨迹跟踪性能：
matlab复制% 仿真参数设置
sim_time = 10; % 仿真时间(s)
dt = 0.01; % 步长(s)
steps = sim_time/dt;
% 初始状态
x = zeros(12,1); % [位置;速度;欧拉角;角速度]
x(3) = 1; % 初始高度1m
% 期望轨迹 (圆形)
t = 0:dt:sim_time;
r = 2; omega = 0.5;
pos_des = [r*sin(omega*t); r*cos(omega*t); ones(size(t))];
att_des = zeros(3,length(t));
% 仿真循环
for i = 1:steps
% 获取当前状态
pos_actual = x(1:3);
att_actual = x(7:9);
% 控制器计算
[omega, direction] = independent_rotor_control(...
pos_des(:,i), pos_actual, att_des(:,i), att_actual, params);
% 动力学更新
F = params.k * sum(omega.^2);
tau = [params.k*(omega(2)^2-omega(4)^2)*params.arm_length;
params.k*(omega(3)^2-omega(1)^2)*params.arm_length;
params.b*(omega(1)^2-omega(2)^2+omega(3)^2-omega(4)^2)];
% 数值积分 (四阶龙格-库塔)
k1 = dynamics(x, F, tau, params);
k2 = dynamics(x+0.5*dt*k1, F, tau, params);
k3 = dynamics(x+0.5*dt*k2, F, tau, params);
k4 = dynamics(x+dt*k3, F, tau, params);
x = x + dt*(k1 + 2*k2 + 2*k3 + k4)/6;
% 可视化
visualize_quadcopter(x(1:3), x(7:9), omega);
end
4.2 独立旋转优势分析
通过与传统固定旋转方向的对比，可以分析独立旋转机翼的优势：
- 能量效率：在非对称机动中，独立旋转可以优化能量分配
matlab复制% 能量消耗比较
energy_fixed = sum(omega_fixed.^2)*dt;
energy_independent = sum(omega_independent.^2)*dt;
fprintf('固定方向能耗: %.2f J\n', energy_fixed);
fprintf('独立旋转能耗: %.2f J\n', energy_independent);
fprintf('节能比例: %.1f%%\n', (energy_fixed-energy_independent)/energy_fixed*100);
- 故障容错：模拟单个旋翼失效时的恢复能力
matlab复制% 模拟第2个旋翼失效
omega(2) = 0;
% 传统系统无法完全补偿
tau_fixed = [params.k*(omega(2)^2-omega(4)^2)*params.arm_length;
params.k*(omega(3)^2-omega(1)^2)*params.arm_length;
params.b*(omega(1)^2-omega(2)^2+omega(3)^2-omega(4)^2)];
% 独立旋转系统可通过调整方向补偿
direction = [1 -1 1 -1]; % 调整旋转方向
omega_independent = omega;
omega_independent(4) = sqrt((tau_des(3)-params.b*(omega(1)^2-omega(2)^2+omega(3)^2))/params.b);
- 机动性测试：比较完成相同机动动作的时间
matlab复制% 快速滚转机动测试
att_des = [pi/4*sin(2*pi*1*t); zeros(size(t)); zeros(size(t))];
% 测量达到期望姿态的时间
settling_time_fixed = find(abs(att_actual_fixed(1,:)-att_des(1,:))<0.05,1)*dt;
settling_time_independent = find(abs(att_actual_independent(1,:)-att_des(1,:))<0.05,1)*dt;
4.3 参数敏感性分析
独立旋转系统的性能受多个参数影响，需要进行敏感性分析：
- 控制分配矩阵精度：
matlab复制% 添加不同级别的噪声测试鲁棒性
noise_level = 0:0.05:0.3;
tracking_errors = zeros(size(noise_level));
for i = 1:length(noise_level)
E_noisy = E_independent + noise_level(i)*randn(size(E_independent));
% 运行仿真并记录跟踪误差
tracking_errors(i) = mean(sqrt(sum((pos_actual-pos_des).^2)));
end
- 电机响应延迟影响：
matlab复制% 一阶延迟模型测试
tau_motor = 0.05; % 电机时间常数
omega_actual = zeros(4,1);
omega_dot_prev = zeros(4,1);
for i = 1:steps
% 电机动态模型
omega_dot = (omega_des - omega_actual)/tau_motor;
omega_actual = omega_actual + (omega_dot + omega_dot_prev)/2*dt;
omega_dot_prev = omega_dot;
% 其余仿真步骤...
end
- 转动惯量不确定性：
matlab复制% 测试±20%转动惯量变化的影响
I_variation = linspace(0.8, 1.2, 5);
performance_metrics = zeros(length(I_variation),3); % 跟踪误差,能耗,稳定时间
for i = 1:length(I_variation)
params.Ixx = 0.034 * I_variation(i);
params.Iyy = 0.034 * I_variation(i);
params.Izz = 0.06 * I_variation(i);
% 运行仿真并记录性能指标
% ...
end
参数分析时建议使用MATLAB的Design of Experiments (DOE)工具箱或Simulink Design Optimization工具箱进行系统化测试。对于复杂参数空间，可以考虑使用机器学习方法建立代理模型。