ARTICLE DETAIL

资讯详情

深耕网站建设与运营推广的一线实战洞察。

固定翼仿真配平实战:用工具箱+脚本稳定求解平衡点

固定翼仿真配平实战:用工具箱+脚本稳定求解平衡点 自己做过固定翼模型飞行仿真的人应该都遇到过这种场景飞机建模、气动系数、舵面效率都搞定了满怀信心跑一次动态仿真结果飞机刚起飞就一头扎下去或者抬头失速。这时候很多人的第一反应是“控制器参数没调好”于是疯狂调PID折腾一晚上其实根本没触到问题的根源——你的飞机连配平点都没找到控制器再强也拉不住一个天生不平衡的对象。我最早接触固定翼仿真程序时就在这一步被卡了将近两周。后来发现配平这件事本质上是一个有约束的多元方程求解问题而现在的仿真工具链里既有专门的航空工具箱也有通用优化工具箱配合脚本完全可以实现系统化的求解。这篇文章就是记录我怎么用“工具箱脚本”的方式把一个固定翼仿真程序从“跑起来”做到“配得平”。内容适合刚开始做飞行仿真、或者被配平计算反复折磨的同行参考。我会先讲清楚配平到底在算什么再给出一个可以落地的仿真程序框架然后一步步演示如何用工具箱和脚本解配平最后聊聊那些代码之外、很容易让人抓狂的坑。1. 配平到底在平衡什么先搞懂力和力矩的账1.1 固定翼的六自由度运动和“稳定直线飞行”的条件固定翼飞机在空中的运动完整描述需要六个自由度三个位置分量和三个姿态角。仿真程序里通常用一组非线性微分方程来描述机身线速度、角速度、欧拉角的位置关系。你不需要时时刻刻全部关注配平关心的是“力的平衡”和“力矩的平衡”。所谓配平就是在给定的飞行状态比如定常平飞、定常爬升、定常转弯下找出一组状态量迎角、侧滑角、角速度和控制量升降舵、副翼、方向舵、油门的组合使得飞机的线加速度和角加速度都为零。用大白话说飞机在该状态下保持稳定不需要额外的动态调整松手也不会立刻飞出原有的飞行轨迹。对固定翼而言最常见的就是对称平面内的纵向配平。它要求升力与重力的分量平衡也就是法向力方程成立推力和阻力平衡也就是切向力方程成立俯仰力矩等于零也就是仰角不再变化。很多人在仿真里给飞机一个初始迎角和舵面偏度然后直接跑时间历程期望它自己收敛到某个平衡状态。这在物理上不是不可能但收敛速度取决于飞机的静稳定性而且初始条件稍微偏一点非线性模型很容易发散。正确的做法是先用配平计算求出那个点再用这个点作为动态仿真的初始条件。1.2 升降舵、油门、迎角之间怎么凑成一个平衡点纵向配平的未知量一般包括迎角 (\alpha)、升降舵偏度 (\delta_e)、油门位置 (\delta_T)。如果只要求定常平飞那这三个未知量要满足三个方程纵向力平衡、法向力平衡、俯仰力矩平衡。看起来简单但气动力是迎角和舵面的非线性函数推力和油门也有对应关系方程并不线性。举个例子很多教材里会给出小扰动线化公式[ C_L C_{L0} C_{L\alpha}\alpha C_{L\delta_e}\delta_e ][ C_m C_{m0} C_{m\alpha}\alpha C_{m\delta_e}\delta_e ]这是线性化的结果用在小扰动稳定性分析很好用。但真实气动数据在大迎角、大舵偏下往往有明显的非线性尤其是无人机或低速模型飞机的升力曲线不是一条直线失速前会有弯曲力矩系数也会随舵偏变化。这时候再用解析式手推配平点精度很难保证。1.3 为什么手动试凑经常发散需要对称求解手动试凑的思路通常是先定一个迎角算升力是否等于重力差多少就去改升降舵让力矩平衡力矩平衡后又发现升力变了再改迎角。这在接近线性区还有可能收敛但进入非线性区后经常来回振荡甚至越调越偏。因为你是在用一个手动的“迭代算法”而且没有考虑交叉耦合。配平问题的本质是求解一组形如 (f(x)0) 的非线性方程组其中 (x) 是状态量和控制量的组合。既然是非线性方程组用数值方法对称求解就比手动试凑可靠得多。程序里常用的手段有三类直接调用航空工具箱自带的配平函数例如MATLAB/Simulink中的trim用优化工具箱构造残差函数让残差平方和最小自己写牛顿-拉夫逊迭代。这三条路我都在项目中试过。最稳妥的组合是用航空工具箱求出初值再用优化工具箱处理边界约束。这也是后面要详细展开的内容。2. 仿真程序怎么搭从气动数据到可配平的飞机模型2.1 用一个最小可用固定翼模型起步不是非得建一个全量的高精度模型才能做配平。实际上配平逻辑和模型复杂度是解耦的只要你的模型能提供受力、力矩对状态和控制的响应就可以配平。我建议先搭一个“最小可用”的纵向模型包含质量与惯量质量 (m)、俯仰转动惯量 (I_y)气动升力、阻力、俯仰力矩系数推力模型重力模型。一般来说升力、阻力和俯仰力矩系数可以用表格形式给出插值计算。最常见的表达式是[ C_L C_L(\alpha, \delta_e, q) ][ C_D C_D(\alpha, \delta_e) ][ C_m C_m(\alpha, \delta_e, q) ]其中 (q) 是俯仰角速率阻尼项对动态稳定性有用但对稳态配平影响不大。把模型搭成Simulink还是写成纯函数取决于后续要做什么。如果只做配平纯函数最方便如果还要做动态仿真Simulink更直观。我的习惯是核心气动和推力写成MATLAB函数外面包一层可配置的数据结构让配平脚本和Simulink模型都调用同一套函数避免两边模型不一致。2.2 气动系数与推力模型气动数据可以从风洞试验、CFD计算或者公开的飞机模型数据中获取。对于学习用途可以用一组典型无人机气动系数比如参数项数值说明(C_{L0})0.2零迎角升力系数(C_{L\alpha})5.0 /rad升力线斜率(C_{L\delta_e})0.4 /rad升降舵效率(C_{D0})0.03零升阻力系数(k) (诱导阻力因子)0.04与升力系数平方相关(C_{m0})0.02零迎角俯仰力矩(C_{m\alpha})-0.8 /rad俯仰静稳定性导数(C_{m\delta_e})-0.6 /rad舵面俯仰力矩效率这些参数不是绝对的关键是要体现“力矩系数对迎角为负斜率、对升降舵有足够权限”这两个基本特征。推力模型我常用简单形式[ T \delta_T \cdot T_{\max} ](\delta_T) 范围是0到1(T_{\max}) 是最大静推力。如果螺旋桨飞机最好再考虑速度对推力的影响不过配平阶段可以先忽略等动态仿真时再细化。2.3 把模型放进仿真环境传感器和激励怎么加如果你的目标只是算配平点模型可以用函数计算。但如果后续要做传感器仿真、控制律验证就建议把飞机模型放到Simulink的标准框架里同时把空速管、陀螺仪、加速度计这些传感器模型也放进去。这里有个容易忽略的点传感器测量的物理量通常是机体坐标系下的数据而配平脚本计算用的可能是风轴系或大地坐标系坐标系转换一定要核对清楚。我见过很多次配平结果在纯代数计算里完全正确但接到Simulink后飞机还是“自己动了”最后发现是传感器反馈的迎角和真实迎角差了半个角度或者是坐标系旋转顺序写反了。所以在这篇文章的语境下我的建议是先在模型内部统一使用机体坐标系和风轴系传感器模型放在外围配平计算直接面向飞机模型的核心函数不要站在传感器输出端反推。3. 工具箱解配平的正确姿势trim函数的思路与优化工具箱的兜底3.1 先用Aerospace工具箱的trim函数如果你在用MATLABSimulink Control Design和Aerospace工具箱里提供了trim函数这是最直接的工具。trim函数的思路是为模型指定状态、输入和输出的约束然后寻找一个平衡工作点使得状态导数为零。使用前你需要先有一个Simulink模型并且明确模型中哪些是状态、哪些是输入、哪些是输出。举一个简化后的调用示例% 定义模型和操作点 model fixedwing_model; load_system(model); % 指定要配平的初始猜测 x0 [25; 0; 0; 0.05; 0; 0; 0; 0; 1000]; % 速度、姿态、角速度、位置等 u0 [0.05; 0.6; 0; 0]; % 升降舵、油门、副翼、方向舵 % 设置哪些状态在配平时保持不变或可变化 ix [1;4;5;6;7;8;9]; % 状态的哪些分量需要匹配 iu [1;2;3;4]; % 输入 iy []; % 输出 % 调用trim [op, report] findop(model, y0, [0; 0; 0; 0], x0, x0, u0, u0, ... states, ix, inputs, iu, outputs, iy);严格说findop会生成一个操作点对象。你需要查看输出report里的残差确认力与力矩平衡是否达到目标精度。问题在于trim函数默认用的是数值优化如果你的模型非线性很强或者初始猜测离真实配平点太远它可能收敛到局部解甚至直接报错。3.2 优化工具箱做配平的几个关键设计当trim不给力时优化工具箱就成了兜底方案。其实配平问题可以转化成一个无约束或带约束的优化问题已知状态和控制量的关系构造残差向量的二范数然后用fsolve或fmincon求解。这里的关键是残差向量怎么定义。对纵向配平我通常定义三个残差function [res] trim_residual(x, aircraft, flight_condition) % x [alpha; delta_e; delta_T] alpha x(1); de x(2); dt x(3); % 计算气动力和力矩 [L, D, M] aero_forces_moments(alpha, de, flight_condition.V); T dt * aircraft.Tmax; theta alpha flight_condition.gamma; % 爬升角为gamma时俯仰角为alphagamma % 纵向力方程机体坐标系 res(1) T * cos(alpha) - D - aircraft.m * aircraft.g * sin(flight_condition.gamma); res(2) T * sin(alpha) L - aircraft.m * aircraft.g * cos(flight_condition.gamma); res(3) M T * 0; % 若推力不产生俯仰力矩则只要求气动力矩为零 res res(:); end注意这里的力方程使用了风轴系和机体坐标系的混合概念实际使用时要严格统一。推力线如果不经过重心还要加上推力力矩项。然后用优化工具求解x0 [0.05; 0.0; 0.5]; options optimoptions(fsolve, Display, iter, ... Algorithm, trust-region-dogleg, ... MaxFunctionEvaluations, 200, ... SpecifyObjectiveGradient, false, ... FunctionTolerance, 1e-8); [x_trim, fval, exitflag] fsolve((x) trim_residual(x, aircraft, flight_condition), x0, options);如果存在约束比如舵面最大偏度、油门最大最小值就改用fmincon目标函数是残差的平方和objective (x) sum(trim_residual(x, aircraft, flight_condition).^2); nonlcon []; lb [-0.2; -0.6; 0.0]; ub [0.3; 0.6; 1.0]; options optimoptions(fmincon, Algorithm, sqp, ... Display, iter, OptimalityTolerance, 1e-8); [x_trim, J_min] fmincon(objective, x0, [], [], [], [], lb, ub, nonlcon, options);3.3 配平脚本框架拿来就能跑为了让整个过程可复现我建议写一个独立的配平脚本模块输入是飞机参数和飞行状态输出是配平点。这个脚本的骨架如下%% trim_script.m clc; clear; close all; % 1. 飞机参数结构体 aircraft.m 13.5; % kg aircraft.g 9.81; aircraft.S 0.82; % 机翼面积 m^2 aircraft.b 2.4; % 翼展 m aircraft.cbar 0.34; % 平均气动弦长 m aircraft.Iy 0.5; % 俯仰转动惯量 kg*m^2 aircraft.Tmax 80; % N % 2. 飞行状态 flight_condition.V 22; % m/s目标速度 flight_condition.rho 1.225; flight_condition.gamma 0; % 爬升角 rad % 3. 气动系数函数 aero init_aero_data(); % 装载插值表 % 4. 配平求解 x0 [0.04; 0.0; 0.5]; % alpha, delta_e, delta_T lb [-0.1; -0.5; 0.1]; ub [0.25; 0.5; 1.0]; options optimoptions(fmincon, Algorithm, sqp, ... Display, iter, StepTolerance, 1e-10, ... OptimalityTolerance, 1e-8, ... MaxFunctionEvaluations, 300); obj (x) sum(trim_residual(x, aircraft, flight_condition, aero).^2); [x_trim, fval, exitflag] fmincon(obj, x0, [], [], [], [], lb, ub, [], options); % 5. 输出 alpha_trim x_trim(1); delta_e_trim x_trim(2); delta_T_trim x_trim(3); fprintf(配平结果: alpha%.4f rad, delta_e%.4f rad, delta_T%.4f\n, ... alpha_trim, delta_e_trim, delta_T_trim); fprintf(残差: %.6e\n, fval);这套脚本的好处是飞机参数和气动数据可以随时替换换一架飞机只需要改数据不需要改求解逻辑。我后来在多个项目中都用这个框架从微型无人机到3米翼展的验证机改一改气动表就能用。4. 非线性配平的一些坑初始化、约束和收敛判定4.1 状态初值选错配平点会跑到天上配平计算本质上是非线性方程求根初值非常重要。你可能会想多给几个初值多点尝试不就行了是但这不是最优策略。更好的策略是先用一个简化的线性模型粗算一个接近的初值再用完整模型精算。比如从线性模型能得到近似迎角[ \alpha_0 \frac{W - T \sin\gamma}{\frac{1}{2}\rho V^2 S C_{L\alpha}} - \frac{C_{L0}}{C_{L\alpha}} ]然后用这个(\alpha_0)计算升降舵初值[ \delta_{e0} -\frac{C_{m0} C_{m\alpha}\alpha_0}{C_{m\delta_e}} ]这两个值是很好的初值。我自己踩过的坑就是直接用(0.1)弧度起步结果高迎角状态下气动表进入了非线性段fmincon收敛到另一个“数学上合理但物理上离谱”的点负舵面配平升力系数异常大看起来残差很小但对应的速度根本不是目标速度。所以在初值设计上永远要先从线性化公式出发再交给优化器。4.2 控制量边界和物理合理性fmincon可以带边界约束这是比单纯fsolve更适合配平的原因。升降舵偏度、油门开度都有物理限制。带上边界后残差为零的点可能落在边界外此时优化器会在边界上找到一个最小残差点而不是强行给出一个不可用的配平值。这里有个容易被忽视的问题如果残差在边界上无法接近零说明目标飞行状态本身不可配平比如要求的平飞速度低于该飞机的最低平飞速度或者超过最大速度。遇到这种情况优化器给的“最佳点”其实是在提示你改飞行状态而不是继续调边界。我建议在脚本里加一个判断if any(abs(x_trim - lb) 1e-6) || any(abs(x_trim - ub) 1e-6) warning(配平点落在控制量边界上请检查飞行状态是否可行); end这个警告能帮你省下很多排查时间。4.3 残差判定、雅可比矩阵和数值噪声配平收敛判定不能只看优化器报“optimal”要看物理残差的实际大小。fmincon里用OptimalityTolerance控制的是一阶最优性条件不是残差本身。我见过残差的模为(10^{-4})的情况看似很小但换算成角加速度可能仍然不可接受。最好在求解完成后单独计算一次力与力矩残差并归一化[ \epsilon \frac{|f(x)|}{|W|} ]如果(\epsilon 10^{-6})说明这个配平点已经足够精确如果只到(10^{-3})建议缩小StepTolerance或更换初值。数值噪声方面气动插值表如果做得很粗糙导数不连续会让优化器很难判断梯度方向。解决方案是使用二阶光滑的插值方法或者在MATLAB里用griddedInterpolant设spline方法。当然spline也可能带来过冲对于气动数据我一般用pchip既保证光滑又不会产生明显假的震荡。5. 从配平到动态仿真验证把脚本结果接回仿真模型5.1 配平输出如何初始化Simulink模型算出配平点只是第一步更关键的是让仿真程序在这个配平点上跑起来。在Simulink里你可能需要把配平得到的迎角、俯仰角、升降舵偏度和油门值换算成初始状态向量和控制输入向量。举个例子如果模型的状态向量是([V, \alpha, \theta, q, x, z])那么平飞配平条件下(V V_{\text{target}})(\alpha \alpha_{\text{trim}})(\theta \alpha_{\text{trim}})因为爬升角为0(q 0)(x, z) 任意控制输入(\delta_e \delta_{e_{\text{trim}}})(\delta_T \delta_{T_{\text{trim}}})把这些值写到工作区Simulink模型用对应的变量名初始化积分器和输入端口即可。最关键的是仿真开始后第一步状态导数应该接近零。5.2 验证配平点是否真的稳定扰动测试配平点只代表“力与力矩平衡”不代表“稳定”稳定性是另一个问题。为了验证配平点的可用性我会做两类扰动测试迎角扰动给一个幅度为(1^\circ)~(3^\circ)的瞬时迎角偏差看飞机能否在几秒内回到配平状态。若发散说明静稳定性不足或配平点在失速边界附近。升降舵阶跃给一个很小的舵面阶跃扰动看俯仰速率是否收敛。若出现持续振荡可能配平点附近阻尼不足。这个步骤能帮助你判断控制律设计的起点。很多情况下配平点算得准但动态响应很糟糕这不属于配平脚本的锅而是需要加入增稳控制。5.3 从平飞配平到爬升、转弯配平平飞配平是最基础的情况。实际飞行任务里定常爬升和定常转弯配平同样重要。定常爬升时爬升角(\gamma \neq 0)力的方程中重力分量不再只沿法向需要额外处理。此时迎角不再是(\theta)而是(\theta - \gamma)。你可以用同样的残差脚本传入不同的(\gamma)值批量扫掠。定常转弯时需要考虑侧向力平衡和偏航力矩平衡未知量会扩展到副翼、方向舵和滚转角。如果你只有一个纵向配平脚本可以把它扩展为完整六自由度配平或者先用平飞配平结果作为初值再放开副翼和方向舵。我推荐后者因为交叉耦合比较强直接全自由度高维求解容易发散。下面是一个简单的转向角扫描示例思路gamma_list deg2rad(0:5:15); trim_points zeros(length(gamma_list), 3); for i 1:length(gamma_list) flight_condition.gamma gamma_list(i); [res, x_trim] solve_trim(aircraft, flight_condition, aero, x0); trim_points(i, :) x_trim; end plot(gamma_list, rad2deg(trim_points(:,1)), -o);这种扫描在项目初期非常有用能快速看到飞机的配平范围判断每个飞行状态是否可达。6. 个人工具箱与脚本的实用建议6.1 把工具箱函数和手写脚本分开管理很多初学者喜欢把所有计算都塞到Simulink模型里觉得可视化方便。但配平计算更适合做成独立脚本因为要反复迭代、调参、批量扫描用脚本比在模型里点界面高效得多。我的做法是建一个trim_utils目录里面放几个固定接口的函数build_aircraft.m定义飞机参数load_aero_tables.m加载气动数据trim_solve.m调用优化工具箱求解trim_validate.m验证残差和边界update_dynamic_model.m把配平结果写回仿真模型。这样每个函数职责单一换飞机、改飞行状态都很方便。配合脚本整个流程可以一键重跑这也是“工具箱脚本”最大的优势可复现、可追溯。检查积分器初值、插值表是否越界时如果都是脚本生成排查会快很多。6.2 提高成功率的几个工具链优化技巧我在实践中常用的几个小技巧在这里一并分享给气动表加上范围检查越界时直接报错而不是外推。静稳定配平点如果出现在插值表边界上多半是数据范围不够需要补充气动数据。使用多个初值批量尝试。写一个循环从初值网格里取多个点做配平最后选择残差最小且不落在边界上的那个解。这比单点求解稳很多。打开优化求解器迭代输出但不要迷信它。看到“iteration”的残差变化可以帮判断是否卡在边界上。把气动系数改成无量纲形式配平计算中保持一致的单位体系否则很容易出现量级差大导致优化失败的尴尬。6.3 后续扩展方向从配平到飞行包线分析配平点一旦可批量计算你就能大大方方地做飞行包线扫描。速度从失速速度到最大速度高度从海平面到升限爬升角从负到正每个状态点都算出一组配平值就可以绘制平飞包线。遇到包线边缘点不收敛回头检查气动数据和控制权限往往能发现模型层面的问题。我后来在几个项目里把配平脚本和后续的传感器仿真、控制律设计串成了一条流水线先配平再做线性化再设计控制器最后投入非线性仿真。每一步都基于同一个模型基准误差来源就很清晰。这个流程看起来朴素但比直接上复杂的自动调参工具要可靠得多。配平这件事说到底是飞行仿真里绕不开的底层能力。工具箱负责把数值方法跑稳脚本负责把流程固化下来两者结合才能让仿真程序从一个“能开机”的演示变成一个能真正支撑飞行控制设计的工具。碰到配平发散时别急着怀疑求解器先从模型定义、坐标系、气动数据范围这些基础查起大部分问题都能找到答案。
返回列表