固定翼六自由度仿真配平工具箱与批量扫描脚本实战

我在调固定翼六自由度仿真程序的时候遇到最多的问题不是控制器参数不是气动数据而是最基础的一步配平。模型建好之后初始状态随手给个迎角油门给个30%升降舵给个0一按运行飞机要么在地面弹两下要么瞬间抬头拉成一个夸张的爬升角。后来我把配平逻辑做成一套独立的工具箱再配一个批量扫描脚本这个问题才算彻底解决。这篇文章就把这套配平工具箱加脚本的完整思路拆开讲一遍从方程怎么写、函数怎么封装、脚本怎么扫描到配平结果怎么往下接尽量都讲透。这里说的工具箱其实有两层含义。一层是MATLAB自带的优化工具箱fsolve、lsqnonlin 这几个求解器就是配平计算的引擎另一层是我自己封装的配平函数库把气动模型、大气环境、残差方程全部模块化后续换飞机模型、换工况都不用改主逻辑。脚本则是用来做批量扫描的比如把速度从 50 m/s 扫到 300 m/s把高度从海平面扫到一万米一次性生成全包线配平数据表。这套组合适合刚把仿真模型搭起来、却卡在“飞机怎么飞都飞不稳”的人也适合想了解配平背后原理、不想只点按钮看结果的人。1. 为什么要做配平仿真稳定飞行的第一道关卡1.1 配平到底在配什么配平这个词听起来很玄但说白了就一句话在给定的飞行条件下找一组状态量和控制量让飞机受到的合外力为零、合外力矩也为零。状态量包括速度、迎角、俯仰角、角速度控制量就是油门位置和舵面偏度。满足这个条件的点就是一个配平点也叫平衡点。打个比方你手里捧着一杯水想让水面完全静止不动手必须保持一个稳定的姿势既不能抖也不能慢慢倾斜。飞机在天空中的定直平飞也是一样升力必须等于重力推力必须等于阻力俯仰力矩必须为零。如果这三个条件不满足飞机就会一直加速、掉高度或者抬头低头姿态根本稳不住。很多初学者容易忽略一个细节平衡和配平不是一回事。飞机在爬升或者转弯的时候力也可能平衡但力矩不为零或者运动轨迹不是定直平飞。配平通常特指定直平飞这种最基础的平衡状态也就是迎角、速度都不随时间变化飞行轨迹是一条直线。1.2 配平在仿真程序里的位置仿真程序建好之后并不是直接丢一组初值就能跑。六自由度运动方程是非线性微分方程组初始状态必须满足力平衡和力矩平衡否则积分开始后飞机就会出现剧烈的瞬态响应。你可能会看到俯仰角速度立刻飙到每秒几十度或者迎角直接发散这时候你根本分不清是气动模型错了、控制器写错了还是只是初始状态没配平。所以配平在仿真流程里的位置是在“模型建模”和“动态仿真”之间。它的输出是一组基准状态后续的非线性仿真从这组状态开始才能观察真实的动态特性。线性化更需要配平点因为在配平点附近做小扰动展开才能得到状态空间矩阵。换句话说配平做不好后面飞控设计、模态分析、操纵品质评估全都无从谈起。2. 配平的数学模型先把方程组写对2.1 纵向运动配平方程组固定翼飞机的完整配平分纵向和横航向。纵向配平处理速度、迎角、俯仰姿态的问题横航向配平处理侧滑、滚转、偏航的问题。大多数飞机对称布局横航向在零侧滑、零副翼、零方向舵的对称条件下自然满足平衡所以纵向配平是第一步也是最核心的一步。纵向定直平飞条件下假设无风、零侧滑、飞机对称飞行沿速度方向和垂直速度方向可以写出两个力平衡方程T * cos(α) - D - W * sin(γ) 0T * sin(α) L - W * cos(γ) 0其中 T 是推力D 是阻力L 是升力W 是重力α 是迎角γ 是航迹倾角。定直平飞时 γ 0方程简化为T * cos(α) - D 0T * sin(α) L - W 0再加上俯仰力矩平衡也就是绕飞机重心的俯仰力矩 M 等于零M M_thrust 0M 是气动俯仰力矩M_thrust 是推力产生的俯仰力矩。如果推力线正好通过重心M_thrust 就是零但很多飞机发动机安装位置有一定力臂这个力矩必须考虑。2.2 气动数据如何进入残差方程升力、阻力、俯仰力矩都来自气动模型。工程中最常见的形式是风洞数据或CFD数据的插值表模型的输入是迎角、侧滑角、舵面偏度以及马赫数、动压等输出是气动系数。为了让这篇文章的代码可复现我这里用一个经典的线性化气动导数模型做演示C_L C_L0 C_Lα * α C_Lδe * δeC_D C_D0 K * C_L^2C_m C_m0 C_mα * α C_mδe * δe然后用标准公式换算成力L qbar * S * C_LD qbar * S * C_DM qbar * S * cbar * C_m其中 qbar 0.5 * ρ * V² 是动压S 是机翼参考面积cbar 是平均气动弦长ρ 是大气密度。C_D 的抛物线阻力极曲线是一种常见的工程近似真实项目里还是应该用插值表但用来讲清楚配平的原理和代码逻辑完全够用。2.3 为什么纵向配平通常是三个未知量给定飞行速度 V 和高度 h 之后大气密度和动压就确定了。再看纵向配平方程组未知量是迎角 α、升降舵偏度 δe、油门 δT正好三个未知量三个方程组成一个完备的非线性方程组。状态量里其实还有俯仰角 θ 和角速度 q但定直平飞时 γ 0α 和 θ 是同一个问题q 等于零。所以配平求解本质上是“给速度、高度求 α、δe、δT”输出就是一组配平状态和控制指令。需要提醒的是油门到底代表多少推力不同发动机差别很大。有的项目直接建模成推力曲线 δT 与转速、速度、高度的函数有的项目简化为最大推力乘以油门杆位置。在示例代码里我采用 T δT * Tmax 的简化形式真实项目把这个地方替换成年发动机推力模型即可。3. 配平工具箱的设计与核心封装3.1 工具箱整体结构与模块划分我先说一下工具箱的模块划分这个结构用在我的几个项目里都比较顺。整体上分四个模块环境模块、气动模块、残差模块、求解主模块。环境模块负责给定高度下的密度、音速、重力加速度最简单的实现是国际标准大气模型低空段直接用指数公式 rho rho0 * exp(-h / 8400) 也能凑合。气动模块负责输出升力系数、阻力系数、俯仰力矩系数输入是迎角和舵面偏度输出可以是系数也可以是插值表结果。残差模块负责把力平衡和力矩平衡整理成 F(x) 0 的形式。求解主模块负责调用 MATLAB 优化工具箱的 fsolve 或 lsqnonlin并把结果包装成结构体返回。这样拆的好处是以后换机型只改气动模块换大气模型只改环境模块配平主逻辑永远不动。3.2 核心残差函数 TrimResidual 的实现残差函数是整个配平程序的灵魂。它的作用是把一组待求解的变量也就是 α、δe、δT映射成三个残差值。三个残差越接近零说明这组解越接近配平状态。文件名字就叫 TrimResidual.mfunction R TrimResidual(x, V, h, m, S, cbar, g, Tmax) % 输入 x [alpha; delta_e; delta_T] alpha x(1); delta_e x(2); delta_T x(3); % 大气密度仅做演示用正式项目抽成单独函数 rho0 1.225; rho rho0 * exp(-h / 8400); qbar 0.5 * rho * V^2; W m * g; % 气动系数占位模型真实项目用插值表 CL0 0.2; CLalpha 5.2; CLde 0.5; CD0 0.025; K 0.045; Cm0 0.02; Cmalpha -1.2; Cmde -1.5; CL CL0 CLalpha * alpha CLde * delta_e; CD CD0 K * CL^2; Cm Cm0 Cmalpha * alpha Cmde * delta_e; % 简化推力模型油门 * 最大推力 T delta_T * Tmax; % 残差无量纲化方便求解器收敛 R zeros(3,1); R(1) (T * cos(alpha) - D) / (qbar * S); R(2) (T * sin(alpha) L - W) / (qbar * S); R(3) Cm; endD 和 L 在残差函数内需要用 qbar * S * CD 和 qbar * S * CL 计算上面的代码为了简洁没有单独拆开实际完整版本可以写成L qbar * S * CL; D qbar * S * CD; M qbar * S * cbar * Cm;有个比例尺度的细节值得单独说。刚写配平程序的时候我用 R(1) 直接放 Tcos(alpha)-D单位是牛顿量级可能是几万牛R(3) 如果放力矩 M单位是牛米量级可能是几十万。fsolve 处理这种跨量级的残差时有收敛变慢的风险所以我习惯把力残差除以 qbarS把力矩残差直接换成力矩系数 Cm让残差都落在零点几的量级。这样求解器跑起来稳定得多。3.3 配平主函数与求解器配置残差函数写好之后配平主函数就简单了。它负责把初值、边界、求解器参数组装起来调用 fsolve再把结果封装成结构体。下面这段代码是配平主函数 trimAircraft.m 的核心逻辑function trim trimAircraft(V, h, m) % 演示用飞机参数真实项目从这里接参数配置 S 30.0; cbar 3.5; g 9.80665; Tmax 80000; % 初值迎角 2.8 度升降舵 2 度油门 0.3 x0 [0.05; 0.03; 0.3]; opts optimoptions(fsolve, ... Display, iter, ... Algorithm, trust-region-dogleg, ... MaxFunctionEvaluations, 500, ... MaxIterations, 200, ... FunctionTolerance, 1e-10, ... StepTolerance, 1e-10); fun (x) TrimResidual(x, V, h, m, S, cbar, g, Tmax); [x, ~, exitflag] fsolve(fun, x0, opts); trim.valid exitflag 0; trim.alpha x(1); trim.delta_e x(2); trim.delta_T x(3); trim.V V; trim.h h; end这里有个容易被忽略的点初值 x0 的选择。fsolve 是局部算法初值离真实解太远很可能不收敛或者收敛到一个物理上不合理的解。比如迎角初值给 0.5 rad也就是 28 度在低速状态下可能接近失速区配平出来的解就完全没有意义。后文会专门讲初值的延续法策略。算法选择上trust-region-dogleg 是 fsolve 处理中小规模非线性方程组的默认推荐一般不用改。如果你手里的配平方程组规模变大比如加入横航向算法可以换 Levenberg-Marquardt两者的区别主要在搜索策略上实际效果需要试。3.4 有约束求解推荐使用 lsqnonlin 而不是 fsolve用 fsolve 的时候还容易遇到一种情况数学上收敛了但结果物理上不可用。比如油门解出来 -0.2或者升降舵偏度解出来 1.5 rad相当于 85 度这种解从数学角度看确实让残差等于零了但真实飞机根本不可能在这种状态下飞行。解决办法是把配平当成有约束的优化问题用 lsqnonlin 替代 fsolve。lsqnonlin 可以给变量设置上下界同时用最小二乘的方式让残差平方和最小即使没有严格零点也能给出一个近似配平点。适合配平使用的边界大致是alpha_lb -0.3; % 迎角下限约 -17 度 alpha_ub 0.5; % 迎角上限约 28 度 de_lb -0.5; % 升降舵下限约 -28 度 de_ub 0.5; dt_lb 0.05; % 油门下限慢车状态一般不是 0 dt_ub 1.0; lb [alpha_lb; de_lb; dt_lb]; ub [alpha_ub; de_ub; dt_ub]; opts_lsq optimoptions(lsqnonlin, ... Display, iter, ... MaxFunctionEvaluations, 500, ... MaxIterations, 200, ... FunctionTolerance, 1e-10, ... StepTolerance, 1e-10); [x, resnorm, residual, exitflag] lsqnonlin(fun, x0, lb, ub, opts_lsq);我在实际项目里基本固定使用 lsqnonlin 做配平。原因很简单配平问题天然带物理意义的边界油门不可能大于 1、升降舵不可能偏转 90 度无视这些边界跑出来一个数学解后续环节全都白做。4. 脚本化批量配平从单点解到全包线数据4.1 为什么必须脚本化单点配平函数写出来只能解决“给定一个状态求一个解”的问题。但飞机设计里经常要的是整条包线看不同速度、不同高度下的配平曲线趋势。升降舵配平偏度随速度怎么变化、油门指令随高度怎么变化、迎角在低速点是不是接近失速区这些信息都是靠批量扫描得到的。如果手动一个个点去配平效率低而且容易出错。尤其是改了一版气动数据之后全包线要重新生成一遍手点根本不现实。写一个脚本把所有要计算的工况排列组合跑一遍一次性输出数据表和趋势图才是正常做法。脚本批处理的另一个好处是方便发现气动模型的异常。正常飞机的配平曲线是光滑的如果扫描结果里某个速度点突然跳变大概率是气动数据在那个点有洞或者是插值表设置有问题。4.2 批量扫描脚本 runTrimSweep 的编写批量扫描脚本的核心逻辑是两层循环外层遍历高度内层遍历速度每个点都调用 trimAircraft 函数。为了加速计算Serial 循环已经够用不需要上并行除非包线点数特别多。下面是 runTrimSweep.m 的关键逻辑以固定海平面高度扫描速度为例V_vec 50:5:300; % 速度范围m/s h_fixed 0; % 高度m m 12000; % 质量kg n length(V_vec); alpha_arr NaN(1, n); delta_e_arr NaN(1, n); delta_T_arr NaN(1, n); exit_arr zeros(1, n); x_prev [0.05; 0.03; 0.3]; for i 1:n Vc V_vec(i); trim trimAircraft(Vc, h_fixed, m); if trim.valid alpha_arr(i) trim.alpha * 180 / pi; delta_e_arr(i) trim.delta_e * 180 / pi; delta_T_arr(i) trim.delta_T; exit_arr(i) 1; % 延续法技巧把这次解作为下一个速度点的初值 x_prev [trim.alpha; trim.delta_e; trim.delta_T]; else warning(V %.1f m/s 配平失败, Vc); end end延续法是这个脚本的精髓。配平计算最怕初值乱给但速度从低到高扫描时相邻两个速度点的配平状态相差很小把上一个点的解直接拿来做下一个点的初值收敛概率大大提升。这里 x_prev 的作用就是传递这个信息比每次都用固定初值稳得多。多高度扫描就是在外面再套一层高度循环比如 0 m、2000 m、4000 m、6000 m然后按高度分组存储结果。如果数组维度变得复杂建议定义一个结构体数组每个高度一个元素内部再放速度扫描结果。4.3 结果可视化与合理性质检跑完扫描之后第一件事不是看数据而是画曲线。至少画三张图配平迎角随速度变化、升降舵配平偏度随速度变化、油门配平指令随速度变化。如果是多高度扫描把不同高度用不同颜色画在同一张图上一眼就能看出趋势。figure; subplot(3,1,1); plot(V_vec, alpha_arr, o-); grid on; xlabel(V (m/s)); ylabel(alpha (deg)); title(配平迎角 vs 速度); subplot(3,1,2); plot(V_vec, delta_e_arr, s-); grid on; xlabel(V (m/s)); ylabel(delta_e (deg)); title(升降舵配平偏度 vs 速度); subplot(3,1,3); plot(V_vec, delta_T_arr, ^-); grid on; xlabel(V (m/s)); ylabel(delta_T); title(油门配平指令 vs 速度);从物理直觉检查低速时动压小要产生足够的升力迎角必然大升降舵配平偏度也通常偏向抬头方向高速时动压大迎角小升降舵配平偏度会偏向低头方向。油门则呈现一种近似 U 型的趋势低速时为了维持升力诱导阻力大推力需求高高速时零升阻力大推力需求也高中速段存在一个最小阻力点对应最省油的巡航速度。如果扫出来的曲线不符合这些基本趋势就要回头检查气动模型、动压计算或者质量参数了。5. 常见问题与排查技巧实录5.1 配平失败与不收敛速查表配平程序跑不收敛是家常便饭把常见的情况整理成了一张速查表排查起来比较快现象常见原因排查方向fsolve 直接报初始点值奇异雅可比矩阵为奇异通常是气动导数为零导致检查气动数据α 和 δe 是否进入方程用有限差分检查雅可比残差一直震荡不下降初值离真实解太远用延续法从巡航点逐步扫描收敛但迎角是负数初值里 α 给太小求解器落入另一个根检查载荷因数是否合理或用带边界 lsqnonlin收敛但油门大于 1无约束 fsolve 跑飞了改用 lsqnonlin 加边界低速点配平失败动压太小需要的 α 超出边界确认边界范围是否覆盖失速前区域检查是否已经超过可用升力高速点配平失败最大推力不足以克服阻力检查 Tmax 和阻力极曲线确认速度是否超过该飞机理论极速5.2 延续法与初值策略配平函数里最容易出问题的就是初值。我见过不少人每个工况都从 x0 [0; 0; 0.5] 开始算低速点运气好能收敛高速点经常直接发散。原因很简单高速配平状态和低速配平状态差别很大一个固定初值不可能覆盖所有工况。延续法是最实用的解决手段。具体做法是先把不可行的工况空间划分成一条路径从已知收敛点出发沿着速度轴或者高度轴一步步推进每一步都用上一步收敛的结果作为新初值。这相当于把一个大范围的搜索问题拆成了一串相邻的小范围搜索问题。配合规模大的时候还可以更进一步用该速度点附近的两个已收敛配平解做线性外插作为初值。比如知道 V 50 和 V 55 的配平解算 V 60 的初值就按斜率外推。这个方法在气动数据线性度好的区域非常稳但穿过非线性较强的区域时效果会打折扣需要结合实际情况判断。5.3 单位制、量纲与力矩参考点说几个我亲手踩过的坑每一个都能让配平结果诡异地离谱。第一个坑是角度单位。气动导数里 C_Lα 的单位是 1/rad如果你传进去的 α 是角度制比如 0.1 而不是 0.0017升力系数会大得离谱。我建议所有内部计算严格使用弧度制只在输入输出时转换角度制代码里用变量名 angle_rad 做区分。第二个坑是力矩参考点。俯仰力矩系数 C_m 的基准是平均气动弦的某个参考点一般是机翼焦点或者四分之一弦点。气动数据表里 C_m 的参考点和算力矩用的力臂必须保持一致否则你还在算配平实际上是拿一个拧着劲的力矩模型在算平衡结果自然不对。第三个坑是推力线。推力如果不是通过重心就需要额外叠加 M_thrust T * z_T其中 z_T 是推力线到重心的垂直距离。很多简化模型假设 z_T 0但真实飞机这个值可能是正负一到两米对配平影响相当大。我建议在气动模型和机体模型里显式配置推力力臂不要默默省略。第四个坑是质量。飞机在飞行中质量一直在变燃油消耗、载荷投放都会让配平点移动。很多算例直接用最大起飞质量跑全包线结果高速点根本不可能配平其实是质量取得太大了。配平常量可以用质量做输入参数方便按飞行阶段切换。6. 配平结果往下怎么用6.1 初始化非线性仿真配平结果最常见的出口就是初始化非线性六自由度仿真。跑动态响应之前把速度设成配平速度迎角设成配平迎角俯仰角设成配平迎角角速度全部归零油门设成配平油门升降舵设成配平偏度然后才开始积分。这样飞机的初始状态就在力平衡和力矩平衡点上动态响应完全由控制输入或者扰动激起。一个细节是俯仰角的设置。配平得到的是迎角 α定直平飞时航迹倾角 γ 0俯仰角 θ 等于 α。如果你初始化的时候把 θ 设成 0那飞机初速度方向就是水平的但机体系相对水平面有一个正迎角开始积分后重力分量会立刻产生一个向下的加速度飞机会快速掉高度。这个错误看起来不明显因为程序不会崩但结果曲线很难看。6.2 线性化与模态分析基准点配平点做线性化的意义在于非线性方程在配平点附近可以用线性系统近似。状态矩阵 A 和控制矩阵 B 就是在这个点对状态变量和控制量求偏导得到的。飞行力学里两个经典模态短周期和长周期就是从这个线性系统的特征值里读出来的。实际操作可以在每个配平点计算一次状态矩阵然后看特征值的分布。如果某个速度点特征值实部突然变正说明该点附近的动态是不稳定的可能是失速或颤振边界的前兆。如果配平点本身不对这些分析全部建立在错误的工作点上结论没有意义。6.3 控制律设计与飞行品质评估现代飞控设计基本都是围绕配平点展开的。LQR 也好H∞ 也好第一步都是选定工作点在工作点附近做线性化再针对线性模型设计控制器。增益调度控制的多个工作点本质上就是一系列配平点的集合。飞行品质评估也依赖配平点。比如 Cooper-Harper 评分对应的响应特性是在特定配平状态施加小扰动之后测出来的。如果配平状态本身不合理比如把迎角配在 25 度接近失速边界那评估出来的飞行品质必然很差但这不是飞机的真实毛病而是你选了一个不应该飞的工作点。在工程流程里我通常的做法是先生成一组速度-高度网格上的配平数据表这组表同时服务于非线性仿真、线性化分析和控制律调度。同一个数据源三处复用避免各算各的造成不一致。7. 最后分享一点配平经验说个我自己的体会。配平函数写起来不难真正难的是把它写进一个复杂的仿真工程里和一堆已有模块衔接好。我一开始想做一个全自动的包线扫描只要气动数据改一版就重新生成整个包线的配平曲线。后来发现自动化越早做越容易把气动数据的错误藏在一堆看起来正常但实际上不合理的曲线里。现在我改成先挑三个点手算验证一个低速大迎角点、一个巡航点、一个高速点三个点都符合物理直觉后再跑全包线扫描。配平这东西不是把 fsolve 调用通了就是会了而是要把单位、量纲、气动系数导数、力矩参考点这些底层细节全部嚼碎了。你先在自己的模型上手动推几步把配平过程的每个环节看清楚后面再去封装工具箱和批量脚本就会顺手得多。