1. 四旋翼无人机建模概述
四旋翼无人机作为一种典型的垂直起降飞行器,其建模与仿真一直是控制领域的热点课题。与传统固定翼飞机不同,四旋翼通过四个独立控制的旋翼实现飞行控制,这种独特的结构使其具有悬停、垂直起降等特殊能力。在MATLAB环境下建立四旋翼模型,不仅可以直观展示其飞行原理,还能为后续控制算法开发提供验证平台。
提示:完整的四旋翼建模需要考虑动力学方程、空气动力学效应、电机响应特性等多个方面,本文将从基础建模开始,逐步构建完整的仿真系统。
1.1 四旋翼基本结构
典型的四旋翼采用十字形或X形布局,四个旋翼对称分布在机体四个端点。根据旋翼旋转方向的不同,可分为"+"型和"X"型两种配置:
- "+"型配置:前后旋翼同向旋转,左右旋翼同向旋转
- "X"型配置:对角线上的旋翼同向旋转
本文采用更常见的X型配置,其优势在于飞行控制时能产生更均匀的力矩分布。每个旋翼由独立电机驱动,通过改变电机转速来控制升力大小。
1.2 坐标系定义
建立数学模型前,需要明确定义坐标系:
- 地面惯性坐标系(E系):固定于地面,遵循"北-东-地"右手定则
- 机体坐标系(B系):固连于无人机,原点在重心
- X轴:指向机头方向
- Y轴:指向右侧
- Z轴:垂直向下
两坐标系间的转换通过欧拉角(φ,θ,ψ)描述,分别对应滚转、俯仰和偏航角。
2. 动力学建模
2.1 运动方程推导
四旋翼的六自由度运动方程可分为平移运动和旋转运动两部分:
平移运动方程:
m(d²x/dt²) = (cosφsinθcosψ + sinφsinψ)U₁
m(d²y/dt²) = (cosφsinθsinψ - sinφcosψ)U₁
m(d²z/dt²) = -mg + (cosφcosθ)U₁
旋转运动方程:
Iₓₓ(d²φ/dt²) = (Iᵧᵧ - Izz)θ̇ψ̇ + l(U₂ - U₄)
Iᵧᵧ(d²θ/dt²) = (Izz - Iₓₓ)φ̇ψ̇ + l(U₃ - U₁)
Izz(d²ψ/dt²) = (Iₓₓ - Iᵧᵧ)φ̇θ̇ + (U₁ - U₂ + U₃ - U₄)
其中:
- U₁ = b(ω₁² + ω₂² + ω₃² + ω₄²) 总升力
- U₂ = b(ω₄² - ω₂²) 滚转力矩
- U₃ = b(ω₃² - ω₁²) 俯仰力矩
- U₄ = d(ω₂² + ω₄² - ω₁² - ω₃²) 偏航力矩
- b: 升力系数
- d: 阻力系数
- l: 旋翼到重心的距离
2.2 电机动力学模型
实际电机响应存在惯性,不能瞬时达到指令转速。常用一阶惯性环节描述:
G(s) = 1/(τs + 1)
其中τ为电机时间常数,典型值在0.05-0.2秒之间。在MATLAB中可通过Transfer Function模块实现。
3. MATLAB实现
3.1 参数初始化
首先定义无人机的基本参数:
matlab复制% 物理参数
m = 1.2; % 质量(kg)
g = 9.81; % 重力加速度(m/s^2)
l = 0.25; % 旋翼到重心距离(m)
% 惯性参数
Ixx = 0.023; % X轴转动惯量(kg·m^2)
Iyy = 0.023; % Y轴转动惯量(kg·m^2)
Izz = 0.046; % Z轴转动惯量(kg·m^2)
% 空气动力参数
b = 5.5e-6; % 升力系数(N·s^2)
d = 1.1e-7; % 阻力系数(N·m·s^2)
3.2 控制分配矩阵
将控制输入转换为各电机转速:
matlab复制% 控制分配矩阵
A = [b b b b;
0 -b 0 b;
-b 0 b 0;
d -d d -d];
% 计算电机转速平方
omega_sq = pinv(A) * [U1; U2; U3; U4];
3.3 可视化实现
使用MATLAB的3D绘图功能实现无人机可视化:
matlab复制function draw_quadrotor(x,y,z,phi,theta,psi)
% 机体尺寸
arm_length = 0.25;
prop_size = 0.1;
% 旋转矩阵
R = rotation_matrix(phi,theta,psi);
% 绘制机身
body = [-arm_length 0 0; arm_length 0 0;
0 -arm_length 0; 0 arm_length 0];
body = (R*body')' + repmat([x y z],4,1);
% 绘制旋翼
props = [arm_length 0 0; -arm_length 0 0;
0 arm_length 0; 0 -arm_length 0];
props = (R*props')' + repmat([x y z],4,1);
% 绘图
plot3(body([1 2],1),body([1 2],2),body([1 2],3),'k-','LineWidth',2);
hold on;
plot3(body([3 4],1),body([3 4],2),body([3 4],3),'k-','LineWidth',2);
for i=1:4
draw_propeller(props(i,:),prop_size,R);
end
axis equal; grid on;
xlabel('X'); ylabel('Y'); zlabel('Z');
view(30,30);
end
4. 独立旋翼控制实现
4.1 旋翼旋转动画
实现旋翼独立旋转的关键是实时更新每个旋翼的方位角:
matlab复制function update_propellers(omega, dt)
% omega: 各旋翼转速(rad/s)
% dt: 时间步长(s)
persistent prop_angle;
if isempty(prop_angle)
prop_angle = zeros(4,1);
end
% 更新旋翼角度
prop_angle = prop_angle + omega*dt;
% 绘制旋翼
for i=1:4
% 计算旋翼位置和方向
pos = get_prop_position(i);
dir = get_prop_direction(i);
% 绘制旋转中的旋翼
draw_rotating_propeller(pos, dir, prop_angle(i));
end
end
4.2 控制策略
基本PID控制算法实现:
matlab复制function [U1, U2, U3, U4] = pid_controller(state, desired, Kp, Ki, Kd)
% state: 当前状态 [x,y,z,φ,θ,ψ,...]
% desired: 期望状态
% Kp,Ki,Kd: PID参数
persistent integral error_prev;
if isempty(integral)
integral = zeros(6,1);
error_prev = zeros(6,1);
end
% 计算误差
error = desired(1:6) - state(1:6);
% 积分项
integral = integral + error;
% 微分项
derivative = (error - error_prev)/dt;
% PID输出
output = Kp.*error + Ki.*integral + Kd.*derivative;
% 更新误差
error_prev = error;
% 分配控制量
U1 = m*(g + output(3))/cos(state(4))/cos(state(5));
U2 = Ixx*output(4);
U3 = Iyy*output(5);
U4 = Izz*output(6);
end
5. 仿真与结果分析
5.1 仿真参数设置
matlab复制% 仿真时间
t_start = 0;
t_end = 10;
dt = 0.01;
t = t_start:dt:t_end;
% 初始状态
x0 = [0; 0; 0; 0; 0; 0; zeros(6,1)];
% 期望状态
x_des = [1; 1; 1; 0; 0; pi/4; zeros(6,1)];
% PID参数
Kp = diag([5, 5, 10, 8, 8, 4]);
Ki = diag([0.1, 0.1, 0.5, 0.3, 0.3, 0.1]);
Kd = diag([3, 3, 5, 2, 2, 1]);
5.2 仿真主循环
matlab复制% 初始化
x = x0;
results = zeros(length(t), length(x0));
for i = 1:length(t)
% 获取当前状态
state = x(:,end);
% 计算控制量
[U1, U2, U3, U4] = pid_controller(state, x_des, Kp, Ki, Kd);
% 动力学方程
dx = quadrotor_dynamics(state, U1, U2, U3, U4);
% 数值积分
x_next = state + dx*dt;
% 存储结果
results(i,:) = x_next';
x = [x x_next];
% 可视化
if mod(i,10) == 0
draw_quadrotor(x_next(1),x_next(2),x_next(3),...
x_next(4),x_next(5),x_next(6));
drawnow;
end
end
5.3 结果分析
仿真完成后,可绘制各状态变量的时间历程曲线:
matlab复制% 位置曲线
figure;
subplot(3,1,1);
plot(t, results(:,1), 'b', t, x_des(1)*ones(size(t)), 'r--');
title('X Position'); legend('Actual','Desired');
subplot(3,1,2);
plot(t, results(:,2), 'b', t, x_des(2)*ones(size(t)), 'r--');
title('Y Position');
subplot(3,1,3);
plot(t, results(:,3), 'b', t, x_des(3)*ones(size(t)), 'r--');
title('Z Position');
% 姿态曲线
figure;
subplot(3,1,1);
plot(t, rad2deg(results(:,4)), 'b', t, rad2deg(x_des(4))*ones(size(t)), 'r--');
title('Roll Angle');
subplot(3,1,2);
plot(t, rad2deg(results(:,5)), 'b', t, rad2deg(x_des(5))*ones(size(t)), 'r--');
title('Pitch Angle');
subplot(3,1,3);
plot(t, rad2deg(results(:,6)), 'b', t, rad2deg(x_des(6))*ones(size(t)), 'r--');
title('Yaw Angle');
6. 常见问题与调试技巧
6.1 仿真发散问题
当仿真结果出现发散时,可检查以下几个方面:
- 时间步长设置:过大的dt会导致数值积分不稳定,建议从0.001s开始尝试
- PID参数调节:先调节P项使系统稳定,再加入I和D项
- 控制量限幅:实际电机有最大转速限制,需在控制算法中加入饱和限制
6.2 姿态控制振荡
姿态控制出现振荡通常是因为:
- 微分增益过大:减小Kd参数
- 传感器噪声:在仿真中加入低通滤波器
- 电机延迟:增加电机动力学模型的时间常数
6.3 可视化性能优化
当动画卡顿时,可以:
- 减少更新频率,如每10步更新一次图形
- 使用轻量级的绘图函数,如line代替plot3
- 关闭不必要的图形属性,如阴影、光照等
注意:在调试过程中,建议先验证各子模块的正确性,如单独测试动力学模型、控制算法等,再集成到完整系统中。