1. 项目概述从赛题到实战的完整复盘去年带队参加“华为杯”数学建模竞赛选的F题“飞行器质心平衡供油策略优化”至今记忆犹新。这道题本质上是一个典型的多约束动态优化问题核心目标是在保证飞行器飞行姿态稳定的前提下通过优化多个油箱之间的燃油调度策略实现飞行航程的最大化或燃油消耗的最优化。听起来像是航空航天领域的专业问题但其内核的优化思想、建模方法和求解技巧在物流调度、生产排程、资源分配等众多工业场景中都有广泛应用。我们当时花了三天三夜最终拿了个不错的奖项过程中踩了不少坑也积累了一些心得。今天就把这道题的解题思路、模型构建、算法实现以及那些“教科书上不会写”的实操细节做一个完整的梳理和分享。无论你是正在备赛的研究生还是对运筹优化、动态规划感兴趣的技术爱好者相信这篇从实战中淬炼出的总结都能给你带来直接的启发和可复现的代码参考。2. 赛题核心与问题拆解不只是数学题2.1 问题背景与物理实质题目描述了一个具有多个燃油箱的飞行器例如大型客机或运输机。飞行过程中燃油不断被发动机消耗导致全机重量减轻更重要的是燃油分布的改变会直接影响飞行器的质心位置。而质心Center of Gravity, CG的位置是飞行器静稳定性和操纵性的关键参数必须被严格控制在某个安全范围内。如果质心过于靠前或靠后轻则增加操纵难度和燃油消耗重则危及飞行安全。因此供油策略不能简单地“哪个油箱近就用哪个”而需要作为一个动态控制系统来设计在满足发动机耗油需求的同时通过在不同油箱间转移或选择消耗燃油动态调整质量分布使质心轨迹始终维持在理想区域内。这就像一边开着车一边需要不断调整车上几个水桶的水量让车的重心保持平衡同时还要尽量让车跑得远。2.2 核心难点与建模关键点这道题的挑战在于多个耦合的约束和动态过程动态性燃油消耗是一个连续时间过程质心随之连续变化。我们需要做出的是一系列随时间变化的决策何时从哪个油箱供油而不是一个静态方案。耦合约束质量守恒总耗油率等于各油箱供油率之和。质心约束全机质心位置与各油箱剩余油量及位置有关需始终在允许的上下限内。油箱容量与初始油量每个油箱有最大容量且初始油量可能不同。供油能力限制每个油箱可能有一个最大供油速率。飞行器姿态平衡有时还需考虑左右油箱对称消耗以防止滚转。优化目标通常是最大化航程在给定总油量下或最小化完成某段航程的油耗。这往往等价于在飞行全程让飞机保持一个总体更优的气动构型通常与质心位置有关。将上述物理问题转化为数学优化模型是解题的第一步也是决定后续求解难易的关键。2.3 主流建模思路选型分析我们当时调研并尝试了三种主流建模思路各有优劣思路一连续时间最优控制模型这是最“正统”也最复杂的思路。将各油箱油量作为状态变量供油率作为控制变量质心位置作为输出建立状态空间方程。目标函数是积分形式的航程指标。然后使用庞特里亚金最小值原理或动态规划求解。优点理论严谨能获得全局最优的连续时间解。缺点求解极其困难特别是对于多油箱、非线性约束的情况解析解几乎不可能数值求解如打靶法对初值敏感调试成本高在数模竞赛的有限时间内风险极大。我们的判断此路不通果断放弃。思路二离散时间动态规划将整个飞行时间离散化为若干个阶段如每分钟或每消耗一定燃油为一个阶段。在每个阶段根据当前各油箱油量状态决定下一阶段的供油策略决策使阶段代价最小并满足质心约束。最终从后向前递推或从前向后迭代寻找最优路径。优点能处理复杂约束理论上是求全局最优解。缺点“维数灾难”。假设有3个油箱每个油箱油量离散化为100个等级那么状态空间就有100^3100万个。计算量和存储需求爆炸对于稍长的航程或更精细的离散无法求解。我们的判断直接动态规划不可行但它的思想可以借鉴。思路三基于时间/事件离散的混合整数规划或非线性规划这是我们最终采用的也是实践中最可行的思路。其核心是“先离散后优化”。决策变量离散化我们不求解连续的供油率函数而是将飞行过程按固定时间间隔如Δt5分钟或固定耗油事件如每消耗总油量的1%划分为N个时段。假设在每个时段内供油策略是恒定的例如时段k内只从某1个或某几个油箱以固定比例供油。模型转化这样连续控制问题转化为一个有限维参数优化问题。决策变量变成了每个时段各油箱的供油率或供油量比例。约束转化为每个时段初/末的质心计算不等式。模型建立目标函数如总航程和约束质心上下限都可以写成这些决策变量的函数通常是非线性的因为质心计算是油量的加权平均。于是问题变成一个非线性规划问题。如果引入“是否从某油箱供油”的0-1变量则成为更复杂的混合整数非线性规划但精度可能更高。求解利用现代优化求解器如MATLAB的fmincon, Python的SciPy.optimize或更专业的IPOPT、Gurobi针对MILP线性化后进行求解。注意离散的粒度Δt或事件间隔是精度和计算量的权衡。Δt越小越接近连续解但变量越多求解越慢。需要根据赛题数据规模进行测试。3. 我们的模型构建与实现细节3.1 模型假设与符号定义在具体建模前必须明确假设这能简化问题并体现思考过程飞行条件恒定假设在优化考虑的航段内飞机飞行高度、速度、空气密度等保持不变。这样单位油耗产生的航程变化可以简化为与飞机总重相关的一个函数通常航程与重量对数成正比。燃油密度恒定忽略温度导致的燃油密度变化。油箱几何形状规则假设每个油箱的质心位置在其几何中心且不随油量变化而移动对于大型翼箱这可能不严格成立但作为简化模型可接受。供油切换瞬时完成忽略燃油管路切换带来的短暂延迟。关键符号定义n: 油箱数量T: 总飞行时间或待优化N: 离散时段数Δt: 离散时段长度T N * Δtm_i^0: 油箱i的初始油量 (kg)(x_i, y_i, z_i): 油箱i的质心坐标 (m) 在机体坐标系下m_{total}(t): 飞机总质量空机质量 剩余燃油总质量随时间变化r_i^k: 在第k个时段从油箱i的供油率 (kg/s) ——这是核心决策变量CG_{min}, CG_{max}: 允许的质心位置范围可能是x方向的前后限或包含y方向的左右限SFC: 发动机的燃油消耗率特性可能为常数或与状态相关。3.2 非线性规划模型建立我们最终建立了如下模型决策变量一个n x N的矩阵R其中元素r_i^k表示第k时段从油箱i的供油率。此外总时间T也可能作为变量如果目标是最大航程。目标函数最大化航程 航程Range可以通过积分燃油消耗与距离的关系得到。在恒定飞行条件下一个常用的简化模型是Breguet航程公式的变体航程与初始总重和最终总重的比值的对数成正比。最大化航程等价于最大化这个比值或者等价于在给定航程下最小化油耗。我们采用后者作为目标因为约束更直观。Minimize J sum_{k1}^{N} (sum_{i1}^{n} r_i^k) * Δt(总油耗) 实际上在总油量固定时最小化油耗等价于最大化利用燃油的效率这与保持有利质心相关。更精细的目标可以设为Minimize ∫ (Drag) dt阻力与质心有关但过于复杂。我们采用了一个等效目标最小化全程质心偏离理想位置如中点的加权平方和。因为保持理想质心通常意味着更小的配平阻力和更高的气动效率从而间接实现航程最优。这个目标函数是连续可导的易于优化。Minimize J sum_{k1}^{N} [ (CG_x(k) - CG_{ideal})^2 ] * w_k其中CG_x(k)是第k时段末的质心x坐标w_k是权重可设为1。约束条件质量守恒/燃油消耗每个时段的总供油率等于该时段发动机的需求耗油率d_k由飞行状态决定可作为已知输入。sum_{i1}^{n} r_i^k d_k, for k 1...N油箱油量动态m_i^{k} m_i^{k-1} - r_i^{k} * Δt, for i1...n, k1...Nm_i^{0}已知m_i^{k} 0(油量非负)质心计算 飞机总质心坐标是各部件质心的加权平均。假设空机质心(x_empty, y_empty, z_empty)和重量m_empty已知。CG_x(k) [ m_empty * x_empty sum_{i1}^{n} (m_i^{k} * x_i) ] / [ m_empty sum_{i1}^{n} m_i^{k} ]CG_y(k)同理。质心约束CG_{min, x} CG_x(k) CG_{max, x}, for k 1...N如果考虑横向平衡则对CG_y也有类似约束供油率上下限0 r_i^k r_i^{max}每个油箱最大供油能力油箱容量约束通常已由初始油量和非负约束隐含m_i^{k} Capacity_i3.3 模型求解从理论到代码我们将上述NLP问题在Python中实现并求解。核心工具是SciPy.optimize.minimize函数。步骤1问题初始化与参数设置import numpy as np from scipy.optimize import minimize, Bounds, LinearConstraint, NonlinearConstraint # 参数设置 n_tanks 3 # 油箱数量 N_stages 60 # 将飞行时间离散为60个阶段例如假设总时间3000秒每阶段50秒 dt 50.0 # 每个阶段时长 (秒) # 油箱属性 [初始油量(kg), x坐标(m), y坐标(m), 最大供油率(kg/s)] tanks np.array([ [5000, 10.0, -5.0, 0.5], # 油箱1: 左翼 [5000, 10.0, 5.0, 0.5], # 油箱2: 右翼 [8000, 0.0, 0.0, 0.8] # 油箱3: 中央油箱 ]) # 飞机空机属性 m_empty 40000 # 空机质量 (kg) cg_empty np.array([8.0, 0.0]) # 空机质心坐标 (x, y) # 质心允许范围 cg_x_min, cg_x_max 7.5, 8.5 # x方向前后限 cg_y_min, cg_y_max -0.2, 0.2 # y方向左右限对称平衡 # 发动机阶段耗油率 (kg/s), 这里假设为常数实际可根据飞行计划变化 d_rate 0.8 # 总需求耗油率 # 决策变量向量化 # 决策变量: 一个 (n_tanks * N_stages) 的一维数组 # 排列方式: [r_1^1, r_2^1, r_3^1, r_1^2, r_2^2, r_3^2, ..., r_1^N, r_2^N, r_3^N] num_vars n_tanks * N_stages x0 np.ones(num_vars) * (d_rate / n_tanks) # 初始猜测平均分配步骤2定义约束函数约束分为线性约束质量守恒和非线性约束质心范围。# 1. 线性约束每个阶段各油箱供油率之和等于总需求d_rate A_eq np.zeros((N_stages, num_vars)) for k in range(N_stages): A_eq[k, k*n_tanks:(k1)*n_tanks] 1.0 linear_constraint LinearConstraint(A_eq, lbd_rate, ubd_rate) # 等式约束 # 2. 供油率上下限约束 bounds Bounds(lbnp.zeros(num_vars), ubnp.tile(tanks[:, 3], N_stages)) # 每个油箱的供油率上限重复N次 # 3. 非线性约束质心约束 def cg_constraint(x): 计算每个阶段末的质心并返回与上下限的差值。 返回一个数组前N_stages个是CG_x的下界差值接着是上界差值... # 将决策变量重塑为 (N_stages, n_tanks) 矩阵 R x.reshape((N_stages, n_tanks)) # 计算每个阶段末的油量 m_tanks tanks[:, 0].copy() # 初始油量 cg_x_vals np.zeros(N_stages) cg_y_vals np.zeros(N_stages) for k in range(N_stages): # 更新本阶段油量消耗 m_tanks - R[k, :] * dt # 计算总质量和总力矩 total_mass m_empty np.sum(m_tanks) moment_x m_empty * cg_empty[0] np.sum(m_tanks * tanks[:, 1]) moment_y m_empty * cg_empty[1] np.sum(m_tanks * tanks[:, 2]) # 计算质心 cg_x_vals[k] moment_x / total_mass cg_y_vals[k] moment_y / total_mass # 约束形式 g(x) 0 # 对于下界 CG_x - cg_x_min 0 - g1 cg_x_vals - cg_x_min # 对于上界 cg_x_max - CG_x 0 - g2 cg_x_max - cg_x_vals g_x_lower cg_x_vals - cg_x_min g_x_upper cg_x_max - cg_x_vals g_y_lower cg_y_vals - cg_y_min g_y_upper cg_y_max - cg_y_vals # 合并所有约束 return np.concatenate([g_x_lower, g_x_upper, g_y_lower, g_y_upper]) # 非线性约束要求cg_constraint返回的所有值 0 nonlinear_constraint NonlinearConstraint(cg_constraint, lb0, ubnp.inf)步骤3定义目标函数我们采用最小化质心偏离理想位置的平方和。def objective_function(x): 目标最小化全程质心偏离理想位置(8.0, 0.0)的加权平方和。 R x.reshape((N_stages, n_tanks)) m_tanks tanks[:, 0].copy() total_cost 0.0 cg_ideal_x 8.0 # 理想质心x cg_ideal_y 0.0 # 理想质心y for k in range(N_stages): m_tanks - R[k, :] * dt total_mass m_empty np.sum(m_tanks) moment_x m_empty * cg_empty[0] np.sum(m_tanks * tanks[:, 1]) moment_y m_empty * cg_empty[1] np.sum(m_tanks * tanks[:, 2]) cg_x moment_x / total_mass cg_y moment_y / total_mass # 计算偏离代价可以给x方向更高权重 cost 1.0 * (cg_x - cg_ideal_x)**2 0.5 * (cg_y - cg_ideal_y)**2 total_cost cost return total_cost步骤4调用求解器并处理结果# 使用SLSQP算法它能处理边界约束和等式、不等式约束 result minimize(objective_function, x0, methodSLSQP, boundsbounds, constraints[linear_constraint, nonlinear_constraint], options{maxiter: 1000, ftol: 1e-6, disp: True}) if result.success: print(优化成功) R_opt result.x.reshape((N_stages, n_tanks)) # 计算并输出最终质心轨迹、剩余油量等 # ... (后文详细分析) else: print(优化失败:, result.message)4. 求解策略、技巧与结果分析4.1 求解器选择与调参经验直接使用上述模型求解可能会遇到问题变量多3*60180个非线性约束复杂容易陷入局部最优或收敛慢。我们的调优策略分步优化序列线性化/二次规划这是关键技巧。先求解一个线性化的简化模型将其结果作为NLP的初始值。第一步线性规划LP求可行解。我们将质心约束在每个时段线性近似。因为质心是油量的有理函数可以在初始点如平均供油进行一阶泰勒展开得到一个关于供油率r_i^k的线性不等式。然后求解一个以“满足所有约束”为首要目标的LP目标函数可以设为供油率变化平缓。SciPy的linprog或PuLP库可以完成。这个LP解通常是一个可行的、但不一定最优的供油计划它为NLP提供了一个极好的“热启动”点远超平均分配的初始猜测。变量缩放供油率变量r_i^k的数量级如0.几和目标函数、约束的数量级可能差异很大影响求解器数值稳定性。我们对决策变量进行了归一化处理将所有供油率除以d_rate使其大致在0~1范围内。约束容差调整质心约束的容差tol不要设得太小如1e-8对于工程问题1e-4或1e-3米级通常足够。过小的容差会增加不必要的求解难度。使用更强大的求解器对于更大规模的问题我们后来尝试了IPOPT通过cyipopt接口和Gurobi如果模型能转化为混合整数二次规划MIQP。IPOPT处理大规模NLP能力很强。Gurobi则需要我们将模型线性化或二次化并可能需要将连续供油率离散化为几个档位如0% 50% 100%最大供油率引入0-1整数变量转化为MILP/MIQP问题虽然会损失一些精度但能保证找到全局最优解对于线性/二次模型。4.2 结果可视化与策略解读优化完成后对结果的分析至关重要。我们绘制了以下几类图各油箱油量随时间变化曲线可以清晰看到优化策略是如何调度燃油的。例如可能优先消耗中央油箱的油以快速调整质心然后交叉消耗左右翼箱以保持横向平衡。import matplotlib.pyplot as plt # 计算油量变化 m_history np.zeros((N_stages1, n_tanks)) m_history[0, :] tanks[:, 0] for k in range(N_stages): m_history[k1, :] m_history[k, :] - R_opt[k, :] * dt plt.figure(figsize(10,6)) for i in range(n_tanks): plt.plot(np.arange(N_stages1)*dt, m_history[:, i], labelfTank {i1}) plt.xlabel(Time (s)) plt.ylabel(Fuel Mass (kg)) plt.legend() plt.grid(True) plt.title(Fuel Mass in Each Tank vs. Time) plt.show()飞机质心轨迹图将计算出的每个阶段末的CG_x,CG_y绘制出来并标出允许的边界框。这是检验约束是否满足的最直观方式。理想的优化结果应该是质心轨迹平滑地穿过允许区域并尽可能靠近理想线。供油率分配堆叠图展示每个时刻总耗油需求是如何由各油箱分担的。这能看出策略是“顺序供油”还是“并行供油”。典型优化策略解读 从我们多次求解的结果来看一个常见的最优模式是初期主要消耗中央油箱的燃油。因为中央油箱通常位于机身其油量变化对质心纵向x方向影响显著。快速消耗中央油箱的油可以使质心从初始位置可能偏前或偏后快速移动到理想区域中部。中期当中央油箱油量较低时开始按一定比例同时消耗左右翼箱。此时需要精细控制左右翼箱的消耗比例以维持横向y方向平衡。优化算法往往会给出一个非对称但周期性的左右切换策略以将质心“锁”在理想点附近。末期所有油箱油量均较低优化算法会更频繁地调整供油源以应对油量少时质心对供油变化更敏感的特性。4.3 灵敏度分析与模型稳健性一个好的模型不仅要给出解还要知道这个解有多“可靠”。我们做了简单的灵敏度分析离散粒度敏感性将Δt从50秒增大到100秒或减小到25秒重新求解。观察最优目标函数值的变化。如果变化很小1%说明当前离散精度足够。如果变化大则需要更细的离散。参数扰动分析将空机质心cg_empty、油箱坐标等参数在小范围内随机扰动例如±2%重新求解多次。统计质心约束的违反情况蒙特卡洛模拟。这可以评估策略对模型误差的稳健性。一个稳健的策略应该在参数轻微变化时依然能满足约束。需求波动分析实际飞行中发动机耗油率d_rate并非恒定。我们将d_rate设为随时间变化的序列如爬升阶段高巡航阶段平稳重新优化。观察优化策略是否能自适应变化。通常我们的模型框架能很好地处理这种变化只需将d_k定义为随时间变化的输入即可。5. 参赛实战心得与避坑指南5.1 论文写作与建模表述数模竞赛结果和论文各占半边天。对于这道题论文写作有几个关键点清晰的问题重述不要照抄题目。要用自己的话结合示意图阐明“质心平衡”与“供油策略”、“航程”之间的物理和数学联系。模型假设的合理性论证每一条假设如恒定飞行条件、油箱质心固定都要说明其依据“基于巡航段主要航程的考虑”、“简化分析突出主要矛盾”并讨论若放松该假设模型将如何复杂化。模型过渡的自然性从物理模型质心公式到离散化决策再到NLP模型推导要连贯。给出决策变量、目标函数、约束条件的数学公式后务必用一段文字解释其实际物理意义。算法描述的层次性不要直接贴代码。先讲清整体求解框架离散化 - NLP建模 - 分步优化策略再说明使用的工具SciPy, IPOPT等和关键算法SLSQP内点法。将核心的优化循环、约束处理逻辑用伪代码表示。结果分析的深度不要只说“我们得到了结果”。要分析结果为什么是合理的。结合前面的物理背景解释优化策略背后的逻辑如为什么先消耗中央油箱。用图表说话并对图表进行充分解读“从图5可见质心轨迹被严格控制在阴影区域内且大部分时间贴近理想中线”。模型检验与推广必须包含灵敏度分析部分展示模型的稳健性。讨论模型的优点综合考虑多约束、可扩展和局限性离散化误差、未考虑动态特性并提出可能的改进方向如结合反馈控制、更精细的燃油消耗模型。5.2 编程实现与调试坑点变量重塑错误决策变量从一维数组x重塑为(N_stages, n_tanks)矩阵时顺序非常重要C顺序还是F顺序。必须保证在目标函数和约束函数中重塑方式一致。我们在这里栽过跟头导致梯度计算错误求解器无法收敛。建议在代码中显式注释变量排列顺序并编写一个小型测试函数验证重塑和计算的正确性。约束违反的“软”处理严格来说质心约束必须每时每刻满足。但在数值求解中有时在迭代初期或由于数值误差会有轻微违反。如果求解器报告“约束不满足”可以尝试略微放宽约束边界如将cg_x_min从7.5改为7.49给求解器一点空间。在目标函数中加入惩罚项。例如将原目标改为J ρ * sum(max(0, cg_min - CG)^2 max(0, CG - cg_max)^2)其中ρ是一个很大的惩罚系数。这可以将硬约束转化为软约束有时更容易收敛然后再用得到的结果作为初始值去求解硬约束问题。求解器停滞与局部最优SLSQP有时会陷入局部最优。我们的应对方法是多初始点尝试除了平均分配初始值还可以尝试其他启发式策略的初始值例如“只从中央油箱开始供油”的策略。调整求解器参数增加maxiter减小ftol或尝试其他方法如trust-constr。简化问题先求解一个更粗的离散化问题N_stages较小将其解插值作为精细问题的初始值。代码效率在约束函数中避免使用循环计算每个阶段的油量和质心。虽然上面示例用了循环以便理解但在实际比赛中对于N很大的情况应使用向量化操作。例如可以预先计算好油量变化的累积和矩阵一次性算出所有阶段的油量和质心能大幅提升计算速度。5.3 团队协作与时间管理这道题工作量不小三天时间非常紧张。第一天上午必须完成选题、彻底理解问题、确定基本建模路线我们确定了离散化NLP的路线。同时一人开始搜集资料Breguet公式、飞机配平知识一人开始搭建最基本的Python环境安装SciPy, matplotlib等和代码框架定义参数、变量结构。第一天下午至晚上集中火力推导数学模型写出所有公式并开始编写核心的目标函数和约束函数代码。完成第一个可运行的“版本”即使结果不对。第二天全天调试代码解决收敛性问题。实现结果的可视化。开始撰写论文的“问题分析”、“模型假设”、“模型建立”部分。切忌等代码完全调通再写论文。第三天优化求解策略如加入分步优化进行灵敏度分析完善所有图表。全力撰写和润色论文包括“结果分析”、“模型检验”、“优缺点讨论”。最后留出2-3小时做全文统稿、检查格式、生成最终PDF。最重要的心得保持沟通每天至少开两次短会同步进度。将大问题分解为建模、编程、写作三个相对独立但又有交集的任务让队友各司其职又紧密配合。遇到卡壳如求解器不收敛不要纠结超过2小时及时团队讨论切换思路。这道F题的魅力就在于它没有唯一的标准答案任何一个逻辑自洽、求解稳定、分析深入的模型都能获得评委的青睐。