☰
多旋翼PID论文复现:Simulink建模与仿真参数整定全流程
2026/10/20 18:40:19 网站建设 项目流程

“复现论文”这件事,做过的人都知道,最耗时间的往往不是看不懂公式,而是公式里的每一项到底对应Simulink里的哪个模块、哪根信号线。特别是多旋翼这种强耦合、非线性、多输入多输出的对象,光是把论文里的动力学方程“翻译”成仿真模型这一步,就能劝退一拨人。这篇内容我打算围绕“PID期刊论文复现”这个主题,把多旋翼飞行器从动力学建模到控制系统设计,再到Simulink仿真落地的完整链路拆开来讲。重点不是复述教材,而是讲清楚那些论文里不会写、复现时必须自己补上的细节——比如电机响应延迟怎么建模、姿态控制内环外环的带宽怎么拉开、论文给的参数和实际仿真对不上时怎么处理。这篇文章适合正在做相关课题、被“论文复现”折磨过的研究生和工程师,也适合想从零入手多旋翼控制仿真的入门者,我会尽量用实际操作过的经验来填平论文和仿真之间的那道沟。

1. 读论文复现之前:先把坐标系约定和电机细节钉死

1.1 坐标系的选择是一切成败的地基

复现任何一篇多旋翼论文,我建议你上来先干一件事:把论文里所有的坐标系定义抄出来,反复确认。这不是小题大做,90%的复现失败案例,根因不在控制器参数,而在坐标系没对齐。

多旋翼领域常见的做法是“机体坐标系 + 地面惯性系”双坐标系。机体坐标系(通常记作B系)原点在飞行器质心,x轴指向机头方向,y轴指向右翼,z轴向下,符合右手定则。地面坐标系(通常记作E系)原点在地面某固定点,x轴指向正北或论文定义的任意水平方向,z轴垂直向下,这也是很多航空类论文的惯例。

这里第一个坑就来了:同一套姿态角,在不同论文里的旋转顺序可能是不同的。有的论文用ZYX欧拉角(先偏航再俯仰再滚转),有的用ZXY,有的干脆用四元数。如果你复现的论文用的是ZYX,但你的Simulink模型里Rotation Matrix模块默认用的是ZXY,姿态角误差会在几分之一秒内迅速放大,最后整个仿真直接发散。

我的习惯是:在建模前先画一张图,把两个坐标系和三个旋转角(滚转角φ、俯仰角θ、偏航角ψ)的旋转关系画清楚。然后用Rotation Matrix模块和Rotation Angles模块反复验证转换矩阵是否为单位阵(转一圈再转回来,结果应该等于输入)。这一步花10分钟,能给你后面省下三天查bug时间。

1.2 电机模型不是简单的惯性环节

复现论文时长出现这样的问题:论文里的控制器参数明明已经给全了,但你的仿真就是稳不住,姿态角来回震荡。这时候别急着怀疑控制器,先看看电机模型是不是偷工减料了。

很多入门者做电机建模,直接用一组一阶惯性环节把四个电机的转速响应糊弄过去,时间常数往0.01秒一放,看起来有延迟,实则完全不满足控制精度需求。真实的多旋翼电机响应里,螺旋桨的转动惯量、电调的信号延迟和换向延迟、转速与升力之间的非线性关系,都会显著影响控制带宽。

如果论文没有明确给出电机模型,我建议用如下二阶模型近似:

  • 一阶惯性环节传递函数:$\omega_{out}(s)/\omega_{cmd}(s) = 1/(\tau s + 1)$,其中$\tau$取0.02到0.05秒之间(具体取值取决于浆叶尺寸和电机KV值)
  • 在Simulink里实现时,要加Saturation模块限制转速上下限(油门0%~100%对应转速范围)
  • 还要加Rate Limiter模块限制转速变化率,否则仿真启动瞬间会产生物理上不可能的巨大加速度

