十年匠心定制 · 商业建站与技术教学双线并行 咨询热线:400-886-1026 service@lmnt.cn
ARTICLE DETAIL

资讯详情

深耕网站建设与运营推广的一线实战洞察。

四旋翼飞行器MATLAB仿真:从动力学模型到Simulink控制调参全流程

四旋翼飞行器MATLAB仿真:从动力学模型到Simulink控制调参全流程 简介四旋翼飞行器MATLAB仿真程序聚焦Simulink环境下的无人机动力学建模与控制系统设计适合无人机、自动化及控制工程方向的学习者和研究人员。压缩包共17个文件含15个MATLAB脚本和2个Simulink模型mdl整体大小约53KB脚本覆盖姿态控制、速度控制、气动模型、参数初始化、结果绘图等功能两个模型则提供系统级仿真框架。已有264人学习下载。通过该资源读者可掌握四旋翼六自由度动力学建模、旋翼升力与扭矩计算、PID或状态反馈控制器设计以及传感器反馈与轨迹规划的实现思路同时还能了解如何利用MATLAB代码生成工具将模型部署到硬件进行在环测试。对于希望快速搭建无人机仿真环境、验证控制算法的人员来说这是一份结构紧凑、可直接运行的参考资料。1. 四旋翼飞行器MATLAB仿真程序里最容易被忽略的一件事拿到「四旋翼飞行器MATLAB仿真程序_包括simulink仿真系统和程序.zip」这类压缩包多数人的第一反应是解压、双击 .slx、点绿色运行键。但真正做过飞控仿真的人都知道第一跑大概率是红的要么报“未定义变量”要么波形飞上天要么干脆卡死。问题往往不在模型本身而在初始化脚本、求解器配置和坐标系约定这三件事上。这篇博文就按我平时搭四旋翼仿真环境的顺序展开——先把动力学和控制律写成 MATLAB 程序再装进 Simulink 仿真系统最后处理参数整定、批量仿真和半实物衔接。适合刚接触无人机仿真的学生也适合想把 MATLAB 里那套仿真程序搬到 PX4 或自研飞控里的工程师。2. 仿真程序先立模型四旋翼动力学与PID控制律的常见写法2.1 坐标系、状态量与单位约定仿真跑不跑得起来的隐藏前提四旋翼仿真程序里 90% 的“灵异现象”都出在单位上。我把自己的做法放在前面位置用米速度用米/秒角度用弧度角速度用弧度/秒推力用牛力矩用牛·米。别小看这一步很多 zip 包里动力学是厘米克秒制控制器却按国际单位写仿真出来的高度误差直接差两个数量级。坐标系我建议采用 NED北东地或 ENU东北天选一个并固定下来。常见飞控固件用 NEDMATLAB 官方 Aerospace 工具箱也偏 NED。但如果只是写脚本验证控制率ENU 更符合直觉位置 z 轴向上悬停时推力等于重力且方向为正。我一般用 NED原因只有一个后续衔接 Pixhawk 日志和真实飞控数据时NED 不需要额外转换。状态量取 12 维位置 ( p [x,y,z]^T )、欧拉角 ( \Theta [\phi,\theta,\psi]^T )、速度 ( v [v_x,v_y,v_z]^T )、机体角速度 ( \omega [p,q,r]^T )。欧拉角方式直观但要注意 90° 俯仰时出现万向锁gimbal lock。仿真程序里做悬停或小角度机动够用做大机动或特技飞行要把姿态改为四元数。2.2 控制律结构双环 PID 而不是直接调四个电机四旋翼的控制律常见套路是串级 PID外环是位置环输出期望姿态角和期望总推力内环是姿态环分角度环和角速度环最终输出四个电机的拉力指令。理由很朴素——四旋翼是欠驱动系统只有四个输入但自由度有六个你只能先控姿态再靠姿态变化产生水平加速度。直接对 x、y 位置输出四个电机转速是控不住的。我一般在 MATLAB 程序里写一个controller_attitude()和一个controller_position()。位置环用 P 或 PD姿态内环用 PID。下面是内环角速度环的离散实现这段代码可以直接放进 MATLAB Function 或普通 m 函数里function u att_rate_pid(omega_des, omega, dt, kp, ki, kd, i_limit) % 角速度环 PIDomega_des 为期望角速度omega 为当前角速度 % 输出 u 为期望力矩指令单位 N*m persistent integral persistent prev_err if isempty(integral) integral 0; prev_err 0; end err omega_des - omega; integral integral err * dt; % 积分限幅防止 windup integral max(min(integral, i_limit), -i_limit); derivative (err - prev_err) / dt; u kp * err ki * integral kd * derivative; prev_err err; end这段代码的关键在于persistent变量跨仿真步长保存积分累积量和上一次误差。i_limit是积分限幅没有它仿真一开始的误差会在前几秒把积分项推到饱和后面姿态响应会出现明显的超调和振荡。2.3 刚体动力学 m 函数仿真系统的核心程序把动力学写成独立函数Simulink 里通过 MATLAB Function 或 Level-2 S-Function 去调比把公式铺成一大片 Gain 和 Sum 模块好维护得多。下面是我常用的quad_dyn.m输入状态x(12×1)、电机拉力向量u(4×1)、参数结构体params输出状态导数function xd quad_dyn(x, u, params) % 状态分解 pos x(1:3); eul x(4:6); % [phi; theta; psi] vel x(7:9); pqr x(10:12); % [p; q; r] phi eul(1); th eul(2); psi eul(3); T sum(u); % 总推力 I params.I; % 惯量矩阵 3x3 mass params.m; g params.g; % 推力产生的线加速度NED 下注意符号 acc [ -(cos(phi)*sin(th)*cos(psi) sin(phi)*sin(psi)) / mass * T; -(cos(phi)*sin(th)*sin(psi) - sin(phi)*cos(psi)) / mass * T; g - (cos(phi)*cos(th)) / mass * T ]; % 机体轴力矩前后电机差产生俯仰左右差产生横滚对角差产生偏航 % 这里假设 u1 前、u2 右、u3 后、u4 左 M [ params.l * (u(4) - u(2)); params.l * (u(1) - u(3)); params.km * (u(1) - u(2) u(3) - u(4)) ]; % 角加速度I * omega_dot omega × (I * omega) M omega_dot I \ (M - cross(pqr, I * pqr)); % 欧拉角速率与机体角速度的转换 W [1 sin(phi)*tan(th) cos(phi)*tan(th); 0 cos(phi) -sin(phi); 0 sin(phi)/cos(th) cos(phi)/cos(th)]; eul_dot W * pqr; xd [vel; eul_dot; acc; omega_dot]; end注意上面M的计算是基于“前、右、后、左”电机编号不同 zip 包里电机编号不一定一致拿到程序先看混控器定义再看仿真结果的姿态响应方向不然你给的俯仰指令可能让飞机向后翻。代码里params是结构体I 是多旋翼惯量矩阵params.km是力矩系数单位 N·m/N。这些参数在后续初始化脚本里统一赋值不要写在函数内部。3. Simulink仿真系统从模型到可跑框图的分层搭法3.1 MATLAB Function 还是纯 Simulink 模块Simulink 仿真系统搭法分两派一派喜欢从 Simulink 库拖 Continuous、Math Operations、Lookup Table把动力学方程用 Block 拼出来另一派像我一样动力学和控制率都写在 m 函数里Simulink 里只放 MATLAB Function 模块、积分器和信号路由。后者的维护成本低很多尤其当你要把控制器代码直接生成 C 代码时m 函数几乎可以原样移植。我的分层习惯是四个 SubsystemReference Generator给定轨迹、Controller双环 PID、Plant动力学、Sensors反馈信号。Controller 内部再分 Position Loop 和 Attitude Loop 两个子系统。层级越清晰哪个环节发散就越容易定位。散成一堆 Gain 和 Integrator 的模型波形飞了根本不知道是控制器问题还是模型问题。3.2 反馈回路和作动器限幅仿真程序“看起来能用”和“真能用”的区别没做过四旋翼仿真的人常在反馈回路里犯一个错仿真里直接把真实状态接到控制器绕过传感器模型。这样跑出来的结果非常漂亮但一旦放上真机就会发现陀螺仪噪声和延时能把姿态响应彻底拉垮。我一般会在 Sensors 子系统里加一个测量延时 Transport Delay 和零均值高斯噪声用 Random Number 模块采样步长设成与飞控主频一致比如 500 Hz即 0.002 s。作动器饱和是另一个被忽视的点。电机推力有上下限电机响应还有一阶延时。常见做法是在 Controller 输出后接一个 Saturation 模块再串一个1/(tau*s 1)传递函数模拟电机功放响应。如果省略这个环节你在仿真里把 PID 增益调到 20 都没问题但真实电机根本响应不过来表现为仿真收敛良好、真机原地抖。3.3 求解器与步长配置为什么同样的模型别人跑得飞起你这里卡死打开 Simulink 模型先看 Simulation 菜单里的 Model Configuration Parameters重点是 Solver 和 步长。四旋翼这种刚体动力学状态变化速度中等但如果控制器里有积分项和较高带宽的角速度环固定步长建议设在 1e-3 到 1e-4 秒之间。太大扛不住积分带来的高频分量太小时长虚高。我一般先用固定步长ode4四阶龙格-库塔步长1e-3跑通后再换成自动步长的ode45提速。注意如果模型里有 Transport Delay 或离散采样模块就不要用变步长 ode45离散模块和变步长配合常常会造成“过零检测”误报步长被不断缩小仿真速度慢到像死机。下面是配置项的速查表配置项推荐值说明Solverode45变步长或 ode4固定步长有离散模块用固定步长Max step size1e-2超过这个值姿态剧烈变化时误差大Fixed-step size1e-3固定步长时折中精度与速度Stop time10 或 20 秒先跑悬停验证够了再拉长时间Tolerance (相对)1e-3 到 1e-4误差要求越低越慢3.4 用船用初始化脚本喂给 Simulink 模型Simulink 仿真系统必须在运行前知道所有参数。常见做法有两种一是在模型的 PreLoadFcn 回调里写一条init_params;二是在模型工作区的 Model Workspace 里绑定脚本。我推荐前者原因是可读性好别人打开模型能看到这行调用不会满头问号地找一个叫params的变量。初始化脚本内容大致是params.m 1.5; % 总质量 kg params.g 9.81; % 重力加速度 NED 向下为正 params.l 0.25; % 机臂长度 m params.I diag([0.023, 0.023, 0.04]); % 惯量矩阵 kg*m^2 params.km 0.01; % 偏航力矩系数 N*m/N params.kp_att [12; 12; 8]; % 姿态角度环 P params.kp_rate [40; 40; 30]; % 角速度环 P params.ki_rate [1; 1; 0.8]; % 角速度环 I params.kt 1e-5; % 电机推力系数用于转速转推力这些数值不是拍脑袋的。params.I的量级可以按质量 1.5 kg 轴距 0.25 m 的常见小四轴估算如果 zip 包里给了实物样本优先用里面的参数。kt和km牵扯到电机转速到推力的映射真实情况是非线性的但仿真程序里常用线性近似因为控制器输出的本身就是归一化油门指令不是转速。4. 仿真程序跑起来初始化脚本、批量参数扫描与结果可视化4.1 用 sim() 函数把 Simulink 仿真系统从 GUI 里解放出来GUI 里点运行只适合验证单次效果。做参数扫描或调参必须用 MATLAB 脚本控制 Simulink 仿真系统。核心就两条命令set_param改模型参数sim运行仿真。% run_sweep.m % 批量运行四旋翼 simulink 模型扫描不同的 kp_rate kp_list [30, 40, 50]; for i 1:length(kp_list) params.kp_rate [kp_list(i); kp_list(i); kp_list(i)*0.75]; % 把参数写进模型工作空间 assignin(base, params, params); out sim(quad_sim.slx); % 从仿真输出对象提取姿态响应 t out.tout; att out.att.Data; % att 是 3 列phi, theta, psi % 记录超调量第一次越过目标的角度 overshoot(i) max(abs(att(:,1))); endassignin(base, ...)把变量写进基础工作区模型运行时PreLoadFcn或直接用基础工作区的params就能读到。如果out.att这种后续. 语法在你的 MATLAB 版本里报错原因是你没有给模型的 To Workspace 模块设置正确的保存格式把输出格式改成 Timeseries 或 Structure with time 后重跑即可。4.2 批量的正确打开方式parfor 与随机初始化上面的 for 循环在参数组数多时效率很低。控制变量法扫描三个参数、每个参数十五个值就要跑 3375 次单次算 0.5 秒累计快半小时。我一般把参数组合预生成好然后用parfor并行kps linspace(20, 60, 9); over zeros(size(kps)); parfor i 1:length(kps) p params; % 每个 worker 独立副本 p.kp_rate [kps(i); kps(i); kps(i)*0.75]; assignin(base, p, p); % 注意不能直接 assignin base, 要确保 worker 有独立空间 out sim(quad_sim.slx); att out.att.Data; over(i) max(abs(att(:,1))); endparfor里每个 worker 有自己的基础工作区params是用p params复制到 worker 里的避免多个 worker 读写同一个变量导致数据竞争。sim函数在parfor里能并行跑这是 MATLAB 并行计算工具箱的基本用法。如果没装工具箱就用普通 for别勉强。4.3 把仿真结果画成看得懂的图三维轨迹和时间波形仿真程序跑完没人想看满屏的向量。我把固定动作写成一对绘图脚本plot_state.m和animate_flight.m。时间波形用plot(t, eul*180/pi)会把欧拉角转成度因为人脑对度的感知远好于弧度。三维轨迹用plot3figure; plot3(pos.Data(:,1), pos.Data(:,2), -pos.Data(:,3), LineWidth, 1.5); grid on; xlabel(x (m)); ylabel(y (m)); zlabel(z (m, 向下为正)); title(四旋翼三维轨迹);-pos.Data(:,3)是把 NED 下的 z 翻转成“向上为正”的显示习惯。另外强烈建议加一张电机推力曲线图在 Simulink 里给四个电机指令各接一个 To Workspace然后画在同一张图。很多参数问题在姿态波形上看不出来但推力曲线会暴露振荡频率和饱和情况。5. 离入手三招离线调参顺序、代数环排错、外部模式5.1 双环 PID 离线调参顺序拿到现成的四旋翼 MATLAB 仿真程序别一把梭把外环位置环的 P 调到 5 就指望跟定轨迹。我的顺序是先锁死位置环给定期望角度为阶跃信号只调内环角速度环。角速度环稳了再开角度环角度环稳了再解锁位置环。每环做三次阶跃实验看超调量、调节时间和稳态误差记录成下面的表被调环阶跃值超调量调节时间(±5%)结论角速度环2 rad/s8%0.12 s合格角度环10°12%0.3 s略大降一点角度环 P位置环2 m20%2.1 s位置环 P 引入过大阻尼导致响应慢这个表是核心判断工具。如果角速度环阶跃已经明显发散那问题在params.I量级大概率是惯量设太大了力矩方向对不上。5.2 代数环排查Simulink 仿真系统里最常见的卡死原因之一是代数环。表现是仿真进度卡在 0% 不动或者报 “Algebraic loop” 警告。四旋翼模型里位置环和姿态环相互耦合如果控制器输出直接通过动力学反馈回到控制器输入中间没有任何延时单元就会形成代数环。解决方法是在反馈路径上串一个Unit Delay或Memory模块模拟飞控的离散采样特性。加延时后注意离散延时意味着控制周期也就是你的仿真步长应该等于真实飞控循环时间不然相位滞后和实际不符。5.3 外部模式与代码生成仿真到半实物衔接的最后一米如果 zip 包里的仿真程序最终要驱动真实飞控终极验证方式是用 Simulink 的 External Mode 或 SILSoftware in the Loop。常见做法是先在 Simulink 里把控制器子系统配置为 C 代码生成目标用 Embedded Coder 生成 C 代码再编译进飞控固件的控制线程。生成前有两个必查项一是控制器中不能有persistent变量吗允许有但 MATLAB Coder 对它有限制最好改成coder.extrinsic(zero)或直接用传入的dt作为采样周期避免代码生成后变量初始化不确定。二是代码生成后仿真数值和原模型数值对不上的常见原因是除法顺序I \ M和inv(I)*M在代码生成里有时产生不同的浮点误差统一写成前者。最后再强调一个容易踩的点Simulink 里显示的simout数据在 External Mode 下很可能不再实时更新因为代码被部署到目标机后 Simulink 只是通过通信协议收数据存储格式和离线仿真不同要检查你选的 To Workspace 模块的输出格式在外部模式下是否兼容不兼容就换成通过网络输出的 Dashboard 或通过串口日志方式回收。把这一层想清楚四旋翼仿真程序才算真正闭环。本文还有配套的精品资源点击获取
返回列表