
简介四旋翼飞行器MATLAB仿真程序包含Simulink仿真系统与配套脚本是一套覆盖建模仿真全流程的学习资源适合自动化、航空航天等专业学生用于毕业设计、期末大作业或课程设计。资源内置悬停控制、轨迹跟踪、线性编队与圆形编队等典型实验并提供状态转换、姿态解算、四元数与欧拉角转换、3D轨迹绘制等基础函数便于从动力学模型到控制验证逐层理解。压缩包共23个文件以m脚本15个和slx模型4个为主辅以交互式mlapp应用、实验报告PDF及说明文档整体仅5.27MB目录清晰、轻量易部署。源码经本地编译运行验证评审分达98分难度适中内容经助教审定可直接用于复现实验或二次开发。当前已有135人学习使用适合需要快速搭建四旋翼仿真环境、深入理解控制算法的学习者。1. 把四旋翼仿真跑在MATLAB里之前先想清楚两件事一份四旋翼仿真程序拿到手先别急着双击.slx点 Run。十个有九个第一次跑会遇到同一件事模型能走、曲线在跳但改任何一个参数都要翻遍初始化脚本、模块对话框和传感器噪声模块三层地方。标题里既然写了“含说明和报告”真正值钱的不是那几张曲线图而是参数、模型、控制律、验证脚本和报告之间那条能复现的链路。这套仿真本质上要解决两类问题一是验证姿态和位置控制算法在连续时间模型下收不收敛二是给后续硬件移植留一套可回归测试的基线。适合做飞控算法验证、课程设计和课题预演的人照着下面的顺序能自己搭出一套带说明文档的完整仿真程序而不只是会拖动 Simulink 模块。2. 四旋翼Simulink建模初始化脚本、动力学方程与子系统划分2.1 先把参数和工作区的关系定下来一份四旋翼仿真程序拿到手第一件事不是打开模型文件而是找到初始化脚本把它跑一遍。常见做法是init_quad.m与quad_model.slx放在同一目录Simulink 模型里所有m、I、kF这类符号都从 MATLAB 工作区读取。这样参数修改只动脚本不需要进模型一个个点开模块对话框也不会出现“报告里写了一组参数、模型里又是另一组”的情况。% init_quad.m 四旋翼与仿真公共参数 m 0.52; % 整机质量(kg) l_arm 0.20; % 旋翼中心到机体质心距离(m) I diag([0.0021, 0.0021, 0.0038]); % 转动惯量(kg*m^2) kF 2.5e-5; % 推力系数(N/(rad/s)^2) kM 6.5e-7; % 反扭矩系数(N*m/(rad/s)^2) g 9.806; % 重力加速度(m/s^2) wm 1500; % 电机最高转速(rad/s) Ts 0.002; % 仿真步长(s)也用于外部模式这些参数里kF和kM没有统一标准值取决于桨径、电机 Kv 和电压。上面给的是 0.5 kg 级小四轴的量级做仿真验证够用如果要对应真实机架用静拉测试台测一下悬停油门反过来推kF比从说明书上抄更可靠。wm后面控制器限幅、报告里电机裕量分析都要用不要只写在一个子系统里。提示模型里报 Undefined function or variable 时多数情况是初始化脚本没跑或 Model Properties Callbacks InitFcn 里没写run(init_quad.m)。直接点 Run 之前先跑一次脚本是这套仿真程序最基本的使用约定。2.2 控制效率矩阵转速的平方才是线性输入四旋翼四个螺旋桨产生的升力是F_i kF * w_i^2反扭矩是M_i kM * w_i^2。注意这里是转速平方不是转速本身所以电机模型后端一定要接u^2运算控制分配矩阵也工作在平方域。把四个转速平方映射成总推力、滚转力矩、俯仰力矩和偏航力矩可以写成矩阵形式。% 控制效率矩阵B * (omega.^2) [Fz, tau_phi, tau_theta, tau_psi] % 电机编号1前 2右 3后 4左相邻电机转向相反 B [ kF, kF, kF, kF; -kF*l_arm, kF*l_arm, kF*l_arm, -kF*l_arm; kF*l_arm, kF*l_arm, -kF*l_arm, -kF*l_arm; -kM, kM, -kM, kM ];第二行是滚转力矩第三行是俯仰力矩第四行是偏航力矩。正负号来自电机位置和转向如果你换成 X 型机架第二三行的符号要重新推不能直接套。控制器输出四维向量[Fz; tau_phi; tau_theta; tau_psi]后用omega sqrt(pinv(B) * u)反解转速再按[0, wm]限幅。这里用伪逆而不是inv(B)是因为当某一轴饱和时伪逆分配出来的四路转速不容易出现负值后面想加约束优化也好扩展。这里有一个常见坑限幅要在omega.^2之后做。对转速本身限幅会引入平方根分支悬停点附近四个转速都不是零如果先限幅再平方相当于把非线性算子放进了控制回路模型容易发飘。2.3 状态方程机体内力矩转到惯性系加速度被控对象状态常用 12 维向量[x y z phi theta psi vx vy vz p q r]其中phi theta psi是欧拉角p q r是机体系角速度。Simulink 里的 Plant 子系统可以是一个 MATLAB Function输入是四个电机转速和当前状态输出是加速度和角加速度。function [acc, angacc] quad_plant(omega, state, p) omega min(max(omega, 0), p.wm); % 转速限幅 w2 omega.^2; u p.B * w2; % 总推力和三轴力矩 Fz u(1); tau u(2:4); pqr state(10:12); angacc p.I \ (tau - cross(pqr, p.I * pqr)); % 角加速度 phi state(4); th state(5); ps state(6); R [cos(th)*cos(ps), sin(phi)*sin(th)*cos(ps)-cos(phi)*sin(ps), cos(phi)*sin(th)*cos(ps)sin(phi)*sin(ps); cos(th)*sin(ps), sin(phi)*sin(th)*sin(ps)cos(phi)*cos(ps), cos(phi)*sin(th)*sin(ps)-sin(phi)*cos(ps); -sin(th), sin(phi)*cos(th), cos(phi)*cos(th)]; acc [0; 0; -p.g] R * [0; 0; Fz/p.m];旋转矩阵是“机体系到惯性系”的 ZYX 顺序符号要和 MATLAB 自带angle2dcm的输入顺序区分开。angacc里的cross(pqr, I*pqr)项不能省大速率滚转时进动力矩会让姿态收敛变慢悬停附近可以忽略但既然做仿真留着更好。欧拉角在theta ±90°时会奇异因为由p q r到欧拉角速率的变换里有1/cos(theta)。如果这个仿真只做悬停和小角度机动用欧拉角没问题报告里也好解释要是做翻转或大机动状态变量换成四元数只在输出端转回欧拉角。2.4 Simulink子系统拆法参考、控制、分配、动力、传感常见做法是把模型拆成五个子系统数据流是单向的子系统输入输出关键点Reference期望位置、期望偏航位置指令、偏航指令梯形限制器防超调Position Control位置、速度反馈期望姿态、期望油门外环输出不直接给电机Attitude Control期望姿态、角速度反馈总推力、三轴力矩内环带宽高于外环Motor Mixer推力、力矩四个转速指令pinv(B)后限幅Dynamics四个转速12维状态连续积分Sensor状态测量值加噪声和零阶保持这样拆的依据是控制信号流向参考量进控制器控制器出力矩力矩经分配变转速转速进动力学状态再回到控制器。按这个方向画模型不会出现反馈绕不开的问题。代数环是这类模型里最容易把新人卡住的问题。路径是力矩tau影响p q rp q r又会被控制器拿来算力矩同一仿真步内形成闭环。Simulink 检测到代数环会用迭代求解结果就是仿真卡顿或者结果抖动。两种常见处理一种在动力学输出端加一个1/(0.02*s1)的滞后等效于电机响应延迟另一种在反馈通道里放一个 Unit Delay模拟传感器采样延迟。两种都在物理上有意义不是作弊。3. 四旋翼姿态控制与simulink PID整定的正确顺序3.1 内外环结构姿态环要快于位置环一个数量级从控制结构看四旋翼位置控制不直接给电机。位置环先算期望加速度再通过小角度近似把水平期望加速度映射成期望横滚和俯仰角姿态环跟踪这个角度。两层都常用 PID但顺序必须是先整内环再整外环。内环带宽一般取外环的 5 到 10 倍外环每次发出姿态指令时希望内环已经收敛外环眼中看到的才是一个近似一阶惯性的对象。如果带宽倒挂横滚角还没建立起来位置环已经输出饱和最后表现就是起飞瞬间左右晃动、高度掉一块再拉回来。映射关系一般写为theta_des ≈ a_y_des / g、phi_des ≈ -a_x_des / g符号取决于坐标系定义报告里必须写清楚否则换个人看模型会怀疑控制方向反了。3.2 从Ziegler-Nichols到可用的PID初值我一般不用工具箱里的自动整定先在 Simulink 里把 PID 控制器的积分和微分关掉P 从 0.5 开始倍增给一个小幅角度阶跃直到系统出现等幅振荡。记录临界增益Ku和临界振荡周期Tu再查表控制器KpKiKdP0.5 Ku--PI0.45 Ku1.2 Ku / Tu-PID0.6 Ku1.2 Ku / Tu0.075 Ku * Tu四旋翼电机转速有上限Ziegler-Nichols 查表值偏激进工程上一般乘 0.7 再试。先对 roll 和 pitch 调固定这两个环的参数后再去调偏航环偏航环的响应可以慢一些因为偏航角不影响水平位置的直接控制。整定过程中输出限幅先放开等参数基本可用后再收回到电机能承受的范围内。3.3 Simulink PID 模块容易忽略的三个参数打开 Simulink 里的 PID Controller 模块大多数人只改P、I、D三个数实际有三项会影响仿真结论。第一项是 Filter coefficientN。纯微分项会把陀螺仪噪声放大所以 PID 模块内部对微分项做了低通近似N越大越接近纯微分越小滤波越强。四旋翼仿真里取 30 到 80 比较常见先设 50看加了噪声之后的高频抖动再调。第二项是 Anti-windup 方式和积分限幅。把 Anti-windup method 选成 clamping并把 Integrator output limit 设成和控制输出限幅一致。否则大角度偏差时积分一直累计转速已经饱和误差恢复后积分还没退下来表现出来就是大幅超调而且每次手动重跑结果还不一样。第三项是 Time domain 选择。纯仿真用 Continuous后续做 C 代码生成时改成 Discrete并填上Ts。两套参数不能混用连续 PID 的积分和微分在离散化后需要重新验证。提示Control System Toolbox 里有pidtune可以对线性化模型自动算初值再用自动参数做 Ziegler-Nichols 的起点。但不要直接对非线性 Simulink 模型做反复自动整定饱和、噪声和代数环会让优化方向失去物理意义。3.4 内外环初值表与批量调参脚本下面是一组 0.5 kg 级四旋翼的 PID 初值适合作为起点不等于最终值控制环PID输出限幅姿态内环 roll/pitch4.00.022.5±0.6 N*m偏航环2.00.041.0±0.3 N*m位置外环 x/y1.20.020.6姿态指令 ±0.3 rad高度环2.50.051.0推力 0 ~ mg1.6调参时要避免每次手动改模块对话框写一个批量脚本循环跑clear; init_quad; kp_list [3.2 3.6 4.0 4.4 4.8]; for i 1:numel(kp_list) set_param(quad_model/Attitude_PID/P, Gain, num2str(kp_list(i))); out sim(quad_model.slx, StopTime, 5); roll out.logsout.getElement(roll).Values.Data(:); ref out.logsout.getElement(roll_ref).Values.Data(end); os (max(roll) - ref) / ref * 100; fprintf(kp%.1f overshoot%.2f%%\n, kp_list(i), os); endset_param的参数路径要对应实际模型里的模块名Gain直接写数值是为了循环方便更规范的做法是把控制器参数定义成工作区变量脚本里用assignin改。注意logsout要先在 Configuration Parameters Data Import/Export 里勾选 Signal logging否则脚本会报找不到数据。用脚本批量调参时把 Simulink 窗口最小化能省掉大量界面刷新时间。4. 仿真验证阶跃指标、噪声注入、外部模式与C代码生成4.1 用 logsout 而不是 Scope 收集数据Scope 看形状可以写说明文档时要数字必须用脚本取数。仿真结束后从out.logsout按信号名取数据而不是在模型里加一堆 To Workspace 模块。out sim(quad_model.slx, StopTime, 15); t out.tout; roll out.logsout.getElement(roll).Values.Data; ref out.logsout.getElement(roll_ref).Values.Data;logsout的信号名要和模型里 Signal Logging 勾选的名字一致命名用roll_ref这种下划线风格不要带空格和中文。取数后先看长度如果t和roll长度不一致检查一下是否开了变步长四旋翼仿真建议固定步长ode4 步长Ts数据分析和报告图表都好处理。4.2 用函数算指标而不是肉眼看图报告里需要的超调量、最终误差、收敛时间写一个小函数一次性算出来比在图上手动取点可靠。function [os, ess, ts] quad_metrics(t, y, r) % 阶跃输入的时域指标 os max(y) - y(end); % 超调量 ess (r - y(end)) / r; % 归一化最终误差 idx find(abs(y - r) 0.02 * abs(r), 1); ts t(idx); % 2% 收敛时间 endos用最大值减终值比减参考值更抗噪ts的 2% 准则比 5% 严格四旋翼报告里建议统一用 2%否则不同人写出来的收敛时间没有可比性。还可以用periodogram画悬停状态的功率谱看低频段有没有 1 Hz 左右的振荡峰那是姿态环带宽不足的典型信号。4.3 传感器噪声注入与带宽取舍模型里加传感器噪声用的模块是 Band-Limited White Noise放在 Sensor 子系统里。这个模块的参数Noise power不是直接填标准差而是填功率谱密度。如果希望输出序列的标准差是sigmaNoise power要设成sigma^2 * Ts。传感器标准差Noise power 参数陀螺仪0.01 rad/s1e-4 * Ts加速度计0.02 m/s^24e-4 * Ts气压高度0.05 m2.5e-3 * Ts加噪后姿态内环的 D 项会明显放大噪声这正是前面把滤波器系数N调到 3080 的原因。如果加噪后 roll 通道高频抖动超过 0.1 rad把 D 降 20% 比继续加大N更有效因为N太小时微分项相位滞后会抵消阻尼作用。4.4 外部模式与C代码生成、carsim联合仿真Simulink 外部模式是模型部署前的一个重要中间步骤。配置好 Communication Interface 后模型可以运行在目标硬件上Simulink 端实时观测状态和修改参数不需要重新生成完整工程。对这套四旋翼仿真程序来说外部模式验证的不是动力学模型而是控制子系统。如果要把控制律迁移到飞控用 Simulink Coder 配置生成 C 代码。生成前先按 CtrlD 执行 Update Model完成模型检查和代码级检查System target file 选ert.tlc。生成时只选控制器和信号处理子系统不要把整个被控对象动力学模型生成进去动力学模型只是桌面验证用的。和 Carsim、AMEsim 这类外部工具做联合仿真的思路也类似把 Simulink 侧的被控对象替换成外部模型通过 S-Function 或标准接口传递 double 类型信号。这里最容易踩的坑是接口信号命名不统一建议所有联合仿真接口都在模型里单独建一个 InBus/OutBus字段名形成版本记录。5. 把仿真程序整理成说明与报告章节、自动导出和复现技巧5.1 报告不是抄模型而是把“为什么”写清楚写报告时最容易出现的问题是拿模型截图贴一遍读者看完还是不知道参数从哪来。我一般会把说明文档固定成几个章节和模型文件一一对应报告章节内容对应文件系统建模坐标系定义、B 矩阵、状态方程init_quad.m、quad_plant.m参数表所有物理量和来源init_quad.m控制律设计内外环结构、PID 整定过程batch_tune.m验证结果阶跃、噪声场景图表figs/目录复现说明MATLAB 版本、求解器、运行顺序README.md复现说明里一定要写执行顺序先run(init_quad.m)再open_system(quad_model.slx)最后sim或点 Run。这三步缺一不可。另外把 MATLAB 版本和求解器设置写进报告因为不同版本对 PID 模块默认参数和exportgraphics的支持不一样。5.2 一次导出全部图和参数扫描exportgraphics 与 Fast Restart仿真结果图不要手工截图用脚本统一导出统一把分辨率设成 300 dpi 放到report/figs/。mkdir(report/figs); figs findobj(Type, figure); for i 1:numel(figs) ax get(figs(i), CurrentAxes); name get(get(ax, Title), String); if isempty(name) name sprintf(fig%d, i); end exportgraphics(figs(i), fullfile(report, figs, [name .png]), ... Resolution, 300); endexportgraphics在 R2020a 之后都支持比print更适合导出带坐标区内容的图线条不会变细。注意图窗标题要提前用title()写好文件名就自动对应信号名。做参数扫描时打开 Simulink 的 Fast Restart改参数不会重新编译模型适合多组 PID 对比。set_param(quad_model, FastRestart, on); for k kp_list assignin(base, Kp_att, k); out sim(quad_model.slx, StopTime, 5); % 收集时域指标 end set_param(quad_model, FastRestart, off);Fast Restart 下改了模型结构参数不会生效所以只调 PID 数值时用一旦改动模型结构先关掉再重跑。把 Fast Restart 和exportgraphics放进同一个run_all.m参数扫描完成时附录图已经自动生成剩下的事是把 README 里的复现步骤再对一遍。本文还有配套的精品资源点击获取