如果你复现的论文对“电机动态”完全没提,但你发现你的仿真在相同PID参数下比论文更容易发散,大概率问题就在这里。我踩过这个坑之后养成了一个习惯:任何多旋翼仿真,电机延迟和限幅必须建模,哪怕论文省略了,我也会自己补上,并在对比结果时注明这一点。

2. 动力学方程拆成Simulink模块的两种思路

2.1 从方程到模块:刚体动力学建模的三种拆分法

拿到论文的动力学方程之后,你需要把它变成Simulink模型。多旋翼刚体动力学通常包括平动方程和转动方程两大块。平动方程描述的是线加速度,转动方程描述的是角加速度,两者通过姿态角耦合在一起。

第一种拆分思路:直接法

直接把加速度积分得到速度,再积分得到位置。力(推力、重力、空气阻力)之和除以质量得到平动加速度。这样的做法符合直觉,也好调试,但存在一个隐蔽的问题——多旋翼的转动惯量随姿态变化而变化,如果你直接按定常转动惯量去仿真,在某些极端姿态下误差会积累。

第二种拆分思路:矩阵法

把转动方程写成矩阵形式,用四元数或旋转矩阵描述姿态。这样做的好处是姿态角奇异性小,适合大幅机动仿真,但矩阵运算模块在Simulink里不够直观,调试信号时要从矩阵里挑元素,麻烦不少。

第三种思路(我推荐的):代数约束法

直接按“角动量定理 + 欧拉方程”来建模。具体的实现方式是:转动惯量矩阵在机体坐标系下表示(对角线元素为Ixx、Iyy、Izz,非对角线元素通常忽略),角速度p、q、r在机体坐标系下描述,然后利用欧拉方程:

$I \cdot \dot{\omega} + \omega \times (I \cdot \omega) = M$

这里的M是合外力矩。把矩阵乘法用Simulink的Matrix Multiply模块实现,陀螺力矩项用Cross Product模块实现,整个结构非常清晰,而且物理意义明确。仿真的时候注意一点:欧拉方程在惯性系和机体系下的形式不同,必须搞清你用的角速度是在哪个坐标系下描述的。误区在于很多人直接用机体系的角速度代入惯性系方程,结果陀螺力矩符号全反了。

2.2 积分器初始状态:和论文算例保持一致

很多论文都会给出初始条件,比如“初始高度2米,初始姿态角全部为0”。但这些信息经常藏在图注或者仿真设置里,很容易被忽略。如果你的初始姿态角和论文不一致,稳态对比自然对不上。

另一个容易忽略的是trim点(配平点)计算。很多论文直接假设四旋翼在悬停状态初始,即油门刚好等于重力补偿,其余控制输入为0。但在Simulink里,如果你给电机的初始转速是0,或者给的初始油门不足以悬停,飞行器会立刻下坠,姿态控制器拼命修正,仿真初期的数据就是一段“垃圾段”。我的做法是先做一次配平计算,用solve脚本求出悬停油门对应的四个电机转速,作为模型的初始状态。

2.3 力的方向和加速度方向千万别弄反

这个坑我印象特别深。有人跑出来的悬停仿真高度不断漂移,查了一晚上,最后发现是推力方向定义反了。机体系z轴向下,螺旋桨产生的推力在机体系中也沿z轴向下(注意这里的“向下”是相对机体而言,实际上它是把机体往上推的力)。但在地面坐标系中,这个力的方向取决于当前姿态,要经过旋转矩阵变换。

如果你在模型的Force模块里直接把推力写成正数,在机体系z轴上向上,那就错了。正确的是:先算出每个电机在机体系下的推力(沿z轴负方向),然后通过旋转矩阵转换到地面系,再叠加重力。这个“负号”问题,是复现过程中最频繁出现、也最磨人的地方。我的建议是在模型里给力的计算单独建一个子系统,把所有坐标转换都集中放在那一块,方便加Scope看方向对不对。

3. PID控制器三层架构里的参数整定,到底在调什么

3.1 串级PID:内环外环的时间尺度分离

