简介一套面向航空路径规划与燃油效率优化的开源代码包聚焦亚特兰大至西雅图的最佳航线求解。项目利用历史航路点数据结合最短路径算法模拟燃料消耗在MATLAB与Python环境下完成天气数据获取、可行网络构建与候选路径评估适合算法研究人员、交通运输专业学生及开源爱好者参考。包内共7个文件2个MATLAB脚本用于航路点抓取和天气数据处理2个Excel表格存放样本航班与气象数据1个Python程序通过二分法构建航路网络另有1个Markdown说明文档和1个PDF论文阐述建模过程整个压缩包仅949KB结构简洁。目前已有461人学习浏览。借助该资源可掌握从FlightAware采集航路点、关联天气条件到运行最短路径算法的完整实践链并依据论文思路复现实验、迁移到其他起讫城市是理解燃油模型与飞行路径优化的不错入门案例。1. 为什么飞行路径优化往往死在燃油模型上先算准油才知道往哪飞做飞行路径优化的人常有一个错觉把搜索算法调得足够精致航迹就能省油。我见过不少照着公开代码跑出来的结果曲线很光滑一细算油耗比直飞还高。问题几乎全出在燃油模型而不是优化算法。飞行路径优化本质上是个最优控制问题目标函数里那点燃油必须由模型算出来模型差一点优化结果就偏一截。这篇文章就按我在类似项目里的习惯把从燃油模型到 FlightPathOptimization 落地这一路拆开MATLAB 里怎么建模、怎么把航迹变成优化变量、参数怎么定、坑在哪。2. 燃油模型的 MATLAB 落地从 Breguet 方程到瞬时燃油流量函数很多人拿到飞行路径优化任务第一步就去挑优化算法这是常见的次序错误。燃油模型才是真正决定优化结果偏不偏的东西同样一条航迹阻力系数差 5%最优高度能差出上千米TSFC 差 10%最优速度直接换一个档位。所以先把燃油模型在 MATLAB 里写稳优化才有意义。这一章就是把这个地基搭起来。我会从三层常见模型里选一个适合优化的方案然后把代码写到能直接跑的程度最后补上标定和验证的环节。2.1 三种燃油模型怎么选路径优化为什么不直接用 Breguet 方程燃油模型在工程上大概分三层每一层都有它的适用位置。第一层是 Breguet 航程方程。它在巡航状态下给出一个漂亮的解析关系飞行航程 (真空速 / (TSFC × g)) × (升阻比) × ln(初始质量 / 最终质量)反过来可以估算航程对应的燃油质量。这个公式用于快速估算一架飞机从 A 到 B 总共烧多少油非常方便但它假设全程定常巡航、升重平衡、TSFC 恒定。实际航迹里爬升、下降、风场变化都会打破这些假设。如果把 Breguet 方程直接塞进路径优化模型会天然偏向“把整段路都当巡航飞”爬升段和下降段直接被忽略优化结果没法看。第二层是准静态燃油流量模型。它不直接算总油耗而是把每个飞行瞬时的需用推力算出来再按发动机的推力燃油消耗率换算成瞬时燃油流量。路径优化时航迹被切成一小段一小段每一段都有不同的高度、速度、重量和推力把这些瞬时流量沿航迹积分就得到总油耗。这一层模型抓住了“高度越高阻力越小但发动机效率也会变、速度越快阻力越大但飞行时间越短”这些核心物理趋势计算量又小最适合作为飞行路径优化的目标函数。第三层是发动机性能表模型直接查厂家给的推力、燃油流量数据表精度最高但数据通常拿不全而且插值表放进优化循环会带来非光滑问题收敛很容易翻车。除非你有真实发动机数据源否则不推荐在路径优化项目里一上来就上这层。我的习惯是优化阶段用第二层等优化出一个候选航迹再用第三层模型做一次高保真复核。这样既保证优化能收敛又避免模型简化带来的系统偏差。2.2 瞬时燃油流量函数的三件套标准大气、升阻力、TSFC准静态模型要落地需要三个 MATLAB 函数标准大气、升阻力计算、燃油流量计算。标准大气函数负责把高度换算成密度、温度、压力这是计算升力和阻力的基础。下面这版覆盖 0 到 20 公里对民航和大型无人机都够用。function [rho, p, T, a] stdAtmo(h) % STDATMO 国际标准大气模型 % 输入 h: 高度(m)支持向量 % 输出 rho: 密度(kg/m^3)p: 静压(Pa)T: 温度(K)a: 音速(m/s) g 9.80665; % 重力加速度 m/s^2 R 287.0528; % 空气专用气体常数 J/(kg*K) T0 288.15; % 海平面温度 K p0 101325; % 海平面静压 Pa hb 11000; % 对流层顶高度 m L 0.0065; % 对流层温度递减率 K/m Tb T0 - L*hb; % 对流层顶温度 K pb p0 * (1 - L*hb/T0)^(g/(L*R)); % 对流层顶静压 Pa beta h hb; % 对流层内为 1平流层为 0 T beta .* (T0 - L*h) ~beta .* Tb; p beta .* (p0 .* (1 - L*h/T0).^(g/(L*R))) ... ~beta .* (pb .* exp(-(h-hb)*g/(R*Tb))); rho p ./ (R*T); a sqrt(1.4*R*T); end逻辑说明高度小于 11 公里时用对流层的线性温度和压力公式高于 11 公里后大气等温压力按指数衰减。如果你只需要做低空无人机优化可以把 hb 那一段去掉但保留它能让函数兼容到两万米省得以后换机型时再写一遍。参数说明L 0.0065 K/m 是标准对流层递减率这个值不要随便改它是国际标准大气的定义值R 用的是干空气常数忽略湿度影响对燃油优化这个精度足够。升阻力函数负责把当前飞行状态换算成需用推力。这里采用准平衡假设认为每一小段内法向力基本平衡升力近似等于重量的法向分量。function [T_req, D] requiredThrust(h, Vt, m, gamma) % REQUIREDTHRUST 准平衡飞行需用推力 % h: 高度(m) % Vt: 真空速(m/s) % m: 当前飞机总质量(kg) % gamma: 航迹角(rad)爬升为正 [rho, ~, ~, ~] stdAtmo(h); S 122.6; % 机翼参考面积 m^2按机型修改 g 9.80665; CL 2*m*g*cos(gamma) / (rho * Vt^2 * S); CD0 0.0160; % 零升阻力系数需标定 K 0.042; % 诱导阻力因子需标定 CD CD0 K * CL^2; D 0.5 * rho * Vt^2 * S * CD; T_req D m*g*sin(gamma); % 水平加速项忽略准平衡近似 end逻辑说明阻力分成零升阻力和诱导阻力两部分诱导阻力用抛物线极曲线近似。典型的喷气式飞机巡航升阻比大概在 15 到 20 之间这段代码里 CL 算出来会落在正常范围。参数说明S、CD0、K 这三个值必须按实际机型改。我一般先把 S 从公开几何数据里拿准CD0 和 K 留给后面的标定环节反推不建议从网上随便抄一组因为不同飞机阻力特性差很多。有了推力和大气环境燃油流量就简单了瞬时燃油质量流量等于 TSFC 乘以需用推力。function mf fuelFlow(h, Vt, m, gamma) % FUELFLOW 瞬时燃油质量流量(kg/s) % 内部先算需用推力再乘推力燃油消耗率 [T_req, ~] requiredThrust(h, Vt, m, gamma); CT0 1.5e-5; % 典型涡扇发动机巡航TSFC单位 kg/(N*s) CT CT0 * (1 0.05 * (Vt - 250)/250); % 一阶速度修正 mf CT * T_req; end逻辑说明TSFC 不是常数速度提高时发动机推进效率会变化这里给一个简单的一阶修正以 250 m/s 为参考点。如果手头有真实发动机的燃油流量表这一步应该替换成查表函数。参数说明CT0 的量级要特别注意。1.5e-5 kg/(N·s) 大约相当于每小时每牛推力消耗 54 克燃油这是现代涡扇发动机巡航状态的典型量级。老发动机或者涡桨发动机会明显偏高高压高速状态下也可能不同标定阶段要再校一遍。2.3 标定与验证拿一段已知航段反推模型系数模型写好了不能直接用。某模拟项目 X 里A 同学用我从网上找来的阻力系数跑优化结果最优巡航高度比真实手册值低了 2500 米原因就是 CD0 抄的是另一款机型的。燃油模型必须在用之前做一次标定。标定思路很简单找一段已知的真实飞行数据包括巡航速度、巡航高度、起飞质量、落地质量、总燃油消耗然后调整模型里的 CD0、K、CT0让模型算出来的总油耗逼近真实值。% 标定示例已知某航段巡航数据 V_cruise 248; % m/s h_cruise 11000; % m m0 55000; % 初始质量 kg mf_total_obs 4100; % 实测总油耗 kg % 待反推参数: x [CD0, K, CT0] costFun (x) abs(predictFuel(x) - mf_total_obs); x0 [0.016, 0.042, 1.5e-5]; x_est fminsearch(costFun, x0);逻辑说明predictFuel 就是 2.2 节三个函数联起来从巡航起点积分到终点。fminsearch 是无导数优化对这种三个参数的小问题足够稳。参数说明CD0 和 K 存在强耦合如果观测数据只有巡航段建议把 K 固定在经验范围只反推 CD0 和 CT0否则会得到一组物理上没意义的参数组合。验证标准我习惯定在总燃油误差 ±5% 以内。如果对不上先检查单位再检查高度是不是用英尺写进了公式最后才怀疑模型结构。最常见的不是公式错而是某个输入忘了从英制转成国际单位。3. 把 FlightPathOptimization 变成 MATLAB 能解的数值问题路径参数化、目标函数与求解器燃油模型就位之后接下来要回答一个关键问题一条连续的飞行路径怎么表示成优化算法能处理的变量连续航迹在数学上是无穷维的优化算法没办法直接处理无穷维向量所以必须做参数化。这一步做得好优化收敛快、结果可飞做得糙后面全是坑。3.1 路径参数化一条航迹怎么压缩成几十个变量常见的做法是把航迹分成两个剖面高度剖面和速度剖面。水平航迹通常事先给定比如某航路或某条固定走廊优化的重点放在每一段飞多高、飞多快。我用控制点参数化沿航迹取 N 个位置每个位置给一个高度值和一个真空速值航迹上的连续剖面用插值生成。变量向量就是这些控制点拼在一起% 航迹参数化示例 N 8; % 控制点个数越大航迹越灵活但越难收敛 s_total 3000e3; % 航段总距离 m s_node linspace(0, s_total, N); % 优化变量 x [h_1...h_N; V_1...V_N] x0 [10500 * ones(N,1); 248 * ones(N,1)]; lb [ 7000 * ones(N,1); 200 * ones(N,1)]; ub [12500 * ones(N,1); 280 * ones(N,1)];逻辑说明控制点之间用插值连起来目标函数每算一步就对整条航迹做一次积分。控制点数量是关键参数N 太小时航迹没有自由度最优高度可能根本不在可行范围内N 太大时变量多、目标函数曲面复杂SQP 容易卡在局部解。我一般先用 N6 到 8 跑通看最优剖面形态再逐步加密。参数说明上下界要尽量贴近机型包线。7000 米下限保证有足够的爬升空间12500 米接近典型喷气客机的升限速度下界 200 m/s 对应巡航马赫数约 0.65上界 280 m/s 约马赫数 0.85具体按机型包线改。3.2 目标函数与非线性约束一段能直接跑的 MATLAB 骨架目标函数就是把燃油流量沿航迹积分。注意这里要用距离做自变量而不是时间这样积分步长不随速度变化数值上更稳。function J totalFuel(x, ac, s_node) % TOTALFUEL 飞行路径优化目标函数总燃油消耗(kg) % x: 优化变量前半段是高度控制点后半段是速度控制点 % ac: 飞机参数结构体 % s_node: 控制点对应的水平距离位置(m) N numel(x) / 2; h_node x(1:N); V_node x(N1:end); h_fun (s) interp1(s_node, h_node, s, pchip); V_fun (s) interp1(s_node, V_node, s, pchip); s_seg linspace(0, s_node(end), 500); % 积分点 W ac.Mass0; fuel 0; for i 1:numel(s_seg)-1 s_mid (s_seg(i) s_seg(i1)) / 2; ds s_seg(i1) - s_seg(i); h h_fun(s_mid); Vt V_fun(s_mid); % 由高度剖面求航迹角 h_next h_fun(s_mid 1e-3); gamma atan((h_next - h) / 1e-3); % 计算燃油流量 [rho, ~, ~, ~] stdAtmo(h); CL 2 * W * ac.g * cos(gamma) / (rho * Vt^2 * ac.S); CD ac.CD0 ac.K * CL^2; D 0.5 * rho * Vt^2 * ac.S * CD; T D W * ac.g * sin(gamma); mf ac.CT0 * T; dt ds / Vt; % 过该小段所需时间 fuel fuel mf * dt; W W - mf * dt; % 重量随燃耗更新 end J fuel; end逻辑说明这个函数每被调用一次就算完一整条航迹的燃油。航迹角用高度剖面的局部斜率近似所以高度剖面必须是连续的用 pchip 而不是 linear 插值就是为了避免折角带来不真实的爬升率跳变。重量更新放在积分循环里这是后面第 4 章要重点讲的耦合项不能省。参数说明积分点 500 个是我调试下来精度和速度的折中如果你发现燃油量随积分点数还在明显变化就往 800 或 1000 加。约束函数负责保证航迹可飞。最少需要三类约束爬升率上下界、速度包线、高度范围。function [c, ceq] nonlcon(x, ac, s_node) % NONLCON 非线性约束可飞性约束 N numel(x) / 2; h_node x(1:N); V_node x(N1:end); h_fun (s) interp1(s_node, h_node, s, pchip); c []; % 1) 爬升率约束|gamma| 0.06 rad约 3.4 度 s_test linspace(0, s_node(end), 200); for i 1:numel(s_test) s0 s_test(i); h0 h_fun(s0); h1 h_fun(min(s0 50, s_node(end))); gamma atan((h1 - h0) / 50); c(end1) abs(gamma) - 0.06; %#ok end % 2) 速度包线约束由边界 lb/ub 已保证这里可补失速约束 ceq []; end逻辑说明爬升率约束如果不加优化器可能造出一条瞬间爬升 3000 米的垂直剖面。0.06 rad 对应约每分钟爬升 900 米对大型客机偏保守但作为约束界可以接受具体值按机型爬升性能改。约束函数里用循环是直观写法性能要求高时可以向量化后面避坑章会讲。3.3 求解器怎么选fmincon 是主力遗传算法只能当替补参数化完成之后优化问题变成一个标准的带约束非线性规划MATLAB 里最顺手的是 fmincon。遗传算法 ga 也能用但它把连续变量当离散基因搜索收敛慢而且没法精细满足非线性约束我一般只在问题高度非凸、fmincon 多起点全部失败时才拿它兜底。求解器优点缺点适合场景fmincon SQP收敛快约束处理强对初值敏感连续航迹优化默认首选fmincon interior-point大规模稳定性好迭代慢变量数超过 50 时ga全局搜索不易陷局部解慢约束弱离散航路点选择、绕飞规划动态规划可处理随机风场维数灾难高度层离散选择问题主程序调用方式options optimoptions(fmincon, ... Algorithm, sqp, ... MaxIterations, 300, ... Display, iter, ... OptimalityTolerance, 1e-6); [x_opt, J_opt, exitflag] fmincon((x) totalFuel(x, ac, s_node), ... x0, [], [], [], [], lb, ub, ... (x) nonlcon(x, ac, s_node), options);逻辑说明SQP 算法是解决平滑非线性约束问题的默认选择收敛快对中等规模变量最稳。OptimalityTolerance 不要设得太小1e-6 已经足够设到 1e-10 只会让求解器在数值噪声里空转。参数说明exitflag 必须检查等于 1 或 2 才算正常收敛如果等于 0说明达到迭代上限结果不能用等于 -2 说明问题不可行先检查约束是不是互相矛盾。4. 让优化结果真的能飞风场、重量耦合和航迹角约束fmincon 跑通一轮之后你大概率会得到一个“数学上最优”的剖面。但飞行路径优化最要命的一步来了这个剖面搬到真实环境里可不可飞、省不省油取决于你有没有把那几个容易忽略的物理项装进模型。这一章讲三件事风场怎么进模型、重量耦合为什么不能省、航迹角约束的物理来源。4.1 风场怎么进模型顺风逆风对最优高度的影响如果不考虑风优化结果是高度越高越好——高空空气稀薄阻力小同样的速度只需要更小推力。但真实航线上高空往往有强风顺风段应该往高处飞去吃风逆风段则要压低高度避开强风。这个权衡不建模优化结果在风大的航段会偏得离谱。风进模型的正确位置是时间计算。燃油流量从真空速和需用推力算但对地速度才决定飞过一段距离花多少时间。计算时间时要用地速% 沿航迹方向风对地速的影响 windAlong interp1(wind_s, wind_along, s_mid, linear); GS Vt windAlong; % 顺风为正逆风为负 dt ds / GS; % 用对地速度算时间 fuel fuel mf * dt;逻辑说明风场数据按距离给出插值到积分点。顺风时 GS 大于 Vt飞同样距离时间更短燃耗自然降低逆风则相反。优化器会自动把顺风段的高度拉高、逆风段压低这才是有实际价值的飞行路径优化。参数说明这里把风简化为只考虑沿航迹方向的分量。实际航路上还存在侧风侧风会导致偏流需要同时优化航向角属进阶内容放到第 6 章。4.2 重量耦合不能省为什么必须一边积分一边减油很多初学者会把飞机质量当成常数放进目标函数理由是燃油质量占总体比例小。这个近似在短航段勉强能用但在两千公里以上的航段会明显失真。一架起飞质量 55 吨的飞机飞完三千公里烧掉五六吨油质量变化接近 10%。质量变小后需要的升力变小诱导阻力下降最优巡航高度会逐渐升高。忽视这个耦合优化器给出的“最优高度”其实是整条航迹前半段的平均而不是每一段真正的最优高度。解决办法已经在 3.2 的代码里写进去了积分循环里每走一小段就执行一次 W W - mf * dt。这里再强调一下它的效果随着 W 变小需用升力系数 CL 下降诱导阻力减小燃油流量随之下降于是整条曲线的油耗不是线性的而是前高后低。如果你把质量固定误差方向一定是低估后半段油耗高估前半段最终总油耗可能偏差 3% 到 8%。这个量级在优化里足以改变最优高度剖面。4.3 航迹角约束的物理来源爬升率是用燃油换来的3.2 的约束函数里设了航迹角绝对值不超过 0.06 弧度这个值不是随便拍的。航迹角的正弦值乘以真空速就是爬升率0.06 弧度在 250 m/s 时对应约每秒 15 米折合每分钟 900 米。对喷气式客机来说这是一个合理的爬升率上界超过这个值机翼要产生更大的诱导阻力而且发动机推力可能不够维持。更重要的是航迹角直接进入需用推力公式T D m g sin γ。爬升时推力要多出一部分去对抗重力分量这部分推力同样烧油。如果约束函数不加优化器会试着用剧烈的锯齿剖面去“投机”——下降段省油爬升段补回高度整体看起来油耗降低实际飞起来完全不可行。约束的物理来源决定了它该怎么取值上限由最大爬升性能决定下限由下降时的飞行品质决定下降率太大会让旅客觉得失重客机一般限制在每分钟 300 到 500 米以下。这个约束宁可保守也不要放宽因为它直接决定优化结果能不能被真实飞行手册接受。5. 飞行路径优化避坑指南MATLAB 实现里的 5 个高频问题这一章我把实际调试中重复踩过的坑按“现象 → 原因 → 解决”写出来。每一条都对应真实翻车经历照方抓药能省出一两天排查时间。5.1 现象fmincon 迭代几步就停说目标函数不再下降原因目标函数值本身量级太大。总燃油动不动几万千克fmincon 默认的步长和容差是基于相对值的它会在一个看起来“已经收敛”的位置提前停下实际离最优还远。有时原因是数值梯度误差目标函数里积分步长不够细导致相邻两点目标值差异被噪声淹没。解决把目标函数缩放一个量级。常见做法是让 J 返回以吨为单位的燃油量比如 J fuel / 1000然后同步调小 FinDiffRelStep 来获得更稳定的有限差分梯度。options optimoptions(fmincon, ... Algorithm, sqp, ... ScaleProblem, true, ... FinDiffRelStep, 1e-7);逻辑说明ScaleProblem 让求解器自动对变量和目标做缩放FinDiffRelStep 控制数值微分步长1e-7 在高度和速度量纲下是一个安全的初始值。如果优化还是早停下一步检查积分点数是否足够500 点对三千公里航段通常足够但如果控制点很密剖面变化剧烈要加积分点。5.2 现象优化出的高度剖面出现锯齿上上下下完全不可飞原因控制点之间用 linear 插值时折点处的航迹角发生跳变目标函数在这个位置不光滑数值梯度误导求解器。另一个常见原因是控制点数量过多自由度太大模型开始拟合数值噪声。解决改用 pchip 插值并用正则项把相邻控制点之间的变化压平。正则的本质是给目标函数加一点“高度变化成本”让优化器只有在真正省油时才会去改变高度。% 对高度控制点加平滑正则 lambda 5e3; smoothPenalty lambda * sum(diff(h_node).^2); J J smoothPenalty;逻辑说明lambda 是经验参数太小没作用太大把真实的高度优化也压掉了。我一般从 1e3 试到 1e4观察最优剖面从“锯齿”变成“平滑”的拐点取刚好消除锯齿的最小值。如果你发现减少控制点后锯齿自然消失那就先减控制点正则项留到后面处理更复杂的约束。5.3 现象换一个初始猜测最优解完全不同原因目标函数是非凸的存在多个局部最优。不同初值落在不同吸引域里SQP 只能保证找到局部解。高空强风、速度包线约束都会加剧非凸性。解决多起点策略这是最可靠也最土的办法。在 lb 和 ub 范围内随机生成 20 组初值分别跑 fmincon取油耗最小的结果。rng(42); bestJ inf; for k 1:20 xTry lb rand(numel(lb),1) .* (ub - lb); [xTmp, JTmp] fmincon((x) totalFuel(x, ac, s_node), ... xTry, [], [], [], [], lb, ub, ... (x) nonlcon(x, ac, s_node), options); if JTmp bestJ bestX xTmp; bestJ JTmp; end end逻辑说明随机初值本身不是科学但它能暴露出目标函数的非凸程度。如果 20 个起点得到的结果分散说明模型里物理约束不够强优先检查速度包线和爬升率约束是否太松。参数说明81 行里的 rng 固定种子很重要否则每次结果不可复现后面调参没法对比。5.4 现象约束函数的单位混用导致路径不可飞现象爬升率约束写成 |gamma| 0.06但高度控制点用的是英尺、距离控制点用的是海里最后算出来 gamma 放大了 1852 倍约束永远违反fmincon 一开始就报无解。原因路径参数化时长度单位没统一。MATLAB 里最隐蔽的单位坑是把速度写成节、高度写成英尺、距离写成海里而标准大气函数期待的是米和秒。解决在入口处加断言强制单位校验。这个习惯能省掉一整类问题validateattributes(h_node, {double}, {column, , 0, , 20000}); validateattributes(V_node, {double}, {column, , 100, , 350});逻辑说明validateattributes 在变量不合预期时直接抛错把问题暴露在程序入口而不是让一个奇怪的优化结果兜圈子。参数说明高度上界 20000 对应标准大气函数适用范围速度上下界按飞机包线写。单位制建议全工程统一用国际单位制只有输入数据和输出报告才做转换。5.5 现象模型算出的总油耗比真实手册高 20%原因最常出现在 CT0 和 K 两个参数上。CT0 如果按最大推力状态取值巡航油耗会被系统性高估K 如果偏大高空巡航时的诱导阻力会被高估导致最优高度偏保守。第二种常见原因是升力公式里忘了乘 cos(gamma)爬升段推力被低估油耗偏低——两种情况方向相反所以光看总油耗数字对不上猜不出是哪个参数的问题。解决不要猜分段验证。先用平飞巡航段单独验证阻力极曲线输入真实巡航速度、高度、质量算出来的平飞需用推力应当和发动机厂家给出的巡航推力接近。推力对上了再验证燃油流量用真实巡航燃油流量反推 CT0。血泪经验是CT0 和 K 的偏差对最优高度的作用方向相反同时调两个参数可以让总油耗误差看起来很小但最优剖面是错的。所以标定必须一个一个来先定阻力再定燃油流量不要两参数一起最小二乘。6. 进阶从高度剖面到三维航迹以及结果验证的两个习惯6.1 三维航迹的最小扩展把航向角加进变量二维剖面只优化每段飞多高、飞多快航迹在地图上是固定的。三维扩展要把航向角也变成优化变量让飞机在侧风环境下通过调整航向抵消偏流从而减少实际飞过的地面距离。最小改动是在变量向量里追加一组航向控制点% 三维扩展变量: x [h_node; V_node; psi_node] x0 [10500*ones(N,1); 248*ones(N,1); 0.0*ones(N,1)];目标函数里把风速分解成沿航迹和垂直航迹两个分量侧风分量通过航向角修正影响对地速度这一步引入之后优化器会自动让飞机在侧风区域调整航向。要注意的是航向角变量会引入额外的非凸性多起点策略从可选变成必选。6.2 验证最优结果的能量视角拿到优化结果先别急着信。第一个快速检查是能量守恒整段燃耗对应的能量必须大于飞机势能和动能增量加上全程阻力耗能。如果总燃油对应的能量比物理下限还低说明积分有 bug 或者单位错乱。第二个检查是看巡航段的最优真空速是否接近理论最佳巡航速度% 理论最佳巡航真空速近似 V_br sqrt(2 * ac.K * W * ac.g / (rho * ac.S * ac.CD0));逻辑说明在定常平飞且忽略 TSFC 随速度变化的假设下这个公式给出最小阻力对应的速度。真实最优速度会因发动机效率偏移但趋势应该一致高空密度小V_br 会变大。如果优化出的巡航速度明显偏离这个值先查 TSFC 的速度修正是不是方向写反了。6.3 一个改不掉的习惯参数全部收敛到结构体里我早期做模拟项目 X 时CD0、K、CT0 散落在三个脚本里某次改了一个脚本里的 CT0另一个脚本没同步优化结果全部偏移查了整整一天。后来所有模型参数一律收进 ac 结构体脚本里只传结构体不再出现裸数字。这个习惯没什么技术含量但它是飞行路径优化项目里投入产出比最高的一个决定。优化问题本身不难难的是让每个参数都有明确出处、每个约束都有物理依据、每个结果都能复现。先算准油再谈优化希望帮到你。本文还有配套的精品资源点击获取
简介一套面向航空路径规划与燃油效率优化的开源代码包聚焦亚特兰大至西雅图的最佳航线求解。项目利用历史航路点数据结合最短路径算法模拟燃料消耗在MATLAB与Python环境下完成天气数据获取、可行网络构建与候选路径评估适合算法研究人员、交通运输专业学生及开源爱好者参考。包内共7个文件2个MATLAB脚本用于航路点抓取和天气数据处理2个Excel表格存放样本航班与气象数据1个Python程序通过二分法构建航路网络另有1个Markdown说明文档和1个PDF论文阐述建模过程整个压缩包仅949KB结构简洁。目前已有461人学习浏览。借助该资源可掌握从FlightAware采集航路点、关联天气条件到运行最短路径算法的完整实践链并依据论文思路复现实验、迁移到其他起讫城市是理解燃油模型与飞行路径优化的不错入门案例。1. 为什么飞行路径优化往往死在燃油模型上先算准油才知道往哪飞做飞行路径优化的人常有一个错觉把搜索算法调得足够精致航迹就能省油。我见过不少照着公开代码跑出来的结果曲线很光滑一细算油耗比直飞还高。问题几乎全出在燃油模型而不是优化算法。飞行路径优化本质上是个最优控制问题目标函数里那点燃油必须由模型算出来模型差一点优化结果就偏一截。这篇文章就按我在类似项目里的习惯把从燃油模型到 FlightPathOptimization 落地这一路拆开MATLAB 里怎么建模、怎么把航迹变成优化变量、参数怎么定、坑在哪。2. 燃油模型的 MATLAB 落地从 Breguet 方程到瞬时燃油流量函数很多人拿到飞行路径优化任务第一步就去挑优化算法这是常见的次序错误。燃油模型才是真正决定优化结果偏不偏的东西同样一条航迹阻力系数差 5%最优高度能差出上千米TSFC 差 10%最优速度直接换一个档位。所以先把燃油模型在 MATLAB 里写稳优化才有意义。这一章就是把这个地基搭起来。我会从三层常见模型里选一个适合优化的方案然后把代码写到能直接跑的程度最后补上标定和验证的环节。2.1 三种燃油模型怎么选路径优化为什么不直接用 Breguet 方程燃油模型在工程上大概分三层每一层都有它的适用位置。第一层是 Breguet 航程方程。它在巡航状态下给出一个漂亮的解析关系飞行航程 (真空速 / (TSFC × g)) × (升阻比) × ln(初始质量 / 最终质量)反过来可以估算航程对应的燃油质量。这个公式用于快速估算一架飞机从 A 到 B 总共烧多少油非常方便但它假设全程定常巡航、升重平衡、TSFC 恒定。实际航迹里爬升、下降、风场变化都会打破这些假设。如果把 Breguet 方程直接塞进路径优化模型会天然偏向“把整段路都当巡航飞”爬升段和下降段直接被忽略优化结果没法看。第二层是准静态燃油流量模型。它不直接算总油耗而是把每个飞行瞬时的需用推力算出来再按发动机的推力燃油消耗率换算成瞬时燃油流量。路径优化时航迹被切成一小段一小段每一段都有不同的高度、速度、重量和推力把这些瞬时流量沿航迹积分就得到总油耗。这一层模型抓住了“高度越高阻力越小但发动机效率也会变、速度越快阻力越大但飞行时间越短”这些核心物理趋势计算量又小最适合作为飞行路径优化的目标函数。第三层是发动机性能表模型直接查厂家给的推力、燃油流量数据表精度最高但数据通常拿不全而且插值表放进优化循环会带来非光滑问题收敛很容易翻车。除非你有真实发动机数据源否则不推荐在路径优化项目里一上来就上这层。我的习惯是优化阶段用第二层等优化出一个候选航迹再用第三层模型做一次高保真复核。这样既保证优化能收敛又避免模型简化带来的系统偏差。2.2 瞬时燃油流量函数的三件套标准大气、升阻力、TSFC准静态模型要落地需要三个 MATLAB 函数标准大气、升阻力计算、燃油流量计算。标准大气函数负责把高度换算成密度、温度、压力这是计算升力和阻力的基础。下面这版覆盖 0 到 20 公里对民航和大型无人机都够用。function [rho, p, T, a] stdAtmo(h) % STDATMO 国际标准大气模型 % 输入 h: 高度(m)支持向量 % 输出 rho: 密度(kg/m^3)p: 静压(Pa)T: 温度(K)a: 音速(m/s) g 9.80665; % 重力加速度 m/s^2 R 287.0528; % 空气专用气体常数 J/(kg*K) T0 288.15; % 海平面温度 K p0 101325; % 海平面静压 Pa hb 11000; % 对流层顶高度 m L 0.0065; % 对流层温度递减率 K/m Tb T0 - L*hb; % 对流层顶温度 K pb p0 * (1 - L*hb/T0)^(g/(L*R)); % 对流层顶静压 Pa beta h hb; % 对流层内为 1平流层为 0 T beta .* (T0 - L*h) ~beta .* Tb; p beta .* (p0 .* (1 - L*h/T0).^(g/(L*R))) ... ~beta .* (pb .* exp(-(h-hb)*g/(R*Tb))); rho p ./ (R*T); a sqrt(1.4*R*T); end逻辑说明高度小于 11 公里时用对流层的线性温度和压力公式高于 11 公里后大气等温压力按指数衰减。如果你只需要做低空无人机优化可以把 hb 那一段去掉但保留它能让函数兼容到两万米省得以后换机型时再写一遍。参数说明L 0.0065 K/m 是标准对流层递减率这个值不要随便改它是国际标准大气的定义值R 用的是干空气常数忽略湿度影响对燃油优化这个精度足够。升阻力函数负责把当前飞行状态换算成需用推力。这里采用准平衡假设认为每一小段内法向力基本平衡升力近似等于重量的法向分量。function [T_req, D] requiredThrust(h, Vt, m, gamma) % REQUIREDTHRUST 准平衡飞行需用推力 % h: 高度(m) % Vt: 真空速(m/s) % m: 当前飞机总质量(kg) % gamma: 航迹角(rad)爬升为正 [rho, ~, ~, ~] stdAtmo(h); S 122.6; % 机翼参考面积 m^2按机型修改 g 9.80665; CL 2*m*g*cos(gamma) / (rho * Vt^2 * S); CD0 0.0160; % 零升阻力系数需标定 K 0.042; % 诱导阻力因子需标定 CD CD0 K * CL^2; D 0.5 * rho * Vt^2 * S * CD; T_req D m*g*sin(gamma); % 水平加速项忽略准平衡近似 end逻辑说明阻力分成零升阻力和诱导阻力两部分诱导阻力用抛物线极曲线近似。典型的喷气式飞机巡航升阻比大概在 15 到 20 之间这段代码里 CL 算出来会落在正常范围。参数说明S、CD0、K 这三个值必须按实际机型改。我一般先把 S 从公开几何数据里拿准CD0 和 K 留给后面的标定环节反推不建议从网上随便抄一组因为不同飞机阻力特性差很多。有了推力和大气环境燃油流量就简单了瞬时燃油质量流量等于 TSFC 乘以需用推力。function mf fuelFlow(h, Vt, m, gamma) % FUELFLOW 瞬时燃油质量流量(kg/s) % 内部先算需用推力再乘推力燃油消耗率 [T_req, ~] requiredThrust(h, Vt, m, gamma); CT0 1.5e-5; % 典型涡扇发动机巡航TSFC单位 kg/(N*s) CT CT0 * (1 0.05 * (Vt - 250)/250); % 一阶速度修正 mf CT * T_req; end逻辑说明TSFC 不是常数速度提高时发动机推进效率会变化这里给一个简单的一阶修正以 250 m/s 为参考点。如果手头有真实发动机的燃油流量表这一步应该替换成查表函数。参数说明CT0 的量级要特别注意。1.5e-5 kg/(N·s) 大约相当于每小时每牛推力消耗 54 克燃油这是现代涡扇发动机巡航状态的典型量级。老发动机或者涡桨发动机会明显偏高高压高速状态下也可能不同标定阶段要再校一遍。2.3 标定与验证拿一段已知航段反推模型系数模型写好了不能直接用。某模拟项目 X 里A 同学用我从网上找来的阻力系数跑优化结果最优巡航高度比真实手册值低了 2500 米原因就是 CD0 抄的是另一款机型的。燃油模型必须在用之前做一次标定。标定思路很简单找一段已知的真实飞行数据包括巡航速度、巡航高度、起飞质量、落地质量、总燃油消耗然后调整模型里的 CD0、K、CT0让模型算出来的总油耗逼近真实值。% 标定示例已知某航段巡航数据 V_cruise 248; % m/s h_cruise 11000; % m m0 55000; % 初始质量 kg mf_total_obs 4100; % 实测总油耗 kg % 待反推参数: x [CD0, K, CT0] costFun (x) abs(predictFuel(x) - mf_total_obs); x0 [0.016, 0.042, 1.5e-5]; x_est fminsearch(costFun, x0);逻辑说明predictFuel 就是 2.2 节三个函数联起来从巡航起点积分到终点。fminsearch 是无导数优化对这种三个参数的小问题足够稳。参数说明CD0 和 K 存在强耦合如果观测数据只有巡航段建议把 K 固定在经验范围只反推 CD0 和 CT0否则会得到一组物理上没意义的参数组合。验证标准我习惯定在总燃油误差 ±5% 以内。如果对不上先检查单位再检查高度是不是用英尺写进了公式最后才怀疑模型结构。最常见的不是公式错而是某个输入忘了从英制转成国际单位。3. 把 FlightPathOptimization 变成 MATLAB 能解的数值问题路径参数化、目标函数与求解器燃油模型就位之后接下来要回答一个关键问题一条连续的飞行路径怎么表示成优化算法能处理的变量连续航迹在数学上是无穷维的优化算法没办法直接处理无穷维向量所以必须做参数化。这一步做得好优化收敛快、结果可飞做得糙后面全是坑。3.1 路径参数化一条航迹怎么压缩成几十个变量常见的做法是把航迹分成两个剖面高度剖面和速度剖面。水平航迹通常事先给定比如某航路或某条固定走廊优化的重点放在每一段飞多高、飞多快。我用控制点参数化沿航迹取 N 个位置每个位置给一个高度值和一个真空速值航迹上的连续剖面用插值生成。变量向量就是这些控制点拼在一起% 航迹参数化示例 N 8; % 控制点个数越大航迹越灵活但越难收敛 s_total 3000e3; % 航段总距离 m s_node linspace(0, s_total, N); % 优化变量 x [h_1...h_N; V_1...V_N] x0 [10500 * ones(N,1); 248 * ones(N,1)]; lb [ 7000 * ones(N,1); 200 * ones(N,1)]; ub [12500 * ones(N,1); 280 * ones(N,1)];逻辑说明控制点之间用插值连起来目标函数每算一步就对整条航迹做一次积分。控制点数量是关键参数N 太小时航迹没有自由度最优高度可能根本不在可行范围内N 太大时变量多、目标函数曲面复杂SQP 容易卡在局部解。我一般先用 N6 到 8 跑通看最优剖面形态再逐步加密。参数说明上下界要尽量贴近机型包线。7000 米下限保证有足够的爬升空间12500 米接近典型喷气客机的升限速度下界 200 m/s 对应巡航马赫数约 0.65上界 280 m/s 约马赫数 0.85具体按机型包线改。3.2 目标函数与非线性约束一段能直接跑的 MATLAB 骨架目标函数就是把燃油流量沿航迹积分。注意这里要用距离做自变量而不是时间这样积分步长不随速度变化数值上更稳。function J totalFuel(x, ac, s_node) % TOTALFUEL 飞行路径优化目标函数总燃油消耗(kg) % x: 优化变量前半段是高度控制点后半段是速度控制点 % ac: 飞机参数结构体 % s_node: 控制点对应的水平距离位置(m) N numel(x) / 2; h_node x(1:N); V_node x(N1:end); h_fun (s) interp1(s_node, h_node, s, pchip); V_fun (s) interp1(s_node, V_node, s, pchip); s_seg linspace(0, s_node(end), 500); % 积分点 W ac.Mass0; fuel 0; for i 1:numel(s_seg)-1 s_mid (s_seg(i) s_seg(i1)) / 2; ds s_seg(i1) - s_seg(i); h h_fun(s_mid); Vt V_fun(s_mid); % 由高度剖面求航迹角 h_next h_fun(s_mid 1e-3); gamma atan((h_next - h) / 1e-3); % 计算燃油流量 [rho, ~, ~, ~] stdAtmo(h); CL 2 * W * ac.g * cos(gamma) / (rho * Vt^2 * ac.S); CD ac.CD0 ac.K * CL^2; D 0.5 * rho * Vt^2 * ac.S * CD; T D W * ac.g * sin(gamma); mf ac.CT0 * T; dt ds / Vt; % 过该小段所需时间 fuel fuel mf * dt; W W - mf * dt; % 重量随燃耗更新 end J fuel; end逻辑说明这个函数每被调用一次就算完一整条航迹的燃油。航迹角用高度剖面的局部斜率近似所以高度剖面必须是连续的用 pchip 而不是 linear 插值就是为了避免折角带来不真实的爬升率跳变。重量更新放在积分循环里这是后面第 4 章要重点讲的耦合项不能省。参数说明积分点 500 个是我调试下来精度和速度的折中如果你发现燃油量随积分点数还在明显变化就往 800 或 1000 加。约束函数负责保证航迹可飞。最少需要三类约束爬升率上下界、速度包线、高度范围。function [c, ceq] nonlcon(x, ac, s_node) % NONLCON 非线性约束可飞性约束 N numel(x) / 2; h_node x(1:N); V_node x(N1:end); h_fun (s) interp1(s_node, h_node, s, pchip); c []; % 1) 爬升率约束|gamma| 0.06 rad约 3.4 度 s_test linspace(0, s_node(end), 200); for i 1:numel(s_test) s0 s_test(i); h0 h_fun(s0); h1 h_fun(min(s0 50, s_node(end))); gamma atan((h1 - h0) / 50); c(end1) abs(gamma) - 0.06; %#ok end % 2) 速度包线约束由边界 lb/ub 已保证这里可补失速约束 ceq []; end逻辑说明爬升率约束如果不加优化器可能造出一条瞬间爬升 3000 米的垂直剖面。0.06 rad 对应约每分钟爬升 900 米对大型客机偏保守但作为约束界可以接受具体值按机型爬升性能改。约束函数里用循环是直观写法性能要求高时可以向量化后面避坑章会讲。3.3 求解器怎么选fmincon 是主力遗传算法只能当替补参数化完成之后优化问题变成一个标准的带约束非线性规划MATLAB 里最顺手的是 fmincon。遗传算法 ga 也能用但它把连续变量当离散基因搜索收敛慢而且没法精细满足非线性约束我一般只在问题高度非凸、fmincon 多起点全部失败时才拿它兜底。求解器优点缺点适合场景fmincon SQP收敛快约束处理强对初值敏感连续航迹优化默认首选fmincon interior-point大规模稳定性好迭代慢变量数超过 50 时ga全局搜索不易陷局部解慢约束弱离散航路点选择、绕飞规划动态规划可处理随机风场维数灾难高度层离散选择问题主程序调用方式options optimoptions(fmincon, ... Algorithm, sqp, ... MaxIterations, 300, ... Display, iter, ... OptimalityTolerance, 1e-6); [x_opt, J_opt, exitflag] fmincon((x) totalFuel(x, ac, s_node), ... x0, [], [], [], [], lb, ub, ... (x) nonlcon(x, ac, s_node), options);逻辑说明SQP 算法是解决平滑非线性约束问题的默认选择收敛快对中等规模变量最稳。OptimalityTolerance 不要设得太小1e-6 已经足够设到 1e-10 只会让求解器在数值噪声里空转。参数说明exitflag 必须检查等于 1 或 2 才算正常收敛如果等于 0说明达到迭代上限结果不能用等于 -2 说明问题不可行先检查约束是不是互相矛盾。4. 让优化结果真的能飞风场、重量耦合和航迹角约束fmincon 跑通一轮之后你大概率会得到一个“数学上最优”的剖面。但飞行路径优化最要命的一步来了这个剖面搬到真实环境里可不可飞、省不省油取决于你有没有把那几个容易忽略的物理项装进模型。这一章讲三件事风场怎么进模型、重量耦合为什么不能省、航迹角约束的物理来源。4.1 风场怎么进模型顺风逆风对最优高度的影响如果不考虑风优化结果是高度越高越好——高空空气稀薄阻力小同样的速度只需要更小推力。但真实航线上高空往往有强风顺风段应该往高处飞去吃风逆风段则要压低高度避开强风。这个权衡不建模优化结果在风大的航段会偏得离谱。风进模型的正确位置是时间计算。燃油流量从真空速和需用推力算但对地速度才决定飞过一段距离花多少时间。计算时间时要用地速% 沿航迹方向风对地速的影响 windAlong interp1(wind_s, wind_along, s_mid, linear); GS Vt windAlong; % 顺风为正逆风为负 dt ds / GS; % 用对地速度算时间 fuel fuel mf * dt;逻辑说明风场数据按距离给出插值到积分点。顺风时 GS 大于 Vt飞同样距离时间更短燃耗自然降低逆风则相反。优化器会自动把顺风段的高度拉高、逆风段压低这才是有实际价值的飞行路径优化。参数说明这里把风简化为只考虑沿航迹方向的分量。实际航路上还存在侧风侧风会导致偏流需要同时优化航向角属进阶内容放到第 6 章。4.2 重量耦合不能省为什么必须一边积分一边减油很多初学者会把飞机质量当成常数放进目标函数理由是燃油质量占总体比例小。这个近似在短航段勉强能用但在两千公里以上的航段会明显失真。一架起飞质量 55 吨的飞机飞完三千公里烧掉五六吨油质量变化接近 10%。质量变小后需要的升力变小诱导阻力下降最优巡航高度会逐渐升高。忽视这个耦合优化器给出的“最优高度”其实是整条航迹前半段的平均而不是每一段真正的最优高度。解决办法已经在 3.2 的代码里写进去了积分循环里每走一小段就执行一次 W W - mf * dt。这里再强调一下它的效果随着 W 变小需用升力系数 CL 下降诱导阻力减小燃油流量随之下降于是整条曲线的油耗不是线性的而是前高后低。如果你把质量固定误差方向一定是低估后半段油耗高估前半段最终总油耗可能偏差 3% 到 8%。这个量级在优化里足以改变最优高度剖面。4.3 航迹角约束的物理来源爬升率是用燃油换来的3.2 的约束函数里设了航迹角绝对值不超过 0.06 弧度这个值不是随便拍的。航迹角的正弦值乘以真空速就是爬升率0.06 弧度在 250 m/s 时对应约每秒 15 米折合每分钟 900 米。对喷气式客机来说这是一个合理的爬升率上界超过这个值机翼要产生更大的诱导阻力而且发动机推力可能不够维持。更重要的是航迹角直接进入需用推力公式T D m g sin γ。爬升时推力要多出一部分去对抗重力分量这部分推力同样烧油。如果约束函数不加优化器会试着用剧烈的锯齿剖面去“投机”——下降段省油爬升段补回高度整体看起来油耗降低实际飞起来完全不可行。约束的物理来源决定了它该怎么取值上限由最大爬升性能决定下限由下降时的飞行品质决定下降率太大会让旅客觉得失重客机一般限制在每分钟 300 到 500 米以下。这个约束宁可保守也不要放宽因为它直接决定优化结果能不能被真实飞行手册接受。5. 飞行路径优化避坑指南MATLAB 实现里的 5 个高频问题这一章我把实际调试中重复踩过的坑按“现象 → 原因 → 解决”写出来。每一条都对应真实翻车经历照方抓药能省出一两天排查时间。5.1 现象fmincon 迭代几步就停说目标函数不再下降原因目标函数值本身量级太大。总燃油动不动几万千克fmincon 默认的步长和容差是基于相对值的它会在一个看起来“已经收敛”的位置提前停下实际离最优还远。有时原因是数值梯度误差目标函数里积分步长不够细导致相邻两点目标值差异被噪声淹没。解决把目标函数缩放一个量级。常见做法是让 J 返回以吨为单位的燃油量比如 J fuel / 1000然后同步调小 FinDiffRelStep 来获得更稳定的有限差分梯度。options optimoptions(fmincon, ... Algorithm, sqp, ... ScaleProblem, true, ... FinDiffRelStep, 1e-7);逻辑说明ScaleProblem 让求解器自动对变量和目标做缩放FinDiffRelStep 控制数值微分步长1e-7 在高度和速度量纲下是一个安全的初始值。如果优化还是早停下一步检查积分点数是否足够500 点对三千公里航段通常足够但如果控制点很密剖面变化剧烈要加积分点。5.2 现象优化出的高度剖面出现锯齿上上下下完全不可飞原因控制点之间用 linear 插值时折点处的航迹角发生跳变目标函数在这个位置不光滑数值梯度误导求解器。另一个常见原因是控制点数量过多自由度太大模型开始拟合数值噪声。解决改用 pchip 插值并用正则项把相邻控制点之间的变化压平。正则的本质是给目标函数加一点“高度变化成本”让优化器只有在真正省油时才会去改变高度。% 对高度控制点加平滑正则 lambda 5e3; smoothPenalty lambda * sum(diff(h_node).^2); J J smoothPenalty;逻辑说明lambda 是经验参数太小没作用太大把真实的高度优化也压掉了。我一般从 1e3 试到 1e4观察最优剖面从“锯齿”变成“平滑”的拐点取刚好消除锯齿的最小值。如果你发现减少控制点后锯齿自然消失那就先减控制点正则项留到后面处理更复杂的约束。5.3 现象换一个初始猜测最优解完全不同原因目标函数是非凸的存在多个局部最优。不同初值落在不同吸引域里SQP 只能保证找到局部解。高空强风、速度包线约束都会加剧非凸性。解决多起点策略这是最可靠也最土的办法。在 lb 和 ub 范围内随机生成 20 组初值分别跑 fmincon取油耗最小的结果。rng(42); bestJ inf; for k 1:20 xTry lb rand(numel(lb),1) .* (ub - lb); [xTmp, JTmp] fmincon((x) totalFuel(x, ac, s_node), ... xTry, [], [], [], [], lb, ub, ... (x) nonlcon(x, ac, s_node), options); if JTmp bestJ bestX xTmp; bestJ JTmp; end end逻辑说明随机初值本身不是科学但它能暴露出目标函数的非凸程度。如果 20 个起点得到的结果分散说明模型里物理约束不够强优先检查速度包线和爬升率约束是否太松。参数说明81 行里的 rng 固定种子很重要否则每次结果不可复现后面调参没法对比。5.4 现象约束函数的单位混用导致路径不可飞现象爬升率约束写成 |gamma| 0.06但高度控制点用的是英尺、距离控制点用的是海里最后算出来 gamma 放大了 1852 倍约束永远违反fmincon 一开始就报无解。原因路径参数化时长度单位没统一。MATLAB 里最隐蔽的单位坑是把速度写成节、高度写成英尺、距离写成海里而标准大气函数期待的是米和秒。解决在入口处加断言强制单位校验。这个习惯能省掉一整类问题validateattributes(h_node, {double}, {column, , 0, , 20000}); validateattributes(V_node, {double}, {column, , 100, , 350});逻辑说明validateattributes 在变量不合预期时直接抛错把问题暴露在程序入口而不是让一个奇怪的优化结果兜圈子。参数说明高度上界 20000 对应标准大气函数适用范围速度上下界按飞机包线写。单位制建议全工程统一用国际单位制只有输入数据和输出报告才做转换。5.5 现象模型算出的总油耗比真实手册高 20%原因最常出现在 CT0 和 K 两个参数上。CT0 如果按最大推力状态取值巡航油耗会被系统性高估K 如果偏大高空巡航时的诱导阻力会被高估导致最优高度偏保守。第二种常见原因是升力公式里忘了乘 cos(gamma)爬升段推力被低估油耗偏低——两种情况方向相反所以光看总油耗数字对不上猜不出是哪个参数的问题。解决不要猜分段验证。先用平飞巡航段单独验证阻力极曲线输入真实巡航速度、高度、质量算出来的平飞需用推力应当和发动机厂家给出的巡航推力接近。推力对上了再验证燃油流量用真实巡航燃油流量反推 CT0。血泪经验是CT0 和 K 的偏差对最优高度的作用方向相反同时调两个参数可以让总油耗误差看起来很小但最优剖面是错的。所以标定必须一个一个来先定阻力再定燃油流量不要两参数一起最小二乘。6. 进阶从高度剖面到三维航迹以及结果验证的两个习惯6.1 三维航迹的最小扩展把航向角加进变量二维剖面只优化每段飞多高、飞多快航迹在地图上是固定的。三维扩展要把航向角也变成优化变量让飞机在侧风环境下通过调整航向抵消偏流从而减少实际飞过的地面距离。最小改动是在变量向量里追加一组航向控制点% 三维扩展变量: x [h_node; V_node; psi_node] x0 [10500*ones(N,1); 248*ones(N,1); 0.0*ones(N,1)];目标函数里把风速分解成沿航迹和垂直航迹两个分量侧风分量通过航向角修正影响对地速度这一步引入之后优化器会自动让飞机在侧风区域调整航向。要注意的是航向角变量会引入额外的非凸性多起点策略从可选变成必选。6.2 验证最优结果的能量视角拿到优化结果先别急着信。第一个快速检查是能量守恒整段燃耗对应的能量必须大于飞机势能和动能增量加上全程阻力耗能。如果总燃油对应的能量比物理下限还低说明积分有 bug 或者单位错乱。第二个检查是看巡航段的最优真空速是否接近理论最佳巡航速度% 理论最佳巡航真空速近似 V_br sqrt(2 * ac.K * W * ac.g / (rho * ac.S * ac.CD0));逻辑说明在定常平飞且忽略 TSFC 随速度变化的假设下这个公式给出最小阻力对应的速度。真实最优速度会因发动机效率偏移但趋势应该一致高空密度小V_br 会变大。如果优化出的巡航速度明显偏离这个值先查 TSFC 的速度修正是不是方向写反了。6.3 一个改不掉的习惯参数全部收敛到结构体里我早期做模拟项目 X 时CD0、K、CT0 散落在三个脚本里某次改了一个脚本里的 CT0另一个脚本没同步优化结果全部偏移查了整整一天。后来所有模型参数一律收进 ac 结构体脚本里只传结构体不再出现裸数字。这个习惯没什么技术含量但它是飞行路径优化项目里投入产出比最高的一个决定。优化问题本身不难难的是让每个参数都有明确出处、每个约束都有物理依据、每个结果都能复现。先算准油再谈优化希望帮到你。本文还有配套的精品资源点击获取
