简介这份Matlab四旋翼仿真资源面向毕业设计、课程设计及无人机控制方向初学者提供可运行的四旋翼建模、控制算法与仿真验证工具。压缩包共280个文件主要包含125个m脚本用于核心算法与仿真逻辑37个csv数据文件对应飞行位置、速度等实验结果记录13个mat与6个fig文件可用于状态数据读取、曲线查看与结果可视化另有msg、srv等辅助文件便于功能扩展与二次开发。整体约6.63MB结构清晰、轻量易用。目前已有132人学习下载所有源码经严格测试可直接运行适合需要快速搭建仿真环境或完成课程项目的人群。通过本资源可获得完整的四旋翼模型、控制算法实现、实验数据及可视化结果帮助理解飞行动力学与控制原理同时也可以作为毕业设计、课程论文或Matlab编程练习的重要参考素材。1. 毕业设计选四旋翼仿真Matlab这套资源在解决什么做毕业设计选四旋翼方向的人大部分时间其实不是耗在写代码而是耗在“模型跑不起来”和“控制器没法调”这两件事上。Matlab平台的优势在于生态完整Simulink里搭控制回路、改参数、看波形都比纯手写C直观得多。如果手里有一套Matlab四旋翼仿真框架把动力学建模、控制律、Simulink仿真串成闭环工作量能省掉一半以上。这套资源正是围绕毕业设计场景打包的从牛顿-欧拉方程到串级PID再到可视化输出基本覆盖了课设和毕设最常见的技术路线。适合那些刚接触四旋翼建模、又想在答辩前拿到稳定仿真结果的同学。我自己当年就是从类似的框架入手的踩了不少坑才把模型从“飞出去”改成“停得住”。2. 四旋翼动力学建模从刚体运动学到状态空间方程2.1 坐标系定义与欧拉角约定四旋翼仿真第一课不是写方程而是先把坐标系钉死。我见过太多成品代码里俯仰角正负号反了或者欧拉角顺序用了两种结果模型一跑就翻跟头。推荐的做法是地面系定义为NED北-东-地或ENU都可以但机体轴系必须右手定则。这里以多数教材采用的ZYX欧拉角为例先偏航角ψ再俯仰角θ最后横滚角φ。如果把机体坐标系到地面坐标系的旋转矩阵写成一个函数Matlab里通常是function R rotation_matrix(phi, theta, psi) Rz [cos(psi) -sin(psi) 0; sin(psi) cos(psi) 0; 0 0 1]; Ry [cos(theta) 0 sin(theta); 0 1 0; -sin(theta) 0 cos(theta)]; Rx [1 0 0; 0 cos(phi) -sin(phi); 0 sin(phi) cos(phi)]; R Rz * Ry * Rx; end这个矩阵的作用是把机体坐标系下的力、速度转换到地面系。注意Matlab的三角运算默认弧度初始化时千万别把角度直接塞进来。我自己习惯在工程入口写一个转换函数deg2rad统一处理。为了验证坐标系没写错可以做个简单测试输入φ0θ0ψpi/2旋转矩阵应该把机体x轴转到地面y轴方向。如果得到的结果是x轴仍然指向原方向说明旋转矩阵顺序错了。这种小测试能省下后面排查的时间。另外还有一个容易忽略的点欧拉角约定必须与你的控制期望一致。如果你在控制器里把期望俯仰角定义为绕机体y轴旋转而旋转矩阵里用的是ZYX顺序那外环计算出来的期望角速度就会对不上。建议在初始化脚本里用注释写明欧拉角顺序并把这个注释同步复制到控制器模块里。2.2 牛顿-欧拉方程的Matlab实现四旋翼的刚体运动学分为平动和转动两部分。平动方程基于牛顿第二定律转动方程基于欧拉方程。这里忽略桨叶挥舞、气动阻力高阶项保留主要力与力矩每个旋翼产生的推力T_i k * w_i^2k是推力系数力矩包括推力产生的滚转/俯仰力矩、反扭矩产生的偏航力矩把这套方程写成Matlab函数输入状态向量与四个电机转速输出状态导数function x_dot quad_dynamics(x, u, param) % x [x, y, z, vx, vy, vz, phi, theta, psi, p, q, r] % u [w1^2, w2^2, w3^2, w4^2] 电机转速平方 g param.g; m param.m; k param.k; % 推力系数 d param.d; % 反扭矩系数 Ixx param.Ixx; Iyy param.Iyy; Izz param.Izz; l param.l; % 机臂长度 phi x(7); theta x(8); psi x(9); p x(10); q x(11); r x(12); T1 k * u(1); T2 k * u(2); T3 k * u(3); T4 k * u(4); T_total T1 T2 T3 T4; % 平动加速度地面系 R rotation_matrix(phi, theta, psi); acc [0; 0; -g] R * [0; 0; T_total] / m; % 力矩机体系X型布局 Mx l * k * (u(1) - u(3)); My l * k * (u(2) - u(4)); Mz d * (u(1) - u(2) u(3) - u(4)); % 角加速度欧拉方程 p_dot (Mx (Iyy - Izz) * q * r) / Ixx; q_dot (My (Izz - Ixx) * p * r) / Iyy; r_dot (Mz (Ixx - Iyy) * p * q) / Izz; % 欧拉角速率注意 theta -90 度时奇异 phi_dot p sin(theta) * (q * sin(phi) r * cos(phi)) / cos(theta); theta_dot q * cos(phi) - r * sin(phi); psi_dot (q * sin(phi) r * cos(phi)) / cos(theta); x_dot [x(4); x(5); x(6); acc(1); acc(2); acc(3); ... phi_dot; theta_dot; psi_dot; p_dot; q_dot; r_dot]; end这里有几个关键点。首先u是电机转速的平方而不是转速本身。这样写的好处是推力与转速的关系变成了线性系数k后续混控矩阵可以直接用。如果你拿到的数据是转速表的记得先平方再做输入。其次力矩公式里的符号。X型布局常见的电机编号是1号右前逆时针2号左后顺时针3号左前逆时针4号右后顺时针。滚转力矩由x轴两侧推力差产生所以我这里写成l * k * (u(1) - u(3))。如果你的布局是十字型符号需要重新推一遍。这个符号错了模型会表现为“给一个正的期望滚转角飞机却往反方向滚”。第三欧拉角速率公式里有一个1/cos(theta)当俯仰角接近90度时数值会发散。毕设阶段如果只做悬停和轨迹跟踪通常不会碰到大角度但如果你要做翻转动作必须换成四元数表示。这个话题在第五章排查里还会展开。2.3 状态空间模型与接口设计为了和控制器衔接把动力学整理成状态空间形式。我们选状态向量X [x, y, z, vx, vy, vz, phi, theta, psi, p, q, r]控制输入U [u1, u2, u3, u4]。将2.2的导数函数写成x_dot f(x, u)。注意这个f是非线性的不是线性状态空间。在悬停点附近线性化后可以写成矩阵A、B用来做LQR或者稳定性分析。在资源包里通常会同时给出非线性模型和线性化矩阵。线性化方法可以直接用Matlab的linmod或符号工具箱的jacobian。我们在工程里更常用的是在悬停姿态φθψ≈0、速度为零、推力等于重力处求数值雅可比% 在悬停点求雅可比 x0 zeros(12,1); u0 [m*g/4; m*g/4; m*g/4; m*g/4]; A jacobian((x) quad_dynamics(x, u0, param), x0); B jacobian((u) quad_dynamics(x0, u, param), u0);注意这里jacobian需要Symbolic Math Toolbox如果没装可以用数值差分法自己算。这一步不是必需的但做了之后你对模型的固有频率、耦合关系会更有底。比如A矩阵里位置和姿态之间的耦合项是哪些直接就能看出来。接口设计方面我习惯把quad_dynamics这个函数封装成Simulink的MATLAB Function模块输入是状态向量和控制量输出是导数。这样后续做S函数或硬件在环时函数体可以原样复用。注意输入输出都用double类型Simulink里默认是double别改成single否则精度损失会导致漂移更严重。状态空间模型还有另一个用途设计观测器。比如很多毕设题目要求“基于状态观测器的姿态控制”这时你需要在仿真里引入状态重构环节。把非线性模型作为被控对象把线性模型作为观测器的标称模型这是标准做法。资源包里如果给了A、B矩阵你要能说明白它是在哪个平衡点线性化的答辩老师最常问的就是这个。3. PID控制器设计与调参串级结构在Simulink里的落地3.1 串级PID结构内环角速度、外环姿态四旋翼控制最常用的是串级PID内环控制角速度(p,q,r)外环控制姿态角(φ,θ,ψ)。外环输出期望角速度给内环内环输出期望力矩再通过混控矩阵生成四个电机的转速指令。这种结构的好处是内环能快速抑制扰动外环可以慢一点抗干扰能力强。在Simulink里可以把内环和外环分开封装成两个子系统便于调试。资源包里一般会有现成的参考设计但最好自己动手连一遍否则答辩时被问到底层逻辑容易卡壳。我见过不少人拿到的模型只给了封装好的子系统双击进去看到一堆常量根本不知道每个数是什么含义。这种黑匣子是毕设答辩的大忌。串级PID与单环PID的区别在于单环直接对姿态角误差输出力矩相当于把角速度环忽略了。角速度环的存在可以极大提升阻尼特性让姿态响应更平滑。对于四旋翼这种开环不稳定的系统强烈建议使用串级结构。如果你发现资源包里给的是单环PID建议改成双环哪怕只加一个P控制器也会有明显改善。3.2 控制器代码与参数初始化典型的姿态控制器可以写成Matlab函数使用位置式PID。下面给一个内环加外环实现的简化版function tau attitude_controller(phi, theta, psi, p, q, r, phi_d, theta_d, psi_d, param) % 外环姿态角误差 - 期望角速度 kp_phi param.kp_phi; kp_theta param.kp_theta; kp_psi param.kp_psi; p_d kp_phi * (phi_d - phi); q_d kp_theta * (theta_d - theta); r_d kp_psi * (psi_d - psi); % 内环角速度误差 - 力矩 kp_p param.kp_p; kp_q param.kp_q; kp_r param.kp_r; Mx kp_p * (p_d - p); My kp_q * (q_d - q); Mz kp_r * (r_d - r); tau [Mx; My; Mz]; end这里只写了P控制为了展示骨架。实际工程中内外环都要加积分和微分。积分项用来消除稳态误差但四旋翼悬停时积分累积很容易超调甚至振荡所以通常加抗饱和限幅。微分项放在测量噪声上会放大抖动一般用低通滤波后的角速度。参数初始化是一个复杂的活。给一个初始参数表参数典型初始值说明kp_phi/kp_theta5~8外环比例过大则姿态回正过快容易振荡kp_psi3~5偏航响应可以慢一些kp_p/kp_q/kp_r0.1~0.3内环比例单位是力矩比角速度ki0.01~0.05积分系数先设0稳定后再加kd0.01~0.05微分系数先设0这些值不是通用的取决于模型里的惯性矩和机臂长度。我用的是一个某小型四旋翼的参数质量1.1kg臂长0.22mIxx约0.012kg·m²。如果你手里是其他机型参数得重新标定。比如质量变成5kg推力系数也要成倍调整控制器比例系数可能就完全不对了。混控矩阵怎么接姿态控制器输出的tau是力矩需要转换成四个电机的转速平方。常用的是% 混控矩阵力矩 总推力 - 电机平方 mix [1 1 1 1; 1 0 -1 0; 0 -1 0 1; 1 -1 1 -1] / 4; % 实际要根据布局正负号调整这个矩阵的意思是把[T_total, Mx, My, Mz]映射到[u1,u2,u3,u4]。混控矩阵的系数需要与动力学里的力矩公式完全对应否则控制器输出对的电机响应却是错的。一个简单的验证方法把四个转速设为某个只产生滚转力矩的组合跑一下动力学函数看输出的力矩是否与期望一致。3.3 调参顺序与响应曲线判读调参常用的顺序是先内环后外环。把外环踩掉直接给角速度阶跃调内环P让响应快且无振荡再加D抑制超调最后加I消除稳态误差。内环满意后再闭合外环调外环P让姿态角阶跃响应达到期望水平。具体操作时我习惯在Simulink里把参考输入设为一个幅值0.3 rad/s的阶跃信号。注意不要太大角速度环的线性区间有限。观察角速度响应曲线如果出现持续等幅振荡说明P过大如果响应慢吞吞说明P太小或者模型里惯性参数设大了。内环调好后外环姿态阶跃可以给5度的阶跃约0.087 rad。如果外环响应超调量超过20%适当减小外环P。四旋翼的滚转和俯仰通道在参数上可以对称偏航通道单独调因为偏航转动惯量通常更大。这里有一个血泪经验调参时千万别同时改多个参数否则振荡来源根本定位不到。我习惯每次只改一个系数跑一轮仿真记录下超调和调节时间形成一张表再判断。这个习惯救了我很多次。调参时还容易犯一个错误在仿真结果里同时看位置和姿态发现位置漂移就以为是控制器问题。实际上位置漂移很可能是初始姿态角没设为零或者模型的重力方向写错了。我建议先只观察姿态角把姿态角稳定性调好之后再观察位置。如果位置还在缓慢漂移检查动力学函数里平动加速度的计算特别是R * [0;0;T_total]这一项是否确实是重力加推进力的方向。数据判读上我一般记录三个指标上升时间从10%到90%、超调量、调节时间进入±2%误差带。把这些指标做成一张表每次调参后更新。这样做的好处是你能直观看到参数变化的趋势而不是凭感觉瞎试。4. 仿真环境初始化与配置从零搭起一个可复现的四旋翼工程4.1 工程目录与脚本结构一个毕业设计要交的代码不只是模型文件还要别人打开就能跑。我见过很多同学把脚本和模型堆在一个目录里命名叫final_new2.m这种交付方式在评阅时很吃亏。推荐的结构是project/ ├── model/ │ ├── quad_dynamics.m │ ├── rotation_matrix.m │ └── attitude_controller.m ├── sim/ │ ├── main_sim.slx │ └── sim_config.m ├── data/ │ ├── results.mat │ └── logs/ ├── doc/ │ └── 说明文档.md └── init.minit.m是总入口里面设置工作区变量并调用sim_config.m。这套资源里通常会有一个类似的模板但目录名不一定一样。强烈建议你保留自己的README写清楚每个文件的作用和运行方法。我看到很多人的项目文件之间是隐式依赖比如某个脚本里调用了另一个脚本里定义的变量但别人找不到这个脚本。这种情况在评阅时容易被扣分。正确做法是所有参数都在init.m中定义模型和控制器脚本只从param结构体读取。这样任何人拿到工程先运行init.m再打开main_sim.slx就能复现结果。另外init.m最好把addpath加好把model和sim目录都包含进去避免依赖当前工作目录。4.2 初始状态与物理参数配置初始化脚本里最重要的就是物理参数。用param结构体统一存放质量、惯性矩、推力系数等不要散落在各个脚本里。例如param.g 9.81; param.m 1.1; param.l 0.22; param.k 1.6e-5; % 推力系数单位 N/(rad/s)^2 param.d 2.3e-7; % 反扭矩系数 param.Ixx 0.012; param.Iyy 0.012; param.Izz 0.023;这些系数怎么来最靠谱的是从旋翼叶素理论估算但毕设阶段常用的是查同尺寸机型的公开参数再微调。注意推力系数是谐振的单位要统一成N/(rad/s)^2不要和rpm混用。如果模型用转速平方作为输入初始化时就必须保证单位一致。我一个同学把转速单位用成了rpm结果推力差了大约1000倍模型直接变成自由落体。初始状态配置也很关键。通常需要设置初始位置x0、速度、姿态角。我建议把初始姿态角设为零向量方便验证控制器是否能让飞机从零姿态稳定悬停。如果初始姿态角设了非零值调试时你分不清是控制器纠偏有效还是初始条件导致的响应。另外初始高度可以设为1米这样能直观看到高度方向的响应。4.3 数值积分步长与仿真时长设置Simulink里求解器设置是门玄学。固定步长和变步长结果差异很大四旋翼这种强耦合系统我一般用固定步长ode4四阶龙格库塔步长1e-3秒。为什么不用变步长因为变步长在刚度大的情况下会自动缩小步长导致仿真时间不可控而且不利于结果复现。步长太大控制更新会变慢步长太小仿真时间成倍增加。1e-3秒是一个工程上比较稳妥的折中。仿真时长方面悬停仿真通常设10秒就够了你看姿态稳定过程完全够。如果是做轨迹跟踪设60秒。但注意Simulink的Stop Time单位是秒别和采样时间搞混。有时候仿真跑到一半曲线不更新了查看是不是Stop Time设成了0。还有一个容易忽略的点控制器和模型如果都放在连续时间域里仿真没问题但如果用离散PID必须设置采样周期比如200Hz。我的建议是先把控制器也做成连续时间跑通后再改成离散版本。很多成品模型直接给离散控制器同学拿去后不知道采样周期是多少结果控制效果完全不一样。另外Simulink的代数环问题在四旋翼仿真里也容易出现。如果你把控制器输出直接反馈到动力学输入而动力学输出又直接决定控制器输入中间没有延迟就会形成代数环。解决方案是给反馈信号加一个Memory模块或单位延迟打破代数环。我见过不少人在初始化时把Save options设置为timeseries但代数环导致的错误提示却看不懂这里顺便提一句。5. 四旋翼仿真的常见问题与排查失稳、漂移与数据跳变5.1 模型失稳现象、原因与解决现象仿真一开始姿态角快速发散几秒后数值变成NaN控制完全失灵。原因最常见的是推力系数或质量不匹配导致悬停推力小于重力飞机直接坠地或者控制器正负号反了输出的是正反馈。我见过一个案例控制器输出力矩的正负号和动力学里力矩公式的正负号相反结果姿态误差越大控制器给的力矩越让误差变大1秒内就翻了。解决先把控制器拆掉给四个电机一个固定输入使总推力等于重力看模型是否能稳定悬停。如果连这个开环测试都发散说明动力学函数或者坐标系有bug。我通常会在初始化脚本里算一下u0 m*g/(4*k)用这个值去测试开环。如果开环下飞机垂直方向保持不动说明平动模型正确如果姿态角开始翻动说明转动模型或旋转矩阵有误。还有一个常见原因是仿真步长太大。我曾经把固定步长设为0.01秒结果显示角速度在每步之间变化剧烈模型发散。对于强非线性系统步长需要小到能捕获系统动态。如果发散时曲线出现锯齿状优先检查步长。5.2 姿态角漂移现象、原因与解决现象控制器明明让期望姿态角为零但实际姿态角长时间稳定在某个小角度比如俯仰角一直偏着2度。原因可能是欧拉角速率公式中混入了错误项或者机体质心和推力中心不重合建模时没考虑。也可能是传感器噪声模型导致的积分漂移但仿真里一般没有传感器所以重点查刚体方程。解决在quad_dynamics函数里单独检查平动和转动部分。比如把转动部分单独跑固定位置不变输入一个常值力矩看角加速度是否符合预期。另外检查欧拉角速率公式是否有1/cos(theta)项缺失。我在某个版本里就是把phi_dot直接写成了p导致大角度时完全对不上。还有一类漂移来自控制器积分项。如果外环积分系数调得过大悬停时一点点姿态误差会被积分放大形成一个小角度偏移。排查方法是把控制器里的积分项全部置零看姿态是否回到零。如果回到零说明是积分参数问题如果仍然有偏移则回到模型检查。5.3 控制量饱和与执行器限幅现象给定的期望速度很大控制器输出力矩也很大但电机转速被限幅后实际出力不够飞机响应滞后。仿真里表现为姿态角回正速度远慢于预期甚至出现缓慢的低频振荡。原因真实电机转速不可能无限大模型里如果不加限幅仿真数据看起来正常但换到实物就会翻车。控制量饱和后积分项继续累积导致“积分饱和”现象。积分饱和会让控制器在退出饱和后仍然保持一个很大的积分输出引起超调甚至振荡。解决在控制器输出端加饱和模块典型转速范围由电机模型决定。同时在内环积分项加上限幅常用做法是给积分器一个输出限制或者用条件积分法当控制量饱和时停止积分。Simulink里可以直接用Saturation模块但要注意把它放在混控之前还是之后。我的习惯是先把四个转速指令限幅再把力矩指令限幅两个都要。另外电机响应延迟也很关键。很多模型把电机视为理想执行器输入转速平方立即输出推力这会导致仿真结果过于乐观。实际电机至少有一阶惯性环节时间常数约20~50毫秒。在仿真中加入一阶惯性环节控制参数才更接近实物。你可以用Simulink的Transfer Fcn模块1/(tau*s1)tau先取0.03试。5.4 数据记录与可视化陷阱现象仿真数据看不出问题但导出到MATLAB工作区后曲线异常比如时间轴重复或数据长度不对。有一次我导出的姿态数据只有最后一个时刻的值前几秒全部丢失。原因Simulink的To Workspace模块默认输出的是带时间戳的结构体如果用普通数组接收可能会被压成二维矩阵或者因为保存格式设置问题导致丢数据。常见的是Save format选了Array而信号本身是向量容易混淆维度。解决设置To Workspace的保存格式为Timeseries然后导出后使用get函数提取时间和数据。我习惯用logsout方式记录信号并勾选“Record all signals”这样数据完整性最好。如果仿真里多个信号同名logsout会自动加下标提取时用getElementByName。还有一个可视化陷阱在Scope里看到的多条曲线如果量纲不同会被压缩到同一坐标轴。比如角速度约1 rad/s姿态角约0.1 rad分不清哪条是哪条。建议把信号分开显示或者用plot脚本绘制时分别设置subplot。这个习惯能大幅减少用眼睛判断的误差。6. 进阶从仿真到硬件在环前的三个验证技巧6.1 用蒙特卡洛跑参数鲁棒性调好一组参数后别急着认为完事了。我习惯把质量、惯性矩、推力系数各加减10%然后随机组合跑50次仿真统计姿态稳定时间。如果某组参数下模型发散说明控制器鲁棒性不足。这个验证在答辩时提到会很加分。具体实现可以写个循环每次用随机扰动后的参数调用sim然后把结果保存到数组里。6.2 把S函数封装成独立模块如果需要把模型移植到其他工程甚至对接硬件飞控建议把quad_dynamics改写成Simulink S函数。用level-2的Matlab S函数封装好处是采样时间可控且可以设置过零检测。封装后整个模型就像一个独立组件在硬件在环测试里可以被调用来模拟被控对象。改写的核心是把quad_dynamics的输出作为S函数的Output把连续状态导数放到Derivatives方法里。6.3 检查仿真与实际飞控的时序差异仿真里的控制器更新频率如果和实物飞控不一致参数移植后大概率翻车。我一般会在仿真里把控制器改成离散时间用200Hz或500Hz采样再看性能退化曲线。如果悬停性能在500Hz下良好2kHz下性能没有明显提升说明参数可以接受。反之如果在500Hz下已经明显振荡就必须重新整定参数而不是指望提高采样频率能解决问题。从那以后我每次做四旋翼毕设都强制走一遍“开环验证—闭环调参—鲁棒性测试”这三步哪怕时间再紧也不跳。这套流程帮我少熬了无数个通宵。希望今天这些经验能帮到你。本文还有配套的精品资源点击获取 
长河君 · 2026/10/20 9:49:00 · 阅读 1.2k 四旋翼Matlab仿真与串级PID控制：从动力学建模到毕设落地
简介这份Matlab四旋翼仿真资源面向毕业设计、课程设计及无人机控制方向初学者提供可运行的四旋翼建模、控制算法与仿真验证工具。压缩包共280个文件主要包含125个m脚本用于核心算法与仿真逻辑37个csv数据文件对应飞行位置、速度等实验结果记录13个mat与6个fig文件可用于状态数据读取、曲线查看与结果可视化另有msg、srv等辅助文件便于功能扩展与二次开发。整体约6.63MB结构清晰、轻量易用。目前已有132人学习下载所有源码经严格测试可直接运行适合需要快速搭建仿真环境或完成课程项目的人群。通过本资源可获得完整的四旋翼模型、控制算法实现、实验数据及可视化结果帮助理解飞行动力学与控制原理同时也可以作为毕业设计、课程论文或Matlab编程练习的重要参考素材。1. 毕业设计选四旋翼仿真Matlab这套资源在解决什么做毕业设计选四旋翼方向的人大部分时间其实不是耗在写代码而是耗在“模型跑不起来”和“控制器没法调”这两件事上。Matlab平台的优势在于生态完整Simulink里搭控制回路、改参数、看波形都比纯手写C直观得多。如果手里有一套Matlab四旋翼仿真框架把动力学建模、控制律、Simulink仿真串成闭环工作量能省掉一半以上。这套资源正是围绕毕业设计场景打包的从牛顿-欧拉方程到串级PID再到可视化输出基本覆盖了课设和毕设最常见的技术路线。适合那些刚接触四旋翼建模、又想在答辩前拿到稳定仿真结果的同学。我自己当年就是从类似的框架入手的踩了不少坑才把模型从“飞出去”改成“停得住”。2. 四旋翼动力学建模从刚体运动学到状态空间方程2.1 坐标系定义与欧拉角约定四旋翼仿真第一课不是写方程而是先把坐标系钉死。我见过太多成品代码里俯仰角正负号反了或者欧拉角顺序用了两种结果模型一跑就翻跟头。推荐的做法是地面系定义为NED北-东-地或ENU都可以但机体轴系必须右手定则。这里以多数教材采用的ZYX欧拉角为例先偏航角ψ再俯仰角θ最后横滚角φ。如果把机体坐标系到地面坐标系的旋转矩阵写成一个函数Matlab里通常是function R rotation_matrix(phi, theta, psi) Rz [cos(psi) -sin(psi) 0; sin(psi) cos(psi) 0; 0 0 1]; Ry [cos(theta) 0 sin(theta); 0 1 0; -sin(theta) 0 cos(theta)]; Rx [1 0 0; 0 cos(phi) -sin(phi); 0 sin(phi) cos(phi)]; R Rz * Ry * Rx; end这个矩阵的作用是把机体坐标系下的力、速度转换到地面系。注意Matlab的三角运算默认弧度初始化时千万别把角度直接塞进来。我自己习惯在工程入口写一个转换函数deg2rad统一处理。为了验证坐标系没写错可以做个简单测试输入φ0θ0ψpi/2旋转矩阵应该把机体x轴转到地面y轴方向。如果得到的结果是x轴仍然指向原方向说明旋转矩阵顺序错了。这种小测试能省下后面排查的时间。另外还有一个容易忽略的点欧拉角约定必须与你的控制期望一致。如果你在控制器里把期望俯仰角定义为绕机体y轴旋转而旋转矩阵里用的是ZYX顺序那外环计算出来的期望角速度就会对不上。建议在初始化脚本里用注释写明欧拉角顺序并把这个注释同步复制到控制器模块里。2.2 牛顿-欧拉方程的Matlab实现四旋翼的刚体运动学分为平动和转动两部分。平动方程基于牛顿第二定律转动方程基于欧拉方程。这里忽略桨叶挥舞、气动阻力高阶项保留主要力与力矩每个旋翼产生的推力T_i k * w_i^2k是推力系数力矩包括推力产生的滚转/俯仰力矩、反扭矩产生的偏航力矩把这套方程写成Matlab函数输入状态向量与四个电机转速输出状态导数function x_dot quad_dynamics(x, u, param) % x [x, y, z, vx, vy, vz, phi, theta, psi, p, q, r] % u [w1^2, w2^2, w3^2, w4^2] 电机转速平方 g param.g; m param.m; k param.k; % 推力系数 d param.d; % 反扭矩系数 Ixx param.Ixx; Iyy param.Iyy; Izz param.Izz; l param.l; % 机臂长度 phi x(7); theta x(8); psi x(9); p x(10); q x(11); r x(12); T1 k * u(1); T2 k * u(2); T3 k * u(3); T4 k * u(4); T_total T1 T2 T3 T4; % 平动加速度地面系 R rotation_matrix(phi, theta, psi); acc [0; 0; -g] R * [0; 0; T_total] / m; % 力矩机体系X型布局 Mx l * k * (u(1) - u(3)); My l * k * (u(2) - u(4)); Mz d * (u(1) - u(2) u(3) - u(4)); % 角加速度欧拉方程 p_dot (Mx (Iyy - Izz) * q * r) / Ixx; q_dot (My (Izz - Ixx) * p * r) / Iyy; r_dot (Mz (Ixx - Iyy) * p * q) / Izz; % 欧拉角速率注意 theta -90 度时奇异 phi_dot p sin(theta) * (q * sin(phi) r * cos(phi)) / cos(theta); theta_dot q * cos(phi) - r * sin(phi); psi_dot (q * sin(phi) r * cos(phi)) / cos(theta); x_dot [x(4); x(5); x(6); acc(1); acc(2); acc(3); ... phi_dot; theta_dot; psi_dot; p_dot; q_dot; r_dot]; end这里有几个关键点。首先u是电机转速的平方而不是转速本身。这样写的好处是推力与转速的关系变成了线性系数k后续混控矩阵可以直接用。如果你拿到的数据是转速表的记得先平方再做输入。其次力矩公式里的符号。X型布局常见的电机编号是1号右前逆时针2号左后顺时针3号左前逆时针4号右后顺时针。滚转力矩由x轴两侧推力差产生所以我这里写成l * k * (u(1) - u(3))。如果你的布局是十字型符号需要重新推一遍。这个符号错了模型会表现为“给一个正的期望滚转角飞机却往反方向滚”。第三欧拉角速率公式里有一个1/cos(theta)当俯仰角接近90度时数值会发散。毕设阶段如果只做悬停和轨迹跟踪通常不会碰到大角度但如果你要做翻转动作必须换成四元数表示。这个话题在第五章排查里还会展开。2.3 状态空间模型与接口设计为了和控制器衔接把动力学整理成状态空间形式。我们选状态向量X [x, y, z, vx, vy, vz, phi, theta, psi, p, q, r]控制输入U [u1, u2, u3, u4]。将2.2的导数函数写成x_dot f(x, u)。注意这个f是非线性的不是线性状态空间。在悬停点附近线性化后可以写成矩阵A、B用来做LQR或者稳定性分析。在资源包里通常会同时给出非线性模型和线性化矩阵。线性化方法可以直接用Matlab的linmod或符号工具箱的jacobian。我们在工程里更常用的是在悬停姿态φθψ≈0、速度为零、推力等于重力处求数值雅可比% 在悬停点求雅可比 x0 zeros(12,1); u0 [m*g/4; m*g/4; m*g/4; m*g/4]; A jacobian((x) quad_dynamics(x, u0, param), x0); B jacobian((u) quad_dynamics(x0, u, param), u0);注意这里jacobian需要Symbolic Math Toolbox如果没装可以用数值差分法自己算。这一步不是必需的但做了之后你对模型的固有频率、耦合关系会更有底。比如A矩阵里位置和姿态之间的耦合项是哪些直接就能看出来。接口设计方面我习惯把quad_dynamics这个函数封装成Simulink的MATLAB Function模块输入是状态向量和控制量输出是导数。这样后续做S函数或硬件在环时函数体可以原样复用。注意输入输出都用double类型Simulink里默认是double别改成single否则精度损失会导致漂移更严重。状态空间模型还有另一个用途设计观测器。比如很多毕设题目要求“基于状态观测器的姿态控制”这时你需要在仿真里引入状态重构环节。把非线性模型作为被控对象把线性模型作为观测器的标称模型这是标准做法。资源包里如果给了A、B矩阵你要能说明白它是在哪个平衡点线性化的答辩老师最常问的就是这个。3. PID控制器设计与调参串级结构在Simulink里的落地3.1 串级PID结构内环角速度、外环姿态四旋翼控制最常用的是串级PID内环控制角速度(p,q,r)外环控制姿态角(φ,θ,ψ)。外环输出期望角速度给内环内环输出期望力矩再通过混控矩阵生成四个电机的转速指令。这种结构的好处是内环能快速抑制扰动外环可以慢一点抗干扰能力强。在Simulink里可以把内环和外环分开封装成两个子系统便于调试。资源包里一般会有现成的参考设计但最好自己动手连一遍否则答辩时被问到底层逻辑容易卡壳。我见过不少人拿到的模型只给了封装好的子系统双击进去看到一堆常量根本不知道每个数是什么含义。这种黑匣子是毕设答辩的大忌。串级PID与单环PID的区别在于单环直接对姿态角误差输出力矩相当于把角速度环忽略了。角速度环的存在可以极大提升阻尼特性让姿态响应更平滑。对于四旋翼这种开环不稳定的系统强烈建议使用串级结构。如果你发现资源包里给的是单环PID建议改成双环哪怕只加一个P控制器也会有明显改善。3.2 控制器代码与参数初始化典型的姿态控制器可以写成Matlab函数使用位置式PID。下面给一个内环加外环实现的简化版function tau attitude_controller(phi, theta, psi, p, q, r, phi_d, theta_d, psi_d, param) % 外环姿态角误差 - 期望角速度 kp_phi param.kp_phi; kp_theta param.kp_theta; kp_psi param.kp_psi; p_d kp_phi * (phi_d - phi); q_d kp_theta * (theta_d - theta); r_d kp_psi * (psi_d - psi); % 内环角速度误差 - 力矩 kp_p param.kp_p; kp_q param.kp_q; kp_r param.kp_r; Mx kp_p * (p_d - p); My kp_q * (q_d - q); Mz kp_r * (r_d - r); tau [Mx; My; Mz]; end这里只写了P控制为了展示骨架。实际工程中内外环都要加积分和微分。积分项用来消除稳态误差但四旋翼悬停时积分累积很容易超调甚至振荡所以通常加抗饱和限幅。微分项放在测量噪声上会放大抖动一般用低通滤波后的角速度。参数初始化是一个复杂的活。给一个初始参数表参数典型初始值说明kp_phi/kp_theta5~8外环比例过大则姿态回正过快容易振荡kp_psi3~5偏航响应可以慢一些kp_p/kp_q/kp_r0.1~0.3内环比例单位是力矩比角速度ki0.01~0.05积分系数先设0稳定后再加kd0.01~0.05微分系数先设0这些值不是通用的取决于模型里的惯性矩和机臂长度。我用的是一个某小型四旋翼的参数质量1.1kg臂长0.22mIxx约0.012kg·m²。如果你手里是其他机型参数得重新标定。比如质量变成5kg推力系数也要成倍调整控制器比例系数可能就完全不对了。混控矩阵怎么接姿态控制器输出的tau是力矩需要转换成四个电机的转速平方。常用的是% 混控矩阵力矩 总推力 - 电机平方 mix [1 1 1 1; 1 0 -1 0; 0 -1 0 1; 1 -1 1 -1] / 4; % 实际要根据布局正负号调整这个矩阵的意思是把[T_total, Mx, My, Mz]映射到[u1,u2,u3,u4]。混控矩阵的系数需要与动力学里的力矩公式完全对应否则控制器输出对的电机响应却是错的。一个简单的验证方法把四个转速设为某个只产生滚转力矩的组合跑一下动力学函数看输出的力矩是否与期望一致。3.3 调参顺序与响应曲线判读调参常用的顺序是先内环后外环。把外环踩掉直接给角速度阶跃调内环P让响应快且无振荡再加D抑制超调最后加I消除稳态误差。内环满意后再闭合外环调外环P让姿态角阶跃响应达到期望水平。具体操作时我习惯在Simulink里把参考输入设为一个幅值0.3 rad/s的阶跃信号。注意不要太大角速度环的线性区间有限。观察角速度响应曲线如果出现持续等幅振荡说明P过大如果响应慢吞吞说明P太小或者模型里惯性参数设大了。内环调好后外环姿态阶跃可以给5度的阶跃约0.087 rad。如果外环响应超调量超过20%适当减小外环P。四旋翼的滚转和俯仰通道在参数上可以对称偏航通道单独调因为偏航转动惯量通常更大。这里有一个血泪经验调参时千万别同时改多个参数否则振荡来源根本定位不到。我习惯每次只改一个系数跑一轮仿真记录下超调和调节时间形成一张表再判断。这个习惯救了我很多次。调参时还容易犯一个错误在仿真结果里同时看位置和姿态发现位置漂移就以为是控制器问题。实际上位置漂移很可能是初始姿态角没设为零或者模型的重力方向写错了。我建议先只观察姿态角把姿态角稳定性调好之后再观察位置。如果位置还在缓慢漂移检查动力学函数里平动加速度的计算特别是R * [0;0;T_total]这一项是否确实是重力加推进力的方向。数据判读上我一般记录三个指标上升时间从10%到90%、超调量、调节时间进入±2%误差带。把这些指标做成一张表每次调参后更新。这样做的好处是你能直观看到参数变化的趋势而不是凭感觉瞎试。4. 仿真环境初始化与配置从零搭起一个可复现的四旋翼工程4.1 工程目录与脚本结构一个毕业设计要交的代码不只是模型文件还要别人打开就能跑。我见过很多同学把脚本和模型堆在一个目录里命名叫final_new2.m这种交付方式在评阅时很吃亏。推荐的结构是project/ ├── model/ │ ├── quad_dynamics.m │ ├── rotation_matrix.m │ └── attitude_controller.m ├── sim/ │ ├── main_sim.slx │ └── sim_config.m ├── data/ │ ├── results.mat │ └── logs/ ├── doc/ │ └── 说明文档.md └── init.minit.m是总入口里面设置工作区变量并调用sim_config.m。这套资源里通常会有一个类似的模板但目录名不一定一样。强烈建议你保留自己的README写清楚每个文件的作用和运行方法。我看到很多人的项目文件之间是隐式依赖比如某个脚本里调用了另一个脚本里定义的变量但别人找不到这个脚本。这种情况在评阅时容易被扣分。正确做法是所有参数都在init.m中定义模型和控制器脚本只从param结构体读取。这样任何人拿到工程先运行init.m再打开main_sim.slx就能复现结果。另外init.m最好把addpath加好把model和sim目录都包含进去避免依赖当前工作目录。4.2 初始状态与物理参数配置初始化脚本里最重要的就是物理参数。用param结构体统一存放质量、惯性矩、推力系数等不要散落在各个脚本里。例如param.g 9.81; param.m 1.1; param.l 0.22; param.k 1.6e-5; % 推力系数单位 N/(rad/s)^2 param.d 2.3e-7; % 反扭矩系数 param.Ixx 0.012; param.Iyy 0.012; param.Izz 0.023;这些系数怎么来最靠谱的是从旋翼叶素理论估算但毕设阶段常用的是查同尺寸机型的公开参数再微调。注意推力系数是谐振的单位要统一成N/(rad/s)^2不要和rpm混用。如果模型用转速平方作为输入初始化时就必须保证单位一致。我一个同学把转速单位用成了rpm结果推力差了大约1000倍模型直接变成自由落体。初始状态配置也很关键。通常需要设置初始位置x0、速度、姿态角。我建议把初始姿态角设为零向量方便验证控制器是否能让飞机从零姿态稳定悬停。如果初始姿态角设了非零值调试时你分不清是控制器纠偏有效还是初始条件导致的响应。另外初始高度可以设为1米这样能直观看到高度方向的响应。4.3 数值积分步长与仿真时长设置Simulink里求解器设置是门玄学。固定步长和变步长结果差异很大四旋翼这种强耦合系统我一般用固定步长ode4四阶龙格库塔步长1e-3秒。为什么不用变步长因为变步长在刚度大的情况下会自动缩小步长导致仿真时间不可控而且不利于结果复现。步长太大控制更新会变慢步长太小仿真时间成倍增加。1e-3秒是一个工程上比较稳妥的折中。仿真时长方面悬停仿真通常设10秒就够了你看姿态稳定过程完全够。如果是做轨迹跟踪设60秒。但注意Simulink的Stop Time单位是秒别和采样时间搞混。有时候仿真跑到一半曲线不更新了查看是不是Stop Time设成了0。还有一个容易忽略的点控制器和模型如果都放在连续时间域里仿真没问题但如果用离散PID必须设置采样周期比如200Hz。我的建议是先把控制器也做成连续时间跑通后再改成离散版本。很多成品模型直接给离散控制器同学拿去后不知道采样周期是多少结果控制效果完全不一样。另外Simulink的代数环问题在四旋翼仿真里也容易出现。如果你把控制器输出直接反馈到动力学输入而动力学输出又直接决定控制器输入中间没有延迟就会形成代数环。解决方案是给反馈信号加一个Memory模块或单位延迟打破代数环。我见过不少人在初始化时把Save options设置为timeseries但代数环导致的错误提示却看不懂这里顺便提一句。5. 四旋翼仿真的常见问题与排查失稳、漂移与数据跳变5.1 模型失稳现象、原因与解决现象仿真一开始姿态角快速发散几秒后数值变成NaN控制完全失灵。原因最常见的是推力系数或质量不匹配导致悬停推力小于重力飞机直接坠地或者控制器正负号反了输出的是正反馈。我见过一个案例控制器输出力矩的正负号和动力学里力矩公式的正负号相反结果姿态误差越大控制器给的力矩越让误差变大1秒内就翻了。解决先把控制器拆掉给四个电机一个固定输入使总推力等于重力看模型是否能稳定悬停。如果连这个开环测试都发散说明动力学函数或者坐标系有bug。我通常会在初始化脚本里算一下u0 m*g/(4*k)用这个值去测试开环。如果开环下飞机垂直方向保持不动说明平动模型正确如果姿态角开始翻动说明转动模型或旋转矩阵有误。还有一个常见原因是仿真步长太大。我曾经把固定步长设为0.01秒结果显示角速度在每步之间变化剧烈模型发散。对于强非线性系统步长需要小到能捕获系统动态。如果发散时曲线出现锯齿状优先检查步长。5.2 姿态角漂移现象、原因与解决现象控制器明明让期望姿态角为零但实际姿态角长时间稳定在某个小角度比如俯仰角一直偏着2度。原因可能是欧拉角速率公式中混入了错误项或者机体质心和推力中心不重合建模时没考虑。也可能是传感器噪声模型导致的积分漂移但仿真里一般没有传感器所以重点查刚体方程。解决在quad_dynamics函数里单独检查平动和转动部分。比如把转动部分单独跑固定位置不变输入一个常值力矩看角加速度是否符合预期。另外检查欧拉角速率公式是否有1/cos(theta)项缺失。我在某个版本里就是把phi_dot直接写成了p导致大角度时完全对不上。还有一类漂移来自控制器积分项。如果外环积分系数调得过大悬停时一点点姿态误差会被积分放大形成一个小角度偏移。排查方法是把控制器里的积分项全部置零看姿态是否回到零。如果回到零说明是积分参数问题如果仍然有偏移则回到模型检查。5.3 控制量饱和与执行器限幅现象给定的期望速度很大控制器输出力矩也很大但电机转速被限幅后实际出力不够飞机响应滞后。仿真里表现为姿态角回正速度远慢于预期甚至出现缓慢的低频振荡。原因真实电机转速不可能无限大模型里如果不加限幅仿真数据看起来正常但换到实物就会翻车。控制量饱和后积分项继续累积导致“积分饱和”现象。积分饱和会让控制器在退出饱和后仍然保持一个很大的积分输出引起超调甚至振荡。解决在控制器输出端加饱和模块典型转速范围由电机模型决定。同时在内环积分项加上限幅常用做法是给积分器一个输出限制或者用条件积分法当控制量饱和时停止积分。Simulink里可以直接用Saturation模块但要注意把它放在混控之前还是之后。我的习惯是先把四个转速指令限幅再把力矩指令限幅两个都要。另外电机响应延迟也很关键。很多模型把电机视为理想执行器输入转速平方立即输出推力这会导致仿真结果过于乐观。实际电机至少有一阶惯性环节时间常数约20~50毫秒。在仿真中加入一阶惯性环节控制参数才更接近实物。你可以用Simulink的Transfer Fcn模块1/(tau*s1)tau先取0.03试。5.4 数据记录与可视化陷阱现象仿真数据看不出问题但导出到MATLAB工作区后曲线异常比如时间轴重复或数据长度不对。有一次我导出的姿态数据只有最后一个时刻的值前几秒全部丢失。原因Simulink的To Workspace模块默认输出的是带时间戳的结构体如果用普通数组接收可能会被压成二维矩阵或者因为保存格式设置问题导致丢数据。常见的是Save format选了Array而信号本身是向量容易混淆维度。解决设置To Workspace的保存格式为Timeseries然后导出后使用get函数提取时间和数据。我习惯用logsout方式记录信号并勾选“Record all signals”这样数据完整性最好。如果仿真里多个信号同名logsout会自动加下标提取时用getElementByName。还有一个可视化陷阱在Scope里看到的多条曲线如果量纲不同会被压缩到同一坐标轴。比如角速度约1 rad/s姿态角约0.1 rad分不清哪条是哪条。建议把信号分开显示或者用plot脚本绘制时分别设置subplot。这个习惯能大幅减少用眼睛判断的误差。6. 进阶从仿真到硬件在环前的三个验证技巧6.1 用蒙特卡洛跑参数鲁棒性调好一组参数后别急着认为完事了。我习惯把质量、惯性矩、推力系数各加减10%然后随机组合跑50次仿真统计姿态稳定时间。如果某组参数下模型发散说明控制器鲁棒性不足。这个验证在答辩时提到会很加分。具体实现可以写个循环每次用随机扰动后的参数调用sim然后把结果保存到数组里。6.2 把S函数封装成独立模块如果需要把模型移植到其他工程甚至对接硬件飞控建议把quad_dynamics改写成Simulink S函数。用level-2的Matlab S函数封装好处是采样时间可控且可以设置过零检测。封装后整个模型就像一个独立组件在硬件在环测试里可以被调用来模拟被控对象。改写的核心是把quad_dynamics的输出作为S函数的Output把连续状态导数放到Derivatives方法里。6.3 检查仿真与实际飞控的时序差异仿真里的控制器更新频率如果和实物飞控不一致参数移植后大概率翻车。我一般会在仿真里把控制器改成离散时间用200Hz或500Hz采样再看性能退化曲线。如果悬停性能在500Hz下良好2kHz下性能没有明显提升说明参数可以接受。反之如果在500Hz下已经明显振荡就必须重新整定参数而不是指望提高采样频率能解决问题。从那以后我每次做四旋翼毕设都强制走一遍“开环验证—闭环调参—鲁棒性测试”这三步哪怕时间再紧也不跳。这套流程帮我少熬了无数个通宵。希望今天这些经验能帮到你。本文还有配套的精品资源点击获取