多旋翼的PID控制器通常不是简单的一层。角度环和角速度环构成内环,位置环构成外环。内环速度快、高频响应好,外环速度慢、对噪声要求高。串级结构的内外环时间尺度要拉开,一般内环响应速度是外环的5到10倍。

这里有个工程经验:先整定内环,再整定外环。内环是角速度环,输入角速度误差,输出力矩指令;外环是角度环,输入期望角度,输出角速度期望给内环。如果你一上来同时调两组PID,出了问题根本没法定位。我的做法是先用“手动输入状态量阶跃”的方式把内环调稳。具体做法是在仿真输入处构造一个角速度阶跃信号,给定1弧度每秒的目标角速度,观察当前角速度能否快速跟踪。稳定后再切到角度环,给定固定角度目标,观察角度跟踪效果和超调量。这一步急不得,我一般会花掉整个调试过程一半的时间。

3.2 位置环的“速度环”到底该放在哪

很多论文的三层串级结构是位置环(外环)→速度环(中环)→姿态/角速度环(内环),但也有些论文把速度环省了,直接位置误差经过一个比例环节就生成姿态指令。两种结构各有道理。

如果论文里明确使用了“位置误差→期望加速度→期望姿态角”这种方案,那就不需要单独的速度环,直接用期望加速度解算期望姿态角就行。具体解算过程是这样的:

期望加速度 = 比例项 × 位置误差 − 微分项 × 速度误差 + 前馈加速度

然后通过一个简化模型,用期望加速度的x、y分量去反解期望滚转角和俯仰角:

$\phi_{des} = (a_{des,y} \cdot \sin\psi_{des} - a_{des,x} \cdot \cos\psi_{des}) / g$

$\theta_{des} = (a_{des,x} \cdot \sin\psi_{des} + a_{des,y} \cdot \cos\psi_{des}) / g$

注意这里的符号和坐标系必须和动力学模型完全匹配,否则内外环的“因果关系”是反的,系统会把姿态指令和位置反馈混为一谈。复现论文时,我最常做的事就是把论文公式里的符号一个个和Simulink里的“期望加速度→期望姿态”模块对齐,宁可多花时间也不能猜。

3.3 高度环为什么不能用纯PID

高度通道在物理上是个积分加二阶振动系统,如果只用纯PID,风干扰和电机饱和会导致高度缓慢漂移。所以我倾向于在高度环中加入积分项,同时对积分项做限幅,防止积分饱和。

具体参数上,高度环的Kp对应的是“高度误差产生垂直加速度”的比例,它决定了系统相对于重力补偿的刚度。高度环的Kd对应“垂直速度抑制”项,它和并把垂直速度反馈引入,形成等效的垂直阻尼。这两个参数和目标高度阶跃跟踪的实验结果关系非常直接,调起来比较直观。

4. 论文里没写的数据,怎么靠物理常识和反向推算补齐

4.1 气动阻尼系数:论文最常“省略”的项

大部分CS或控制方向的多旋翼论文,动力学方程都有这么一项:$\dot{v} = g + (T/m)R \cdot e_3 - (B_1/m)v$。这里的$B_1$是平动阻尼系数。但你去翻论文的仿真部分,很大概率看不到$B_1$的具体值和来源。因为作者可能没做真机实验,这个值就是“当前设计的”或“经验假设的”。

复现时你不可能因为论文没给就不建模,不然后面控制器的鲁棒性测试就是空谈。我的做法是:先取一个数量级合理的初值,比如$B_1=0.1$(对应低速飞行时的线性气动阻尼)。然后在模型里做一次“无控制直升机坠落”测试,观察速度响应是否合理。如果仿真出的水平速度衰减过快,说明B1偏大:如果速度变化平滑且衰减缓慢,说明量级还行。再用表格对比不同B1下的阶跃响应差异,做一次简单的敏感性分析。

4.2 旋翼的转动惯量怎么估算

刚体转动惯量可以直接从论文给定的质量、轴距等参数推出来。但旋翼自身的转动惯量往往被忽略。它影响的是电机响应的延迟特性——转动惯量越大,转速变化越慢,系统稳定裕度越小。

