☰
固定翼六自由度仿真配平工具箱与批量扫描脚本实战
2026/10/2 12:05:26 网站建设 项目流程

我在调固定翼六自由度仿真程序的时候,遇到最多的问题不是控制器参数,不是气动数据,而是最基础的一步:配平。模型建好之后,初始状态随手给个迎角,油门给个30%,升降舵给个0,一按运行,飞机要么在地面弹两下,要么瞬间抬头拉成一个夸张的爬升角。后来我把配平逻辑做成一套独立的工具箱,再配一个批量扫描脚本,这个问题才算彻底解决。这篇文章就把这套配平工具箱加脚本的完整思路拆开讲一遍,从方程怎么写、函数怎么封装、脚本怎么扫描,到配平结果怎么往下接,尽量都讲透。

这里说的工具箱其实有两层含义。一层是MATLAB自带的优化工具箱,fsolve、lsqnonlin 这几个求解器就是配平计算的引擎;另一层是我自己封装的配平函数库,把气动模型、大气环境、残差方程全部模块化,后续换飞机模型、换工况都不用改主逻辑。脚本则是用来做批量扫描的,比如把速度从 50 m/s 扫到 300 m/s,把高度从海平面扫到一万米,一次性生成全包线配平数据表。这套组合适合刚把仿真模型搭起来、却卡在“飞机怎么飞都飞不稳”的人,也适合想了解配平背后原理、不想只点按钮看结果的人。

1. 为什么要做配平:仿真稳定飞行的第一道关卡

1.1 配平到底在配什么

配平这个词听起来很玄,但说白了就一句话:在给定的飞行条件下,找一组状态量和控制量,让飞机受到的合外力为零、合外力矩也为零。状态量包括速度、迎角、俯仰角、角速度,控制量就是油门位置和舵面偏度。满足这个条件的点,就是一个配平点,也叫平衡点。

打个比方,你手里捧着一杯水,想让水面完全静止不动,手必须保持一个稳定的姿势,既不能抖,也不能慢慢倾斜。飞机在天空中的定直平飞也是一样,升力必须等于重力,推力必须等于阻力,俯仰力矩必须为零。如果这三个条件不满足,飞机就会一直加速、掉高度或者抬头低头,姿态根本稳不住。

很多初学者容易忽略一个细节:平衡和配平不是一回事。飞机在爬升或者转弯的时候,力也可能平衡,但力矩不为零,或者运动轨迹不是定直平飞。配平通常特指定直平飞这种最基础的平衡状态,也就是迎角、速度都不随时间变化,飞行轨迹是一条直线。

1.2 配平在仿真程序里的位置

仿真程序建好之后,并不是直接丢一组初值就能跑。六自由度运动方程是非线性微分方程组,初始状态必须满足力平衡和力矩平衡,否则积分开始后飞机就会出现剧烈的瞬态响应。你可能会看到俯仰角速度立刻飙到每秒几十度,或者迎角直接发散,这时候你根本分不清是气动模型错了、控制器写错了,还是只是初始状态没配平。

所以配平在仿真流程里的位置,是在“模型建模”和“动态仿真”之间。它的输出是一组基准状态,后续的非线性仿真从这组状态开始,才能观察真实的动态特性。线性化更需要配平点,因为在配平点附近做小扰动展开,才能得到状态空间矩阵。换句话说,配平做不好,后面飞控设计、模态分析、操纵品质评估全都无从谈起。

2. 配平的数学模型:先把方程组写对

2.1 纵向运动配平方程组

固定翼飞机的完整配平分纵向和横航向。纵向配平处理速度、迎角、俯仰姿态的问题,横航向配平处理侧滑、滚转、偏航的问题。大多数飞机对称布局,横航向在零侧滑、零副翼、零方向舵的对称条件下自然满足平衡,所以纵向配平是第一步,也是最核心的一步。

纵向定直平飞条件下,假设无风、零侧滑、飞机对称飞行,沿速度方向和垂直速度方向可以写出两个力平衡方程:

T * cos(α) - D - W * sin(γ) = 0

T * sin(α) + L - W * cos(γ) = 0

其中 T 是推力,D 是阻力,L 是升力,W 是重力,α 是迎角,γ 是航迹倾角。定直平飞时 γ = 0,方程简化为:

T * cos(α) - D = 0

T * sin(α) + L - W = 0

再加上俯仰力矩平衡,也就是绕飞机重心的俯仰力矩 M 等于零:

M + M_thrust = 0

M 是气动俯仰力矩,M_thrust 是推力产生的俯仰力矩。如果推力线正好通过重心,M_thrust 就是零,但很多飞机发动机安装位置有一定力臂,这个力矩必须考虑。

2.2 气动数据如何进入残差方程

升力、阻力、俯仰力矩都来自气动模型。工程中最常见的形式是风洞数据或CFD数据的插值表,模型的输入是迎角、侧滑角、舵面偏度,以及马赫数、动压等,输出是气动系数。为了让这篇文章的代码可复现,我这里用一个经典的线性化气动导数模型做演示:

C_L = C_L0 + C_Lα * α + C_Lδe * δe

C_D = C_D0 + K * C_L^2

C_m = C_m0 + C_mα * α + C_mδe * δe

然后用标准公式换算成力:

L = qbar * S * C_L

D = qbar * S * C_D

M = 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.m:

function 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; end

D 和 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 调用通了就是会了,而是要把单位、量纲、气动系数导数、力矩参考点这些底层细节全部嚼碎了。你先在自己的模型上手动推几步,把配平过程的每个环节看清楚,后面再去封装工具箱和批量脚本,就会顺手得多。

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

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

立即咨询