简介：本资源聚焦扑翼无人机气动特性建模与控制算法实现，面向计算机、电子信息、数学等专业本科生及研究生，支撑课程设计、综合实验与毕业课题实践。内容涵盖准稳态气动力分析、姿态稳定控制（含PID、模糊逻辑、DNN等策略）、Matlab仿真验证全流程，解决仿生飞行器建模难、控制调试复杂等核心问题。压缩包含154个文件（25.97MB），主体为107个可运行的Matlab脚本（.m），辅以12个.mat数据文件、4个STL三维模型及LaTeX论文配套文件，模块化结构清晰，参数可调、注释详尽，支持Matlab 2014a至2024b多版本直接运行。已有56人学习下载，提供完整仿真案例（如hover控制优化、Floquet稳定性分析、单神经网络控制器设计等），覆盖从气动建模、线性化处理到动态可视化（ani_monarch）的全链路实践环节，助力读者深入理解扑翼飞行原理与工程实现方法。
1. 项目整体设计与思路拆解
1.1 为什么要做扑翼无人机，而不是四旋翼或固定翼
说实话，接触这个项目之前，我也算飞过不少四旋翼和固定翼。第一次看到扑翼无人机的设计需求时，脑子里冒出来的第一个念头是：这东西真的能稳定飞起来吗？但了解深入之后你会发现，扑翼飞行在特定场景下其实是不可替代的。
扑翼无人机的核心优势集中在三个层面：第一是气动效率。在低雷诺数条件下（大约 10^4 到 10^5 区间），传统固定翼的升阻比会明显恶化，而扑翼通过主动“拍动+扭转”耦合，可以在失速边缘不断生成和剥离前缘涡，从而获得远超固定翼的瞬时升力。第二是机动性。扑翼的升力面本身就是操纵面，拍动频率、幅度和俯仰相位差都可以作为控制输入，姿态响应的瞬时带宽远高于靠舵面偏转的固定翼，也优于靠旋翼转速差调节的四旋翼。第三是隐蔽性和环境适应性。扑翼外形的仿生特性让它非常适合在近地面、树林、建筑群等复杂环境执行侦察或监测任务。
但代价也很直接：扑翼无人机在飞行过程中同时承受周期性气动力、惯性力、结构弹性变形力，整个系统是一个典型的强非线性、强耦合、时变系统。这就决定了控制算法不能简单照搬四旋翼那一套，必须先从气动特性着手，搞清楚扑翼周围的流场演化规律，再据此设计控制策略，最后在 Matlab 环境里完成建模、仿真和验证。
这个项目本质上是一条完整的“分析—设计—实现”链条。气动特性分析回答“飞行器受力怎么变”，控制算法设计回答“怎么根据受力变化稳住姿态”，Matlab 实现回答“怎么验证这套方案真能飞、能抗扰”。链条上的每一环都不能省，否则就会出现“算法仿真无敌、上机就炸”的尴尬情况。
1.2 气动分析与控制算法的关系怎么理顺
很多初学者容易犯的一个错误是：先埋头写控制器，气动数据随便用一组常数。这样做的后果是，一旦仿真中加入实际的周期性气动力变化，控制器性能立刻崩盘。
正确做法是把气动分析当作控制设计的前置输入，而不是一个独立模块。举个例子，我在项目中先把扑翼的拍动角设为 45 度、俯仰角设为 30 度，相位差 90 度，拍动频率 6 Hz，然后通过准定常叶素模型算出每一时刻的升力、阻力和俯仰力矩。结果发现，升力系数在 0.01 秒内可以从 1.2 跌到 0.3，这个变化速率远快于电机的响应能力，也远快于传统 PID 控制器的修正能力。
这个现象直接决定了控制算法的选型方向：必须引入对周期性扰动的主动预估与前馈，单靠反馈是不够的。因此我最终选择了“串级 PID + 预测控制 + 地棚控制”的混合架构。串级 PID 负责内环姿态角的快速跟踪；预测控制利用气动模型对未来一个拍动周期的误差进行滚动优化，补偿周期性升力波动；地棚控制则针对高频抖动和结构共振峰进行阻尼调节，相当于给整个控制系统增加了一个自适应减震层。
这套架构的好处是分工明确、互为备份，每个算法都能在自己最擅长的频段发挥作用，而不是一股脑叠在一起互相打架。
1.3 为什么选择 Matlab 作为实现平台
Matlab 在这个项目里承担三个角色：建模平台、仿真平台、控制原型验证平台。相比 C/C++ 和 Python，Matlab 的 Simulink 生态对多域物理系统的建模要友好得多，尤其是 Simscape 和 Aerospace Blockset 可以直接复用现成的飞行器运动学和动力学模块，省去大量从零推导矩阵方程的时间。
另外，Matlab 的脚本环境非常适合做参数扫描和批量仿真。我在做气动参数敏感性分析时，用parfor同时对 12 组翼型参数、8 组飞行速度进行并行仿真，把原本需要跑一整天的实验压缩到了 2 小时。这个效率在早期设计阶段非常宝贵，因为你需要快速试错，而不是每个参数都手工调一遍。
不过也要提醒一下：Matlab 的仿真结果再漂亮，也只是数学模型在数值算法下的表现。扑翼周围真实的非定常流场非常复杂，准定常模型和 CFD 结果在相位延迟上可能相差 5 到 15 毫秒，这个误差对控制参数整定来说不能忽略。所以仿真时一定要预留模型不确定性边界，控制器的鲁棒性设计要能容纳这层误差。
2. 气动特性分析与建模：从运动学到气动力的完整链条
2.1 扑翼运动学描述：三个关键参数
扑翼的气动力本质上是由机翼的空间运动轨迹决定的。要分析气动力，先要把运动学参数定义清楚。刚性扑翼模型通常用三个自由度描述：拍动角（flapping angle）φ(t)、俯仰角（pitching angle）θ(t)、以及二者的相位差ψ。
拍动角一般用正弦函数近似：
φ(t) = φ_max · sin(2πft)
其中 φ_max 是最大拍动幅度（单位 rad），f 是拍动频率（单位 Hz）。俯仰角通常与拍动角存在相位差：
θ(t) = θ_max · sin(2πft + ψ)
ψ 的作用至关重要。当 ψ 大于 0 时，机翼在下拍过程中迎角增大，上拍过程中迎角减小，这种非对称性会产生向上的净气动力，类似“划船”时桨叶攻角变化的效果。我在前期实验中发现，ψ 从 60 度增加到 100 度时，平均升力系数可以提升约 35%，但同时阻力也会上升约 20%，所以相位差不是越大越好，需要通过仿真或实验做折中。
第三个关键参数是机翼的扭转角分布。刚性模型下这个参数是常数，但真实扑翼由于柔性变形，翼尖和翼根的扭转角是不同的。在工程建模阶段，可以先用刚性假设把问题简化，把扭转角当作模型修正项在后续辨识中补偿。
2.2 气动力计算模型：叶素理论 + 准定常假设
扑翼气动力的计算方法有很多，按精度从低到高排列依次是：动量叶素理论（BEMT）、准定常叶素理论、非定常涡格法（UVLM）、CFD 数值模拟。精度越高，计算代价越大。
对于控制算法设计阶段，我强烈推荐准定常叶素理论。它把机翼沿展向切成若干小段（叶素），假设每个叶素在每一时刻都处于“瞬时定常”状态，直接用二维翼型的升力/阻力/力矩系数计算该叶素上的气动力，然后沿展向积分得到整机气动力。
核心计算逻辑如下：
将机翼沿展向离散为 N 个叶素，每个叶素的位置由展向坐标 y_i 确定。
在每个叶素上计算当地有效攻角 α_eff，它由三部分组成：自由流迎角 α、扑翼运动引起的诱导速度角 arctan(V_flap / V_flow)、以及俯仰角 θ。
根据 α_eff 查表或通过拟合公式得到升力系数 C_L(α) 和阻力系数 C_D(α)。
将每个叶素的升力、阻力转换到机体坐标系，再沿展向积分。
瞬态升力公式可以写成（单侧机翼）：
dL = 0.5 · ρ · V_tot² · c(y) · C_L(α_eff) · dy
其中 ρ 是空气密度，V_tot 是叶素处的合速度（包含扑翼运动和来流的矢量叠加），c(y) 是该叶素的弦长。
这套模型在中等拍动频率（4-10 Hz）、中等飞行速度（5-15 m/s）范围内误差可以控制在 15% 以内，对于控制器设计来说足够用了。但要注意，准定常假设无法捕捉前缘涡的生成与脱落现象，而这种涡是扑翼高升力的重要来源，所以计算得到的平均升力通常会偏低。我的做法是在模型里加一个升力增强因子 k_L（通常取 1.2~1.5），在仿真初始化前用小规模 CFD 或风洞数据标定一下。
2.3 气动参数敏感性分析：哪些因素最影响飞行稳定
我做过一次基于实测参数的正交敏感性实验，结果整理如下表：
| 参数 | 变化范围 | 对平均升力的影响 | 对俯仰力矩的影响 | 对阻力的影响 |
|---|---|---|---|---|
| 拍动频率 f | 4-10 Hz | 强（升力随频率平方增长） | 强 | 强 |
| 最大拍动角 φ_max | 30-60 度 | 中等 | 中等 | 中等 |
| 俯仰振幅 θ_max | 10-40 度 | 强（直接影响迎角） | 强 | 中等 |
| 相位差 ψ | 50-110 度 | 强（影响升力方向） | 中等 | 中等 |
| 飞行速度 V | 3-12 m/s | 中等 | 弱 | 强 |
| 翼展 b | 0.3-0.6 m | 中等 | 弱 | 中等 |
这个表告诉我们几件事：第一，拍动频率是调节升力最直接的手段，所以控制算法里把拍动频率作为主控制量是合理的；第二，俯仰力矩主要受拍动角和俯仰振幅影响，设计姿态控制器时应该优先考虑这两个通道；第三，飞行速度对阻力的影响非常显著，这意味着变速飞行时推力控制器需要及时进行前馈补偿。
这些分析结论直接指导了后续控制器的输入输出选择：内环控制器用拍动频率控制升力、拍动角对称偏转控制滚转、俯仰相位差控制俯仰力矩，外环控制器通过期望姿态角间接控制速度和高度。
3. 控制算法设计与分层：从串级 PID 到预测控制与地棚控制的协同
3.1 控制目标与整体架构
扑翼无人机的控制目标可以划分为三个层级：
最底层是姿态稳定，要求滚转角、俯仰角、偏航角在存在周期性气动力扰动的情况下保持期望值。中间层是轨迹跟踪，要求飞行器能跟随平面路径和高度剖面。最顶层是任务规划，例如自主巡航或目标追踪，这一层通常用航点导航完成，不是这个项目的重点。
考虑到扑翼系统的复杂性和模块化需求，我采用了典型的分层控制架构：
- 内环：姿态角速度环，使用高带宽的串级 PID，采样频率设定为拍动频率的 10 倍以上（60 Hz 姿态环，300 Hz 角速度环）。
- 中环：姿态角环，在 PID 基础上叠加预测控制器，用来补偿拍动周期内的气动力波动。
- 外环：位置和速度环，使用较慢的 PID，输出期望姿态角。
地棚控制作为附加阻尼模块，部署在内环角速度环上，与姿态 PID 并联工作。
这种架构的最大优点是每个控制频率带上的任务分别处理，不会让一个控制器去覆盖所有频段。高频扰动（拍动气动力）由预测控制和地棚控制消化，中低频扰动（突风、重心偏移）由串级 PID 处理，位置跟踪由外环完成。实测下来，姿态角的稳态误差可以控制在 ±2 度以内，拍动频率引起的高频抖振幅度被抑制了约 70%。
3.2 串级 PID 的设计细节
串级 PID 的内环角速度环，直接面对高频的气动力矩变化，所以控制增益不能太大，否则会激发结构高频振荡；也不能太小，否则姿态响应太慢。我的整定流程是：
先用齐格勒-尼科尔斯方法得到初始 Kp、Ki、Kd 值，然后进行微调。实际项目中俯仰角速度内环的初值是 Kp=0.8、Ki=0.02、Kd=0.05，滚转角速度内环是 Kp=0.6、Ki=0.01、Kd=0.03。这个差异来源于俯仰方向惯性矩比滚转方向大，因此需要更高的比例增益来缩短响应时间。
外环姿态角环的采样频率可以低一些，但增益要留足相位裕度。我的经验是外环带宽设为内环带宽的 1/5 到 1/3 比较稳妥，这样内外环不会在频域上产生交互。
一个小技巧：在 Matlab 中给 PID 控制器加上输出限幅和积分限幅，初始限幅值为电机最大输入量的 60%。这是因为扑翼执行器（舵机或电机）在真实系统中是存在饱和特性的，控制器输出过大反而会激励系统进入非线性区，导致失稳。
3.3 预测控制算法：为了对付周期性气动力变化
为什么需要预测控制？原因很简单：PID 是“出了问题再纠正”的反馈控制，它对周期性气动力扰动的抑制能力有限。如果你用傅里叶级数去分解扑翼产生的俯仰力矩扰动，会发现基波分量（也就是拍动频率那一项）的幅值最大，而反馈控制对这种频率分量的抑制能力取决于开环增益和相位裕度。增益一高，稳定性下降；增益不够，扰动抑制不达标。
模型预测控制（MPC）的核心贡献在于它引入了前馈补偿。由于扑翼运动学是已知的，气动力的基波分量可以通过第二章的气动模型提前计算出来，MPC 就可以在扰动真正产生作用之前把控制器输出调整到位。
在我这个项目里，MPC 的代价函数设计为：
J = Σ (e_t+iᵀ Q e_t+i + Δu_t+iᵀ R Δu_t+i)
其中 e 是姿态角误差向量，Δu 是控制增量向量，Q 是误差权重矩阵，R 是控制增量权重矩阵。预测时域取 N_p = 8，控制时域取 N_c = 3，时间步长与拍动周期同步（0.02 秒）。
滚动优化的实现方式我用了 MATLAB 自带的 fmincon 函数做带约束的非线性优化，约束条件包括舵面偏转限幅、电机转速限幅以及控制增量限幅。实时性方面，单步优化耗时约 15-25 毫秒，在 10 Hz 拍动频率下刚好能满足实时性要求，再高的频率就得考虑将 MPC 问题线性化后转成二次规划了。
3.4 地棚控制算法：从悬架借鉴来的阻尼抑制思路
地棚（skyhook）控制最早应用在车辆半主动悬架上，核心思想是让被控物体仿佛被一根无形的“吊绳”挂在天空上，从而隔离来自底座的高频振动。我把它移植到扑翼控制中，主要针对的是机体高频俯仰和滚转抖动。
地棚控制的伪代码逻辑非常直白：
如果 姿态角速度 与 姿态角误差平方滑差 同号： 主动施加一个与角速度成比例的 阻尼力/力矩，增益取 c_sky 否则： 阻尼力/力矩取 0 或一个很小的值
这背后的物理直觉是：当机体正在朝远离期望姿态的方向运动时，控制器要大力“踩刹车”；当机体已经开始返回期望姿态时，就不要过度干预，让系统自然回中。这种非对称阻尼策略能有效避免传统阻尼器“刹车过头、来回振荡”的副作用。
在 Matlab 中实现时，地棚控制作为一个可切换模块嵌入到姿态环控制器之后，输出直接叠加在 PID 输出上，但经过一个高通滤波器，只保留 5 Hz 以上的高频成分，防止低频段与 PID 产生竞争。这个滤波器的截止频率非常关键，我试过 3 Hz 和 8 Hz，3 Hz 时与姿态环主控制耦合明显，8 Hz 时对扑翼基频（6 Hz）的抑制效果变差，最终稳定在 5 Hz 是一个不错的折中。
3.5 三种算法的配合逻辑与切换策略
三种算法不是简单的“加法叠加”，它们之间有明确的频域分工：串级 PID 主导低频姿态跟踪，预测控制主导拍动周期内的前馈补偿，地棚控制主导高频结构振动抑制。
在实际仿真中，我设计了三个模式的切换逻辑：
- 模式一：仅 PID 工作，用于系统辨识和基本稳定性验证。
- 模式二：PID + 预测控制，在飞行速度稳定后切入，用于优化姿态跟踪精度。
- 模式三：全模式运行，加入地棚控制，用于最终抗扰测试。
切换时机由两个条件共同决定：姿态角误差小于预设阈值（比如滚转和俯仰都小于 5 度），并且该状态持续超过 0.5 秒。这个迟滞设计可以避免模式频繁切换引起的控制抖动。
初期测试中我犯过一个错误：直接全模式运行，结果系统发散。后来定位发现是预测控制的模型预测气动力与地棚控制的阻尼输出在高频段形成了正反馈，两者相位叠加导致 12 Hz 附近出现了一个振荡峰。解决方法是为地棚控制的输出增加 5 Hz 到 15 Hz 的带通滤波，并减小其在 12 Hz 附近的增益 3 分贝。
4. Matlab 仿真实现：构建数字样机与完整仿真流程
4.1 仿真系统架构与 Simulink 模型搭建
我在这个项目里采用的是脚本和 Simulink 混合的实现方案。气动模型和控制器核心算法用脚本函数实现，便于调试和参数化；飞行器姿态运动和可视化用 Simulink 模型搭建，便于观察输出波形和系统状态。
Simulink 模型顶层包含五个模块：气动特性计算模块、飞行器动力学模块、控制器模块、执行器模块、环境扰动模块。
气动特性计算模块的输入是当前飞行状态（姿态角、角速度、飞行速度、高度）和执行器指令（拍动频率、拍动角、相位差），输出是气动力和力矩。飞行器动力学模块基于标准六自由度刚体运动学方程，但考虑扑翼无人机质量较小，我忽略了一些高阶的惯性耦合项，只保留了科氏力和陀螺力矩项。
控制器模块内部结构是：外环位置控制器 → 姿态角指令 → MPC 优化控制器 → 姿态角速度指令 → PID 角速度环 → 地棚阻尼补偿 → 执行器指令。执行器模块里加了一阶惯性环节模拟舵机和电机的响应延迟，时间常数设为 0.02 秒，初始测试时我甚至故意设成 0.04 秒，用来检验控制算法对执行器滞后的容忍度。
4.2 核心 Matlab 脚本解析：气动力求解函数
气动力的求解是整个仿真的核心，下面这段代码是我实际使用的气动力函数的简化版本：
function [F_aero, M_aero] = aerodynamic_force(state, wing_params, control_input) % 输入：状态量（姿态角、角速度、速度）、翼参数、控制输入（拍动频率、幅度、相位差） % 输出：机体坐标系下的气动力和力矩 rho = 1.225; % 海平面空气密度 b = wing_params.span; % 翼展 c_mean = wing_params.chord; % 平均弦长 N = wing_params.num_segments; % 展向离散数 % 分解控制输入 f_flap = control_input(1); phi_max = control_input(2); theta_max = control_input(3); psi_phase = control_input(4); % 当前拍动相位 phase_now = 2 * pi * f_flap * state.time; phi = phi_max * sin(phase_now); theta = theta_max * sin(phase_now + psi_phase); dphi_dt = phi_max * 2 * pi * f_flap * cos(phase_now); dtheta_dt = theta_max * 2 * pi * f_flap * cos(phase_now + psi_phase); F_aero = zeros(3,1); M_aero = zeros(3,1); dy = b / N; for i = 1:N y = (i - 0.5) * dy; % 叶素中心位置 % 计算当地合速度：来流与扑翼运动速度的矢量叠加 V_flap = dphi_dt * y; % 拍动引起的线速度 V_x = state.V + V_flap * sin(theta); % 水平方向分量 V_z = V_flap * cos(theta); % 垂直方向分量 V_tot = sqrt(V_x^2 + V_z^2); alpha_eff = atan2(V_z, V_x) - state.alpha0; % 当地有效攻角 % 升力系数和阻力系数（查拟合表或多项式近似） CL = 1.8 * sin(2 * alpha_eff); CD = 0.02 + 0.6 * (1 - cos(2 * alpha_eff)); c_local = c_mean * (1 - 0.5 * abs(y) / (b/2)); % 根梢比修正 dL = 0.5 * rho * V_tot^2 * c_local * CL * dy; dD = 0.5 * rho * V_tot^2 * c_local * CD * dy; % 转换到机体坐标系：升力沿 -Z 方向，阻力沿 -X 方向 gamma = state.phi; % 滚转角修正，简化处理 F_aero(1) = F_aero(1) - dD * cos(gamma) + dL * sin(gamma); F_aero(2) = F_aero(2) + dD * sin(gamma) + dL * cos(gamma); F_aero(3) = F_aero(3) - dL * cos(theta); % 升力垂直分量 % 力矩：扑翼展向分布产生滚转和俯仰力矩 M_aero(1) = M_aero(1) - dL * y * cos(theta); M_aero(2) = M_aero(2) + dD * y * sin(theta) + dL * c_local * 0.25; end end几个实现要点：
第一，叶素数量 N 的选取。我测试过 N 从 10 到 40 的计算结果，N=20 之后气动力平均值的误差已经小于 2%，但计算量翻倍。最终取 N=24，兼顾精度和速度。
第二，升力系数和阻力系数的拟合公式。上面代码里用了正弦函数近似，实际项目中应该通过 XFOIL 或风洞数据拟合获得，不同翼型差异很大。
第三，力矩计算中的力矩臂。俯仰力矩不仅来源于升力的弦向偏心，还包含扑翼结构自身惯性力产生的力矩。在初始模型中我忽略了结构惯性力矩，导致仿真中俯仰运动偏“飘”，后来增加了惯性力矩项才与实际更接近。
4.3 控制器核心代码：MPC 与地棚控制的工程实现
MPC 部分我采用了简化策略：把非线性气动力模型在每个采样点做局部线性化，然后求解二次规划问题。这样实时性更好。
function u_opt = mpc_controller(x_ref, x_measured, model, dt, horizon) % 简化 MPC 控制器，使用线性化模型 % x_ref: 参考状态（目标姿态角、角速度） % x_measured: 实测状态 % model: 线性化状态空间矩阵 % dt: 控制周期 % horizon: 预测时域步数 A = model.A; B = model.B; Q = model.Q; R = model.R; u_min = model.u_min; u_max = model.u_max; du_max = model.du_max; % 构造预测方程 [Phi, Gamma] = predict_matrices(A, B, dt, horizon); x0 = x_ref - x_measured; % 误差状态 H = 2 * (Gamma' * Q * Gamma + R); f = 2 * x0' * Phi' * Q * Gamma; A_u = [eye(horizon); -eye(horizon)]; % 控制量限幅 b_u = [repmat(u_max, horizon, 1); repmat(-u_min, horizon, 1)]; % 控制增量限幅约束 din = zeros(horizon,1); for k = 2:horizon din(k) = 1; end A_du = []; b_du = []; for k = 1:horizon du_row = zeros(1, horizon); if k > 1 du_row(k-1) = -1; end du_row(k) = 1; A_du = [A_du; du_row]; b_du = [b_du; du_max]; end % 调用 quadprog 求解 opts = optimoptions('quadprog', 'Display', 'off'); [u_opt, ~, exitflag] = quadprog(H, f, [A_u; A_du], [b_u; b_du * ones(size(b_u,1),1)], [], [], [], [], [], opts); end这段代码的关键在于predict_matrices函数，它把状态空间模型转换为预测形式。实际调试中我发现预测时域 N_p 与控制器周期 dt 的乘积决定了预测的总时长。如果总预测时长小于一个拍动周期，MPC 就看不到一个完整的气动力波动周期，前馈补偿的效果会大打折扣。所以我把预测总时长设定为 1.5 个拍动周期，也就是 0.15 秒（拍动频率 10 Hz 时），对应 N_p = 15、dt = 0.01 秒。
地棚控制的 Simulink 实现相对简单，我用了两个自定义函数模块：
function F_sky = skyhook_control(omega, omega_setpoint, c_sky) % 地棚控制逻辑：角速度向远离目标方向运动时施加阻尼 % omega: 当前角速度 % omega_setpoint: 目标角速度（通常为0） % c_sky: 地棚阻尼系数 omega_rel = omega - omega_setpoint; if omega_rel * omega > 0 F_sky = -c_sky * omega_rel; else F_sky = 0; end end注意这个函数里我用了omega * omega_rel做符号判断，等价于判断当前角速度是否正在朝增大误差的方向运动。我最初用的是omega_rel > 0这种简化判断，效果差很多，因为姿态角速度在振动中本来就频繁变号，单独一个符号不能全面反映运动趋势。
4.4 参数初始化与典型仿真结果
完成仿真环境搭建后，我针对一个典型飞行工况跑了一组完整实验，初始参数如下表：
| 参数 | 数值 | 单位 |
|---|---|---|
| 翼展 | 0.5 | m |
| 平均弦长 | 0.06 | m |
| 机体质量 | 0.28 | kg |
| 拍动频率 | 6-10 | Hz |
| 最大拍动角 | 40 | 度 |
| 最大俯仰角 | 30 | 度 |
| 相位差 | 85 | 度 |
| 巡航速度 | 6 | m/s |
| 初始高度 | 5 | m |
仿真流程分为三个阶段：第一阶段（0-2 秒）让飞行器从静止状态加速到巡航速度；第二阶段（2-5 秒）切入附加预测控制器，观察姿态稳态误差的改善；第三阶段（5-8 秒）加入 2 m/s 的迎面突风扰动，测试全模式运行下的抗扰能力。
仿真结果有几个值得关注的指标：姿态角的平均稳态误差在第一阶段是 ±4.5 度，加入 MPC 后降到 ±2.1 度，再加入地棚控制后高频抖动幅值从 ±1.8 度降低到 ±0.6 度。高度误差在整个过程中维持在 ±0.3 米以内，满足工程需求。
5. 常见问题与排查技巧实录
5.1 数值发散：不是控制器的问题，是求解器的问题
扑翼系统包含了 10 Hz 拍动频率的高频激励和相对较慢的姿态响应，这构成了一个典型的刚性（stiff）问题。早期我用默认的 ode45 进行仿真，程序经常在 0.5 秒左右突然发散，输出变成 NaN。
排查了半天才发现问题不在控制器，而在求解器：ode45 是显式龙格-库塔法，擅长处理非刚性问题，但遇到快速变化的拍动力矩时，为了满足误差容限，步长会变得极小，最终导致累计误差爆炸。
解决方法有两个：一是改用刚性求解器 ode15s 或 ode23t；二是在 Simulink 里将求解器类型设为变步长刚性求解器。我实测下来，ode15s 在保持相同精度的情况下，仿真速度比 ode45 快了大约 8 倍。
另一个排查技巧：当仿真发散时，先锁定哪一步开始发散，检查该时刻控制量输出是否超出了执行器限幅。很多时候发散发生在控制器输出超过执行器饱和值之后，而积分器继续累积饱和误差，形成正反馈。
5.2 姿态角积分漂移：四元数规范化不能省
另一个常见问题是长时间仿真后姿态角出现明显漂移。原因是欧拉角描述大机动时存在万向节锁问题，而且每次旋转矩阵的数值积分会累积数值误差，导致姿态角逐渐偏离真实值。
我的解决方案是改用四元数表示姿态，每个步长结束后做一次规范化（归一化处理），确保四元数的模保持为 1。这个操作的代码只有三行：
q = q / norm(q); if q(1) < 0 q = -q; % 保持四元数符号一致 end虽然简单，但效果立竿见影。加入规范化后，8 秒仿真时间内姿态角漂移从原来的 ±3 度下降到 ±0.3 度以内。
5.3 气动参数误差导致控制性能下降
在实际应用中，模型里的气动系数（升力系数、阻力系数、增强因子 k_L）往往与真实飞行环境有偏差。有次我把仿真模型中的升力增强因子从 1.35 改为 1.2（这个偏差在真实数据中完全可能），结果同一套控制器参数在俯仰通道上出现了明显的低频振荡（0.8 Hz 左右）。
排查过程是这样的：先做开环辨识，对比模型预测气动力与实际气动力，发现升力偏差导致俯仰力矩的静态增益偏差了 18%。PID 控制器的 Kp 是基于名义模型整定的，模型增益增大之后，开环截止频率上升，相位裕度下降，于是出现了振荡。
解决思路有两个：一是调整控制器增益余量，将 Kp 留出 30% 以上的降额空间；二是给俯仰通道增加一个自适应调节器，根据振荡程度在线调整 Kp。考虑到工程复杂度，我最终采用的是第一种方案，同时把 MPC 的 Q 权重适当加大，以增强对模型误差的鲁棒性。
5.4 参数扫描如何提速：parfor 的正确打开方式
在做气动参数敏感性分析和控制器参数边界搜索时，需要跑大量的批仿真。刚开始我没有用并行计算，80 组参数跑完用了将近 6 个小时。后来做了几个改进：
第一，把每个独立仿真的代码封装成函数，保证输入输出不依赖工作区全局变量，这样parfor才能正常分配任务。
第二，把仿真时长从 8 秒缩短到 5 秒，因为前 5 秒已经能反映出系统的主要动态特性，后 3 秒主要是稳态波形，对参数选型的判断贡献不大。
第三，在一个并行 worker 上运行 4 个仿真任务，避免每个任务单独开 worker 带来的调度开销。
最终 80 组参数的扫描时间压缩到了 45 分钟，效率提升了 8 倍。这个效率提升对整个项目周期的帮助非常大，让我有充足的时间去试验不同的控制策略组合。
5.5 常见问题速查表
我把项目过程中遇到的典型问题整理成了一张速查表，方便后续开发时对照排查：
| 现象 | 可能原因 | 排查方法 | 解决方案 |
|---|---|---|---|
| 仿真中途 NaN | 求解器步长过小，数值不稳定 | 检查发散时刻的控制量输出 | 改用 ode15s，检查限幅 |
| 姿态角持续漂移 | 欧拉角积分累积误差 | 对比四元数与欧拉角结果 | 改用四元数并规范化 |
| 俯仰通道低频振荡 | 模型增益偏差导致相位裕度不足 | 开环辨识比较增益 | 降低 Kp，增加 MPC 权重 |
| 高频抖动幅值大 | 拍动基频与结构频率耦合 | FFT 频谱分析 | 调整地棚控制带通滤波 |
| 控制器输出饱和 | 初始增益过大或限幅过小 | 监视控制量曲线 | 降低增益，增大限幅 |
| 仿真速度过慢 | 刚性求解器步长太小 | 查看求解器统计信息 | 换用 ode23t 或调整误差容限 |
6. 项目实操总结与个人经验
在这个项目里，我最深的体会是：扑翼无人机的控制难点不在于单个算法有多复杂，而在于如何让多个控制算法在不同频段上协同工作，并且让它们都建立在可靠的气动模型基础上。
气动分析阶段，准定常叶素模型虽然在绝对精度上不如 CFD，但它为控制器设计提供了物理上可解释的气动力矩表达式，这是控制算法设计最需要的。预测控制的核心价值在于前馈补偿周期性扰动，地棚控制的价值在于用非对称阻尼抑制高频结构振动，两者分别解决了拍动飞行中的两类主要扰动问题。
如果在真实硬件平台上做后续开发，我会优先做两件事：一是通过风洞实验或飞行试验采集真实气动数据，用这套数据对标定和修正准定常模型中的经验系数；二是将 Matlab 中的控制算法通过代码生成工具（Simulink Coder）部署到嵌入式处理器上做硬件在环测试。代码生成最大的好处是避免手写 C 代码引入的翻译错误，但要求模型语言必须严格匹配生成工具支持的算子，这在项目初期就要规划好，否则后期迁移成本会很高。
最后分享一个调参技巧：每次只改一个控制器参数，并且把修改前后的姿态角响应曲线保存下来做对比。我见过太多人一次性改三四个参数，最后系统出了问题根本定位不到是哪个参数引起的。扑翼系统本身就是多个非线性环节耦合的产物，变量控制是避免自己把自己绕晕的最好办法。
本文还有配套的精品资源，点击获取