如果论文没有给旋翼转动惯量,我会用这样的近似方法:把旋翼看成一根宽度极窄的细杆(长度等于旋翼直径,质量等于旋翼质量),则绕中心垂直轴的转动惯量近似为$J_r = (1/12) m_{rotor} d^2$。从某型号常见的质量为12g、直径0.25m的旋翼计算,$J_r$大约在$6.25\times10^{-5} kg·m^2$这个数量级。把旋翼惯量和电机时间常数结合起来,就能得到一个完整的电机+旋翼动态模型。

4.3 推力系数和扭矩系数:用悬停测试反推

论文中常见的电机模型是$T = k_T \cdot \omega^2$,$M = k_M \cdot \omega^2$。其中推力系数$k_T$和扭矩系数$k_M$来自桨叶数据或风洞实验。如果没有,就需要你自己估。

一个比较靠谱的做法是“悬停反推法”:假设四旋翼悬停时总升力等于重力,单个电机的升力就是$mg/4$,此时把4个电机的转速设成已知悬停转速$\omega_h$,就能反推$k_T = mg/(4\omega_h^2)$。扭矩系数也用类似思路:悬停时每个电机输出扭矩和反扭矩,数值上等于$\tau_m = k_M \omega_h^2$。估完这些后,用Matlab的反向推算脚本(比如利用Simulink的Parameter Estimation模块或手动优化)把你的模型输出和论文的阶跃响应曲线做拟合,微调参数。

这个“先估算再拟合”的流程,是我觉得整个复现过程中最“科研”的部分——不看答案,通过物理规律和系统响应反推未知参数,也是以后做真机测试时非常有用的能力。

5. 从“看着像”到“稳得住”:仿真环境搭建与测试建议

5.1 第一步验证:3自由度静浮台测试

冒然把全6自由度模型跑起来,姿态和位置全耦合,出了bug根本分不清是控制器还是动力学的问题。我建议把搭建好的模型拆成两步验证。

第一步:把平动方程和转动方程解耦。固定飞行器位置,只让姿态变化,模拟“3自由度静浮台”实验。具体实现方法是在Simulink模型里把平动通道的积分反馈断开,或者手动态锁定位置状态变量。这样你就相当于在一个固定点上做姿态控制实验,可以专心调整姿态环参数。第二步:把位置环解开,加入平动反馈,做全6自由度仿真。

我在实际工作中还喜欢加一个更严格的测试:在仿真中设定一个“初始姿态偏差”,比如把初始滚转角设成30度,观察整个系统能否快速恢复水平姿态。如果姿态控制器的带宽或阻尼不足,这个测试会非常灵敏地暴露问题。调参的时候多跑几组这样的测试,比盯着悬停曲线反复调要有用得多。

5.2 传感器噪声和故障注入:论文不写,但你要会设计

学术论文的仿真相。

当于在理想环境里,传感器几乎不加噪声,控制器的鲁棒性测试也很少做。但你做复现,不只是为了和论文曲线完全一样,更是为了把系统验证到一个“看得见使用价值”的程度。所以我强烈建议,在做完论文算例复现之后,额外加入以下测试:

  • 传感器噪声:给姿态角测量加白噪声(幅值在0.01~0.05弧度之间),给位置测量加0.1~0.5米的高斯噪声,观察控制器还能不能保持稳定。如果姿态噪声直接导致了大幅度高频运动,说明你的内环滤波不足或微分项增益过大。
  • 电机失效模拟:随机断开一个电机的输出,观察系统姿态是否有恢复能力。大部分论文根本不做这个测试,但这是工程上多旋翼是否能继续安全降落的底线能力。
  • 风扰仿真:在平动方程里加一个沿x轴的恒定风力(比如5m/s)和随机阵风分量,观察控制器能否把位置误差收敛到可接受范围。

这样测试下来,你才会真正理解论文里每一个控制器参数在“发生了什么”时到底起什么作用,而不是只会在Simulink里拖模块、调增益。

