★ FEATURED ARTICLE 去年我替一位做毕设的学弟把关论文复现工作他卡在“把一篇期刊论文里的四旋翼PID控制器搬进Simulink”这件事上整整一周。模型搭出来了仿真一跑就飞走PID参数照着论文抄也没用最后查出来是欧拉角旋转顺序和论文里默认的不一致。这篇文章我就把完整的复现路径重新梳理一遍从选论文、抄方程、搭Simulink模型到串级PID结构、参数整定和结果比对每一步都讲清楚“为什么要这样做”。它在做的事就是把一篇四旋翼PID控制论文从纸面变成可运行的仿真系统。适合读这篇的人有三类正在做四旋翼相关课题的研究生、准备无人机方向毕设的本科生、以及那些拿到论文但苦于“模型搭不出来、曲线对不上”的工程师。整篇内容不挑平台Simulink也好其他仿真工具也好核心思路完全可以平移。1. 复现前的三个认知选论文、选模型、选工具1.1 什么样的论文适合用Simulink做PID复现不是所有论文都值得花时间复现。我见过不少人拿到一篇优化算法改进类的论文就开始搭模型搭到一半才发现——论文根本就没给动力学方程控制器参数也是个黑盒。这种论文只能“参考思路”谈不上“复现”。适合用Simulink做PID复现的论文至少得满足三个条件第三论文里有可对比的响应曲线。阶跃响应、轨迹跟踪曲线、扰动响应曲线都行最好是给了一组明确的仿真工况比如“給0.5 m/s的阶跃速度指令”或者“施加0.1 N·m的常值力矩扰动”。没有这些你复现完也不知道自己到底成了没有。技术类论文的条件其实不难记忆方程完整、参数表齐全、有响应曲线。三条都满足就可以动手缺一条复现难度会成倍增加。1.2 拿到论文后的第一件事不是写代码而是抄方程很多人上来就打开Simulink拖模块这是最大的误区。Simulink模型是“方程的可视化”方程没吃透模型必定出问题。我的习惯是准备一张表格把论文里的关键信息全部摘出来信息类别要抄什么抄的时候注意什么坐标系定义惯性系方向、机体轴方向是NED还是ENU很多论文不写默认NED欧拉角顺序姿态角定义顺序ZYX还是XYZ直接影响仿真曲线状态量位置、速度、角速度、四元数/欧拉角状态向量各下标含义后面接线要用动力学方程力方程、力矩方程有没有空气阻力项、陀螺力矩项执行器模型电机升力系数、反扭矩系数、时间常数单位统一到国际单位制常有坑控制器结构串级还是单环各环节PID形式位置环输出什么姿态环输入什么PID参数表P、I、D数值和量纲有的论文参数是归一化的不能直接抄仿真设置采样时间、求解器、初始状态论文常不写需要按经验补齐表格填完你会发现自己对论文的理解已经上了一个台阶。这个步骤花不了半天但能省下后面一周的调试时间。1.3 为什么选Simulink做这件事在Simulink里复现四旋翼PID控制有几个实打实的优势。第一可视化程度高。控制器结构、反馈路径、饱和模块、混控矩阵全部以模块图的形式呈现。写论文汇报时把模型截图一贴审稿人不用看代码就能理解你的控制结构。第二调试手段丰富。Scope波形实时查看PID参数开个旋钮模块就能在线调整不需要改代码重新跑。这对调参阶段极其友好。第三和实物的衔接性好。控制器调好后可以用自动代码生成工具把控制逻辑导出成C代码直接烧进飞控。很多做实物项目的人选择Simulink就是看中这个“仿真和实物共用一套模型”的链路。但它也有短板。Simulink的模块封装层级一多信号追踪就变得头疼如果直接把动力学方程用几百个Gain、Product模块堆出来出了问题很难定位。我的实践结论是用MATLAB Function模块把动力学方程集中封装控制器部分才用细颗粒度的PID模块搭建。这样既能发挥Simulink的控制器可视化优势又避开了“纯模块搭建”的调试噩梦。2. 从动力学方程到Simulink模块搭建六自由度仿真模型2.1 坐标系与欧拉角约定先统一再动手四旋翼建模的第一步是确定坐标系约定。这个步骤极其枯燥但几乎决定后续所有正负号是否正确。大多数四旋翼控制论文采用NED惯性坐标系加前右上机体坐标系姿态用欧拉角描述。欧拉角的旋转顺序论文里经常不写清楚但仿真结果对顺序很敏感。我的经验是如果论文没有特别注明优先采用ZYX顺序也就是先偏航、再俯仰、最后横滚。这个顺序也是飞行器领域最常用的约定。从机体坐标系到惯性坐标系的旋转矩阵在ZYX约定下写出来是R_bi [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θ]这个矩阵在Simulink里通常是藏在MATLAB Function内部而不是用模块搭因为公式太长用模块搭容易写错下标。还有一个容易被忽略的点状态量的积分参考系。位置方程在惯性系积分速度方程在机体系积分姿态欧拉角由角速度通过坐标转换矩阵积分而来。这三类状态量在不同坐标系下彼此之间的耦合关系不一样接线之前脑子里先过一遍哪里用坐标系转换哪里直接用加减法别搞混了。2.2 六自由度动力学方程的组装方式拿到论文的动力学方程之后我的标准做法是先把方程整理成状态空间形式也就是把所有状态量的导数写成状态量加控制输入的函数然后再翻译成MATLAB Function。四旋翼的刚体动力学通常拆成两部分平动方程在机体坐标系表达dot(u) (F_thrust_x / m) - g * sinθ r*v - w*q dot(v) (F_thrust_y / m) g * sinφ*cosθ w*u - r*u?这里我不再逐项展开因为不同论文符号定义差异很大直接套方程容易出错。核心思路是来自论文的方程符号完全照抄只做单位统一。别自作聪明“优化”人家的方程那是后面的事。转动方程绕质心转动J_x * dot(p) (J_y - J_z) * q*r M_x J_y * dot(q) (J_z - J_x) * p*r M_y J_z * dot(r) (J_x - J_y) * p*q M_z其中J_x、J_y、J_z是三个主轴的转动惯量M_x、M_y、M_z是机体轴力矩。注意陀螺力矩项(J_y - J_z)qr很多简化版的论文会忽略它但复现时如果论文里写了就必须保留否则高速机动时的响应曲线会有可见差异。在Simulink里我推荐把整个六自由度模型写成一段MATLAB Function输入是四个电机的升力和当前状态输出是状态导数。这样从方程到代码的翻译是一一对应的哪里错了也容易对照。以下是模型主体的一段示意代码实测中我最常用的结构是这样的function [state_dot] plant_6dof(state, F, M) % state [u v w p q r phi theta psi x y z] % F为总升力M为三轴力矩 u state(1); v state(2); w state(3); p state(4); q state(5); r state(6); phi state(7); theta state(8); psi state(9); % 转动惯量从论文参数表里读出来 Jx 0.004; Jy 0.004; Jz 0.0084; % 转动方程 p_dot (Jy - Jz)/Jx * q * r M(1)/Jx; q_dot (Jz - Jx)/Jy * p * r M(2)/Jy; r_dot (Jx - Jy)/Jz * p * q M(3)/Jz; % 姿态运动学欧拉角速度与机体角速度的关系ZYX约定 phi_dot p tan(theta)*sin(phi)*q tan(theta)*cos(phi)*r; theta_dot cos(phi)*q - sin(phi)*r; psi_dot sin(phi)/cos(theta)*q cos(phi)/cos(theta)*r; % 平动方程省略重力项和推力项的具体展开 state_dot [...]; end这段代码看起来简单但有两点值得提醒一是姿态运动学方程在θ接近±90°时会奇异如果仿真工况要做大角度机动建议换四元数二是代码里的黏性阻尼项、陀螺力矩项一定要和论文逐一核对论文里有就保留论文里没有就别自己加加了会让响应曲线和论文对不上。2.3 电机与执行机构模型别在细节上栽跟头论文里对电机模型的写法五花八门有的写“电机动态忽略”有的给一阶惯性环节有的甚至用二阶振荡环节。这块差异是对复现结果影响最大的细节之一。最常用的电机模型是一阶惯性环节加二次方升力关系F_i k_f * ω_i² Q_i k_m * ω_i² τ_m * dot(ω_i) ω_cmd - ω_i这里的τ_m是电机时间常数典型值在0.02秒到0.2秒之间具体看论文给什么。如果一篇论文说“电机动态忽略”那意味着仿真里直接把ω_i当作指令值没有延迟如果你复现时自己加了0.1秒的时间常数上升时间会比论文长不少然后就陷入“为什么我的曲线和论文不一样”的困惑。我的建议是宁可让模型比论文更简单也不要擅自加复杂度。复现的优先级永远是“和论文一致”高于“模型更真实”。先把论文的简化模型复现到位再考虑加入更多气动细节。电机转速指令到PWM的换算如果论文没给就用简化方式——把期望转速范围映射到0到1油门量就行这一步不影响控制器核心逻辑只要注意在混控分配后加一个转速上限饱和就好。2.4 Simulink模型的总线划分决定后面调试的幸福感模型结构上我最终采用的布局是这样的第一层信号源期望轨迹 → 控制器子系统 → 电机子系统 → 刚体动力学子系统 → 反馈/Scope控制器子系统和动力学子系统之间用总线Bus信号传递方便后期新增信号而不改动整体连线所有PID参数用一个初始化脚本统一加载或者直接在模型里用常量块引用工作区变量。强烈推荐后者这样调参不用每次改模型文件。关于求解器做四旋翼仿真不要用变步长默认求解器。固定步长ode4步长1e-3就在绝大多数场景下足够。如果论文里不是离散控制器Simulink里的PID模块可以直接用连续时间版本简单直接不要为了“数字控制”而做离散化那是后续做实物才需要考虑的问题。还有一件必做不可的事模型初始化时把状态量设到悬停平衡点。四旋翼悬停时油门指令不是0而是和重力平衡的某个值。很多第一次做的人把初始速度、角度全部设0控制器输出也从0开始结果一仿真飞行器直接坠落或者疯狂翻转他们还以为PID调错了。实际上就是初始状态没配平。3. PID控制器的分层设计与串级结构搭建3.1 为什么四旋翼必须用串级而不是单环PID四旋翼控制领域有个常识没有人用单环PID直接控制位置。原因要从物理因果链说起。四旋翼的拉力永远沿着机体轴方向想让无人机水平移动必须先倾斜机体让拉力在水平方向产生分量。也就是说位置控制的执行链条是期望位置偏差 → 期望横滚/俯仰角 → 期望角速度 → 期望力矩 → 电机转速。这条链路里有好几层积分关系如果只用一个位置环PID直接输出电机指令相位延迟会非常大稍有一点扰动系统就振荡。打个比方单环PID就像你开车时只看终点位置、直接踩油门打方向盘完全不管车身当前的姿态和转向速度路稍微弯一点就失控。串级PID则像是分了两层外层导航只给出“该往哪个方向并多少度”内层负责把车身姿态精确控制在那个转角上。内层快外层稳整体才可控。标准的四旋翼串级结构是三层控制叠起来但从实现上看位置环不是必须和姿态环放在同一采样频率。工程上常见的是两层嵌套结构是这样的位置控制器外环p_ref → 期望横滚角φ_cmd、期望俯仰角θ_cmd 姿态控制器中环φ_cmd/θ_cmd/ψ_cmd → 期望角速度p_cmd/q_cmd/r_cmd 角速度控制器内环p_cmd → 期望力矩M_x等等在Simulink里这三个控制器各用一个子系统封装中间的连接线就是期望角信号。检查信号流向比检查代码容易这也是我坚持控制器用模块图搭建的原因。3.2 内环和外环参数之间的“互锁关系”串级PID的整定有个底层逻辑内环的闭环响应速度必须显著快于外环。工程上这个倍数一般取5到10倍。为什么必须这样因为外环是通过内环间接控制被控量的如果内环太慢外环看到的被控量变化总是“跟不上指令”它就会不断增大输出最后导致振荡。反过来如果内环太快会放大测量噪声系统能耗无谓增加。这个互锁关系可以直接换算成带宽。假设角速度内环带宽是ω_inner姿态外环带宽是ω_outer那么ω_inner ≈ (5~10) * ω_outer带宽和PID增益之间存在简单的近似关系。比如角速度环控制器输出是力矩M被控对象是转动惯量J那么M Kp_omega * (ω_cmd - ω) dot(ω) M / J (Kp_omega / J) * (ω_cmd - ω)这个闭环就是一阶惯性环节带宽等于Kp_omega / J。所以角速度环的P增益的估算起点是Kp_omega ≈ J * ω_inner姿态外环的情况也类似。如果控制器把期望姿态角映射成期望角速度ω_cmd Kp_phi * (φ_cmd - φ)且内环足够快那么姿态环闭环带宽近似等于Kp_phi。这一层逻辑理清后参数整定就不再是盲目试凑而是有物理依据的。3.3 控制分配期望推力和力矩如何换算成四个电机转速控制器算出期望的总推力U1和三个轴的力矩U2、U3、U4之后还要把它们分配到四个电机上。这个步骤叫控制分配Control Allocation四旋翼里就是用一个4×4的混控矩阵。以X型机架为例假设四个电机位置在机体对角线方向力和力矩的分配关系为U1 F1 F2 F3 F4 U2 (F1 - F2 - F3 F4) * L_arm / sqrt(2) U3 (F1 F2 - F3 - F4) * L_arm / sqrt(2) U4 (Q1 - Q2 Q3 - Q4)其中F1到F4是电机升力Qi是反扭矩L_arm是电机到质心的力臂。正负号取决于电机旋转方向和机架类型这一步是复现中正负号最容易出错的地方。混控矩阵符号弄反最典型的现象是推油门正常但一给滚转角指令飞行器反向翻转看起来就像模型发散了。在Simulink里混控矩阵可以做成一个Gain块矩阵常量放在工作区。代码一个坏处是调试时看不见中间量D就在混控矩阵的输出端口接几个Scope测试时看四个电机的转速指令是否合理。3.4 限幅、积分抗饱和和微分滤波这些“小部件”决定成败很多复现的论文里画控制框图时只画PID的简洁形式但实际工程中PID旁边几乎都跟着一串辅助模块。缺少它们仿真里的小毛病就会被放大成“模型错了”的误判。限幅位置环输出的期望角要限制在比如±30°以内否则大偏差时飞行器会翻过去角速度环输出的期望力矩也要限幅对应电机的转速上限。没有限幅的话控制器计算出巨大指令电机会立刻饱和系统非线性表现完全跑偏。积分抗饱和这是最容易忽略的一块。电机转速有上限当误差持续偏大时积分项会不断累积导致控制器输出长期卡在上限、误差已经反向了积分还没退出来。这就是典型的积分饱和windup。最简单的做法是积分分离——误差超过阈值就冻结积分更平滑的做法是条件积分只在输出未饱和且误差较小时才积分。微分滤波先说明一下四旋翼内环的角速度信号如果来自数值差分噪声本身就很大再对误差做微分会把噪声进一步放大。因此D项通常不用在误差上而是只对测量值如角速度做负反馈相当于阻尼项。如果需要真正的微分项一定要在它后面串一个低通滤波器截止频率取采样频率的四分之一到八分之一之间。这些细节看起来不起眼但我在实际对比中发现有没有积分抗饱和直接决定大角度机动后悬停恢复是否会出现两三秒的振荡。论文里就算没画你也得补上否则很难复现出稳定的响应曲线。4. 参数整定路径先内环后外环的实用调参法4.1 论文给的参数不能直接抄原因在这里很多复现者一开始都会把论文的参数表原封不动填进模型然后发现曲线对不上——有的超调偏大有的干脆发散。这不一定是你的模型错了更常见的原因是论文参数的“隐含环境”和你的模型不一致。同一个Kp值论文用的电机时间常数是0.05秒你用了0.1秒带宽就明显下降论文的控制器采样时间是1kHz你在连续仿真里用的是变步长ode45两边的离散化延迟完全不同。这些差异累积起来参数自然不能直接通用。正确的姿势是把这组参数当作“参考量级”。比如论文说角速度环Kp0.15你根据上一节带宽公式估算出来的起点如果是0.12那说明两者在一个量级内这就说明你的模型没有大的结构错误。以此为基础再做微调而不是死磕“为什么论文的参数跑不出论文的曲线”。4.2 从力学关系估算PID初值附一个可复用的算例用公式估算初值比我见过很多人在调参工具里随机乱试要高效得多。我把这套估算方法整理成可复用的步骤。第一步确定转动惯量。从论文参数表抄出Jx、Jy、Jz的值。假设Jx Jy 0.01 kg·m²。第二步确定外环带宽。姿态环带宽取1~2 rad/s想快一点就取2。那么姿态外环增益Kp_phi ≈ 2。第三步确定内环带宽。内环要快5~10倍取ω_inner 15 rad/s。那么角速度环增益Kp_omega J * ω_inner 0.01 * 15 0.15。第四步估算积分项。积分项的作用是消除静差一般取Kp的十分之一到五分之一。先取Ki_omega Kp_omega / 10 0.015。如果后面发现稳态误差始终偏大再按比例往上加。第五步微分项先置零。四旋翼仿真里D项通常不是必须的先把P和I调好再看响应是否有振荡。这套算例的目的不是告诉你“参数就是0.15和0.015”而是让你明白估算路径从期望带宽倒推增益比随机试凑有据可依得多。实测下来这个初值通常在合理范围的50%到200%之内剩下的交给阶跃响应微调。4.3 阶跃响应下的逐环调参流程参数整定必须遵循“由内向外”的顺序先调角速度环再闭合姿态环最后闭合位置环。每一步都只能调当前这一层的参数否则多层耦合在一起你根本不知道哪个环节在振荡。具体流程我是这样走的断开姿态外环直接给角速度内环一个方波或阶跃参考信号比如p_cmd从0到0.5 rad/s。调Kp_omega直到角速度响应快速跟踪且无明显振荡。判断标准没有持续振荡超调可接受。如果角速度静差明显逐步加Ki_omega。注意观察超调积分项加多了会让阶跃响应超调变大此时可以稍微回调Kp或者引入一点角速度负反馈微分项。闭合姿态外环给定10°到15°的小角度阶跃比如θ_cmd从0到15°。调Kp_phi观察角度响应。姿态角超调一般要求小于15%如果因为内环跟踪不够快而导致超调过大优先提高内环带宽而不是继续加外环Kp。姿态环稳定后闭合位置环。先只调水平位置给定小位移阶跃比如向前0.5米。位置环的Kp先给一个保守值逐步加大观察位置轨迹有没有超调。最后做整体联动测试看给定一个路径点时水平和垂直两个方向的响应是否协调。每一步调完后用Scope里保存的波形记录当前参数和响应指标。调参全过程一定要做记录别偷懒。我在实际项目里因为偷懒不记录参数“调好的参数一改就再也回不去”的情况出现过不止一次。4.4 我实测的几组参数变化规律打包送给你为了给后来的人一个直观参考我把自己调参过程中观察到的规律整理成了下表。这些规律在大多数简化模型和中等保真度模型上都适用参数偏小的现象偏大的现象调整方向角速度环Kp响应慢指令跟踪滞后高频振荡甚至diverging加Kp到出现轻微振荡后回退20%~30%角速度环Ki常值扰动下角速度静差明显超调增大低频振荡从小往大加配合观测阶跃超调姿态环Kp姿态角跟踪缓慢转弯滞后角度超调伴随轻微抖动保持内环带宽足够前的上限位置环Kp位置收敛慢像“飘过去”位置来回震荡油门变化剧烈加到轻微震荡再回退或加角饱和微分项噪声放大力矩毛刺系统迟滞明显响应变钝优先用于角速度测量值负反馈慎用误差微分最大的体会是四旋翼串级PID里内环带宽是总瓶颈。内环带宽打不上去外环参数怎么调都别扭。遇到“外环一加参数就振荡”的情况先回头检查内环响应够不够快而不是死磕外环参数。5. 复现结果分析与差异排查为什么我的曲线和论文不一样5.1 复现成功的量化标准不看图形看指标复现不是说“看起来差不多”就算成功。我见过太多人对着两条曲线说“差不多”实际上上升时间差了40%。复现成功的判断要用数值指标说话。至少要对以下几个指标做对比上升时间响应从10%到90%稳态值的时间超调量峰值超出稳态值的百分比调节时间进入稳态值±5%或±2%误差带后不再超出的时间稳态误差最终值和目标值的偏差轨迹跟踪误差如果论文给了路径跟踪曲线算RMSE我把和论文的对比结果整理成一张三列表格左边是论文数值中间是我的仿真数值右边是偏差百分比。指标论文值我的复现值偏差阶跃超调量12%13.5%1.5%上升时间0.48s0.52s8%稳态误差0.1°0.08°-20%偏差在合理范围内的比如20%以内且曲线的整体形状一致基本可以判断复现成功。如果偏差超过50%就需要进入排查流程了。5.2 六个最容易导致结果偏离的模型细节我在复现和帮别人排查的过程中发现大部分“曲线对不上”的问题都出在下面这几个地方按命中频率排序欧拉角旋转顺序不一致。论文用ZYX你却在Simulink里用了XYZ或者ZXZ姿态响应在中等角度以上就会出现明显差异而且这种差异很难通过调PID参数消除。电机动态模型的差异。论文忽略电机动态你为了“真实”加了0.1秒时间常数上升时间立刻变慢。这是我见过最多的偏差来源。混控矩阵的正负号或机架类型错配。X型机架和型机架的分配矩阵不同旋向布局不同分配出来的转矩方向完全相反。初始条件没有配平。没有给悬停油门初值仿真从零油门开始控制器要花很长时间把机体拉起来前一两秒的曲线肯定和论文对不上。求解器或步长不合适。变步长求解器在非线性强的时候会自动缩小步长看不出来问题但如果你用固定步长且步长过大比如0.01秒电机时间常数又是0.02秒那么电机动态已经是欠采样状态响应曲线的形态都会变异。重力补偿缺失。有些论文的控制器里带重力补偿前馈项有些放在控制分配里。复现时如果把前馈弄丢了悬停时节流阀不在预期位置小范围扰动倒是还能稳住大机动时就明显不如论文平滑。5.3 用单一变量原则一步步锁定偏差根因排查“曲线对不上”问题最忌讳的是同时改好几个地方。我推荐的排查方法是单一变量原则一次只改一处改完重新仿真对比论文曲线。具体流程可以这样走第一步先做“看不见信号的排查”。核对模型架构、坐标系定义、欧拉角顺序、混控矩阵正负号。这些是一票否决项错一处整个曲线都会乱。第二步对比悬停状态。给定配平初始状态闭合控制器看悬停时的油门输出、姿态角能不能稳住。如果这步都稳不住先解决模型初始化再看别的。第三步给一个小角度的姿态阶跃比如10°对比论文的阶跃响应。用游标卡尺量上升时间、超调量、调节时间记录下来。每改一个参数记录一次不要把数据存在脑子里。第四步如果上升时间偏大优先怀疑电机时间常数把它减小向论文靠近如果超调偏大优先怀疑内环带宽不足或积分项过强如果稳态误差不为零优先检查积分项和重力补偿。大多数偏离问题走完这四步就能找到答案。5.4 分层复现的思路先简化再保真我在几次复现中养成了一个习惯不要一上来追求高保真模型建立两个模型的思路会更高效。第一个是“设计验证模型”忽略空气阻力、忽略电机动态、忽略陀螺力矩项用连续PID控制器直接在Simulink里跑控制算法验证。这个模型简单调试时一眼就能看透控制系统行为适合处理控制器结构和参数整定问题。第二个是“全保真模型”加上电机一阶动态、线性空气阻力、可能存在的传感器延迟和带限噪声让曲线更贴近论文里的“完整仿真结果”。这个模型用来做最终的对比验证。两个模型共用同一套控制器参数初始化脚本切换只需改一个配置常量。这样做的好处是万一全保真模型跑出来的曲线和论文对不上你还能退回简化模型确认“控制系统本身有没有问题”。如果简化模型是好的问题多半出在外围模型参数上如果简化模型都发散那就回去检查控制器结构和运动学方程。这套“先简单后复杂”的思路和论文复现本身无关但它是保证复现工作不至于陷入泥潭的重要方法。我见过太多人一开始就建全保真模型然后花大量时间调阻力系数、调传感器模块参数最后才发现纯粹是控制结构错了。本末倒置是最可惜的时间浪费。
阅读完成 · 觉得有帮助？