5.3 混控矩阵总距饱和问题

四个电机转速都不可能无限大,转速下限是零,上限是电调和电池能给的最大转速。在风扰或大机动时,如果某个电机一直顶在饱和上限,这个电机的控制力矩就没了,姿态环会试图加大另一个电机的转速来补偿,结果整个系统进入非线性区域,控制器性能断崖式下降。

针对这个问题,我在混控矩阵做了一次归一化处理。混控矩阵把期望的总升力和三轴力矩分配到4个电机的转速指令上,如果某一个电机的期望转速超出上下限,就等比缩放到限幅边界,同时其它三个电机也按同一比例缩放,保证总升力方向不变。这种做法称为“归一化饱和处理”,算是个不大不小的技巧,但缺了它,大机动仿真很容易莫名发散。

6. 模型一致性校验与参数复现的通用方法论

6.1 和论文曲线对不齐时,先查这三件事

我见过太多人复现论文时,第一反应是“我是不是PID参数抄错了”,然后花一整晚调参。其实参数对齐只是最后一步,更基本的前置检查项是这三件:

  • 单位统一:论文里的力是用牛顿还是千克力?角速度是度每秒还是弧度每秒?转速单位是rad/s还是RPM?这个步骤一旦错,后面全错。
  • 采样时间一致性:论文没说控制器和动力学的采样时间是多少,你可能拿默认的连续求解器直接跑。但真实环境下飞控以500Hz到1000Hz离散采样,控制器是离散实现的。如果你在Simulink里用了连续PID模块,却在离散仿真步长下跑,结果可能完全不一样。建议显式设定控制器采样时间(比如0.002秒)和动力学采样时间(用变步长求解器)。
  • 控制器输出限幅:你的PID控制器输出是否加了限幅?论文里多半没写,如果电机转速指令被限幅,而你的模型是无限幅的,初期小扰动差异不大,但大机动下差异立刻拉开。

6.2 模块化设计:同一套代码复现多篇论文

复现模式基于一套“可替换”的Simulink架构。我的做法是把模型拆成四个独立子模块:

  1. 输入生成模块(目标位置、目标姿态信号)
  2. 控制器模块(完全独立,PID参数留成工作区变量)
  3. 混控/执行器模块(电机+电调+限幅)
  4. 刚体动力学模块(包含气动阻尼、重力、推力)

每一篇要复现的论文,只需要改控制器子模块里的结构和参数,其它模块保持公用。这种模块化设计在复现不同论文时特别省力,也方便做横向对比实验。如果一个动力学模型能在同一控制器下跑出多篇论文的姿态和位置响应,说明你的动力学模型本身可靠性很高。

6.3 最后再提两个提高复现质量的工具技巧

第一,Simulink的线性化分析工具很好用。在模型达到悬停状态后,用Model Linearizer在悬停配平点做线性化,得到状态矩阵,再用eig函数看特征值分布。理论上,悬停时闭环系统特征值不应有正实部;如果正实部存在,控制器参数一定有隐患。这个方法比单纯看仿真曲线更快发现系统不稳定根。

第二个技巧是用脚本批量扫参。把PID参数定义为MATLAB工作区变量,用parfor循环批量跑仿真,每组参数记录上升时间、超调量、稳态误差,汇总成一张表。我用这个方法在半个小时内扫了上百组参数,直接找出最接近论文结果的组合。手动在Simulink里调参,很难做到这种效率。

复现论文,最后复现的其实不是公式,而是公式背后的一整套完整逻辑——坐标系怎么定、动态怎么建模、限幅怎么加、控制器和执行的带宽怎么匹配。把这个逻辑真正跑通的仿真模型,才有可能成为你之后做研究、做实验、甚至做进一步算法改进的可靠基础。希望这篇内容能帮你少踩几个我曾经历过的坑。

需要专业的网站建设服务?

联系我们获取免费的网站建设咨询和方案报价,让我们帮助您实现业务目标

立即咨询