简介这是一份数学建模竞赛与课程设计常用的经典案例课件面向需要掌握优化建模、参数估计与MatLab求解的本科生和建模参赛者。内容围绕“降落伞的选择”实际工程问题展开从空投2000kg物资、落地速度不超过20m/s的约束出发建立以总费用最小为目标的有约束优化模型并通过试验数据估计伞面价格与空气阻力系数等参数最后利用优化工具箱求解得到最优半径与伞数并完成落地速度验证。全套资源以1个PPT演示文稿呈现文件类型简洁、容量约1.79MB便于直接用于课堂讲解或赛前自学。课件包含问题提出、模型假设、模型建立、参数估计、模型求解与结果验证的完整逻辑链并给出了目标函数、阻力方程、非线性最小二乘拟合等关键数学推导适合希望快速理解“从实际问题到数学模型再到数值求解”全过程的读者。该资源已有1489人学习浏览是数学建模入门与竞赛训练中值得参考的优化案例材料。1. 降落伞的选择一道包装成物理题的整数规划“降落伞的选择”是数学建模竞赛里出现频率极高的经典题型。题目通常会给出若干种不同半径的降落伞、各自的价格和载重能力要求为一批物资设计空投方案使得在满足安全落地速度、载重等约束的前提下总费用最低。很多人看到题面先去做受力分析算出每顶伞的落地速度然后卡在“不同半径的伞怎么混着买”这一步上。实际上这道题真正在考察的是两件事第一能不能把一个物理过程抽象成可计算的数学模型第二能不能把一个带约束的费用最小化问题写成可求解的优化问题。前者对应微分方程和参数拟合后者对应整数规划或枚举算法。无论国赛、华为杯还是五一赛凡是“设备选型费用优化”字样的题目底层套路都与此一致。这篇文章就按最常见的解决方案展开先用运动学微分方程描述开伞后的下落过程用试验数据拟合阻力系数再建立费用模型并求解最后给出论文里最容易被扣分的敏感性与参数验证部分。2. 从牛顿第二定律到运动学微分方程确定开伞后的运动规律2.1 为什么不能用匀加速公式直接算落地速度中学物理里自由落体的速度公式是 (v \sqrt{2gh})这个公式只在忽略空气阻力时成立。降落伞问题的核心恰恰在于“空气阻力不可忽略”而且不是简单的一次线性阻力。降落伞打开后伞衣面积大空气对伞面的阻力主要来自压差阻力其大小与速度的平方成正比即[ f k v^2 ]其中 (k) 是阻力系数量纲为 (\mathrm{kg/m})。这样下落过程中物体同时受到重力 (mg) 和空气阻力 (kv^2)且阻力随速度增大而增大最终会达到二力平衡的状态即所谓“收尾速度”。如果用匀加速公式算得到的落地速度会远远偏大导致错误的结论明明伞能安全落地却被判定为超速。正确的做法是建立牛顿第二定律的微分方程。设物资和伞的总质量为 (m)从开伞时刻开始计时下落距离为 (y(t))速度为 (v(t))则运动方程为[ m\frac{dv}{dt} mg - k v^2 ]初始条件为 (v(0) 0)。注意这里的 (t) 是从开伞瞬间开始算的不是从飞机投放时刻开始算。题目里如果给出“离地高度500m时打开降落伞”之类的条件就等于给出了初速度为零的初始条件同时也给出了积分终止的条件(y(t) 500)。2.2 这个方程能不能解出解析式虽然阻力与速度平方成正比的模型在数学上存在解析解形式为[ v(t) \sqrt{\frac{mg}{k}} \tanh\left(t\sqrt{\frac{kg}{m}}\right) ]但在实际竞赛中我建议不要依赖这个解析式。原因有三个。第一如果题目给的是分段阻力模型比如开伞前后阻力系数不同解析解就要分段拼接非常繁琐。第二实际中阻力系数 (k) 不是直接给的而是需要通过试验数据拟合得到而试验数据往往是离散的高度或速度观测值直接套解析式反而绕弯路。第三你最终需要的不是速度函数本身而是“落地瞬间速度是否不超过安全阈值”这需要一个数值积分过程解析式反而不好处理落地时刻。因此更通用的做法是直接用数值方法求解微分方程常见的是四阶龙格-库塔法RK4。把速度方程和位移方程联立为一个一阶常微分方程组[ \begin{cases} \frac{dv}{dt} g - \frac{k}{m}v^2 \ \frac{dy}{dt} v \end{cases} ]其中 (y) 是下落距离。当 (y) 到达给定的空投高度时停止积分读取当时的 (v) 值。2.3 用 RK4 数值求解运动方程的最小代码以 Python 为例下面这份代码可以直接替换参数后用于任意规格的降落伞import numpy as np def parachute_landing_speed(m, k, H500, dt0.01): 用 RK4 求解开伞后的下落过程 m: 物资伞的总质量, kg k: 阻力系数, kg/m H: 空投高度, m即开伞时离地高度 dt: 时间步长, s 返回: 落地瞬间速度 m/s g 9.8 v 0.0 # 初速度开伞瞬间视为静止开始 y 0.0 # 已下落的距离 def dvdt(v): return g - (k / m) * v * v while y H: # RK4 积分速度 k1 dvdt(v) k2 dvdt(v 0.5 * dt * k1) k3 dvdt(v 0.5 * dt * k2) k4 dvdt(v dt * k3) v_new v (dt / 6.0) * (k1 2*k2 2*k3 k4) # 对位移做梯形积分避免速度突变引起误差 y_new y 0.5 * dt * (v v_new) v, y v_new, y_new # 防止陷入死循环 if y H or v 100: break # 线性插值修正跨过 H 的那一步 if y_new H: frac (H - y) / (y_new - y) v v frac * (v_new - v) return v # 示例空投 500kg 物资使用阻力系数 k30 的伞 speed parachute_landing_speed(m500, k30, H500) print(f落地速度: {speed:.2f} m/s)代码逻辑分三部分。首先是dvdt函数它对应运动方程右端的加速度表达式 (g - (k/m)v^2)注意这里已经对质量做了归一化k 和 m 的单位要匹配。其次是 RK4 积分循环每步先算四个斜率再合成新速度这是标准的四步龙格-库塔格式位移用梯形公式而不是再套一次 RK4是因为位移是速度的积分梯形公式在步长足够小的时候精度已经满足需求。最后是跨过终点的线性插值修正这一步很容易被忽略当 (y) 在某一瞬间从 499.9 跳到 500.2 时直接取该时刻的速度会带来一个步长量级的误差插值后误差被控制在可忽略范围。参数说明dt的选取直接决定精度和耗时的平衡。对于这种下落十几秒的物理过程dt0.01已经足够如果嫌慢可以放大到 0.02误差变化在 2% 以内。H是空投高度不同题目给出的开伞高度不同有的不是从投放点算起需要先判断题目里“开伞时离地高度”这个条件是否已经隐含了“开伞前自由落体”的加速过程。如果开伞前有一段自由落体初速度就不为零需要先用 (v \sqrt{2g(H-H_0)}) 算出开伞时刻的速度再传入函数。3. 用最小二乘法从试验数据拟合阻力系数 k3.1 试验数据的形态与拟合目标阻力系数 (k) 在题目里通常不会直接给出而是以一个试验数据表的面目出现。常见的表格形式是给定某规格降落伞记录不同时刻或不同高度对应的下落速度。例如时间 (t) (s)03691215速度 (v) (m/s)018.530.137.241.042.5这个速度序列有一个明显特点增速逐渐减慢最终趋向一个稳定的极限值。这个极限值就是收尾速度 (v_T \sqrt{mg/k})。拟合的目标是找到一个 (k)使得运动方程解出的速度曲线尽可能贴合这些观测点。需要说明的是有的题目给的不是速度-时间表而是速度-高度表甚至给出的是“不同高度测得的速度”。这时拟合思路仍然一样只是要把观测数据从“以时间为自变量”转换成“以高度为自变量”再比较。最常见的做法是假设阻力模型为 (kv^2)把表格数据代入运动方程构造残差平方和再用最小二乘求最优 (k)。3.2 把拟合问题写成可计算的损失函数这里有一个容易踩的坑如果直接用观测速度 (v_i) 和数值解在对应时刻的速度 (v(t_i)) 做差那么你需要先确定“对应时刻”的数值解而这依赖于一个尚未确定的 (k)造成了嵌套循环。更简单的做法是不比较速度而是比较“加速度”或“速度变化率”。由运动方程 (m dv/dt mg - kv^2)变形可得[ k \frac{mg - m\frac{dv}{dt}}{v^2} ]如果试验表足够密可以用数值微分近似 (dv/dt)代入每个观测点算出对应的 (k_i)再取平均或中位数。但数值微分对噪声很敏感数据稍微抖一下结果就会偏离很远。所以更稳健的做法还是用最小二乘直接搜索 (k) 使得模型预测的速度序列与观测序列的误差平方和最小。写成数学形式就是[ \min_k \sum_{i1}^{n} \left( v(t_i; k) - v_i \right)^2 ]其中 (v(t_i; k)) 是给定 (k) 后用 2.3 节的 RK4 方法积分出的理论速度。这本质上是一个一维优化问题可以用scipy.optimize.minimize_scalar求解。3.3 Python 实现用 scipy 拟合阻力系数import numpy as np from scipy.optimize import minimize_scalar # 试验观测数据时间与速度 t_obs np.array([0, 3, 6, 9, 12, 15]) v_obs np.array([0, 18.5, 30.1, 37.2, 41.0, 42.5]) # 物资总质量需要结合题目给的伞物资重量 m 500 # 重新定义求解函数返回速度序列而不是落地速度 def velocity_profile(k, m, t_obs, dt0.01): g 9.8 v 0.0 t 0.0 velocities [] idx 0 while idx len(t_obs): if t t_obs[idx]: velocities.append(v) idx 1 if idx len(t_obs): break # RK4 积分同前省略重复代码实际实现可直接调用前文函数 k1 g - (k / m) * v * v k2 g - (k / m) * (v 0.5 * dt * k1)**2 k3 g - (k / m) * (v 0.5 * dt * k2)**2 k4 g - (k / m) * (v dt * k3)**2 v_new v (dt / 6.0) * (k1 2*k2 2*k3 k4) v v_new t dt return np.array(velocities) def loss(k): v_pred velocity_profile(k, m, t_obs) return np.sum((v_pred - v_obs)**2) # 搜索 k 的区间物理上通常在 5~100 之间 result minimize_scalar(loss, bounds(5, 100), methodbounded) k_fit result.x print(f拟合得到的阻力系数 k {k_fit:.3f} kg/m)这段代码的思路是把 RK4 积分器和损失函数包成两个独立的函数velocity_profile负责在观测时间点采样模型速度loss负责计算残差平方和。minimize_scalar用有界搜索法在 5 到 100 之间找最优 k。选这个区间是因为(k) 太小意味着阻力可忽略不符合降落伞场景(k) 太大则收尾速度过低落地速度会远小于安全阈值也不符合典型题目设定。这里有个细节值得注意velocity_profile里用的是if t t_obs[idx]这种判断来采样而不是先在时间轴上生成密集网格再插值。这样做的好处是积分步长和采样点解耦不会因为观测时间点不在积分网格上而丢失精度。坏处是当多个观测点的时间间隔小于积分步长时会漏采样但实际题目给的时间间隔通常是秒级dt0.01完全不会触发这个问题。3.4 数据质量差时怎么办如果题目给的试验数据明显稀疏比如只有 4~5 个点最小二乘拟合出的 k 置信度很低。此时有几个处理策略。第一个策略是优先使用“收尾速度”数据。从 (v_T \sqrt{mg/k}) 反解出 (k mg/v_T^2)其中 (v_T) 可以直接从速度时间序列的末段取平均得到。这个方法只需要最后几个数据点对前段波动不敏感。第二个策略是引入多个规格的伞做联合拟合题目通常会给 2 种以上规格降落伞的试验数据而这些伞的差异主要在伞衣面积理论上 (k) 与伞衣有效面积 (A) 成正比即 (k \rho C_d A / 2)。你可以对不同规格的数据同时拟合得到统一的阻力模型再由此推算未在试验中出现的新规格伞的阻力系数。这个“跨规格建模”的技巧在评分里通常是加分项因为它体现的是物理建模能力而不是单纯的曲线拟合能力。4. 费用最小化从枚举到整数规划求解4.1 费用函数与决策变量经过第 2 章和第 3 章你已经拥有了一个核心工具给定伞的规格即阻力系数 k、物资质量 m 和空投高度 H可以精确计算落地速度。这个工具现在要喂给优化模型。题目里第 4 章的核心问题通常是有若干种规格的降落伞可选每种伞有自己的半径、价格和最大载重。需求是完成指定质量物资的空投任务要求每顶伞的落地速度不超过安全阈值通常取 20 m/s总费用最低。决策变量是每种伞的采购数量记第 (i) 种伞的采购数量为 (x_i)单位为顶。费用函数是最直观的线性加和[ \text{总费用} \sum_{i1}^{n} p_i x_i ]其中 (p_i) 是第 (i) 种伞的单价。注意这里没有把“物资质量”作为决策变量因为每顶伞装载多少物资通常也是一个决策给定总物资质量 (M)分配到各伞之间的载荷 (m_j) 可以是不同的。这会导致问题从单纯的“买几顶伞”升级为“买几顶伞每顶伞装多重”后者是一个带连续变量的混合优化问题复杂度显著提升。最基本的做法是忽略载荷分配差异假设每顶伞装载质量相同或者按最大载重装载。以经典题目为例总物资 500 kg四种伞的售价分别为 500、650、900、1100 元最大载重分别为 100、150、200、250 kg安全落地速度 20 m/s。假设每种伞都装到最大载重那么需要满足[ \sum 100x_1 150x_2 200x_3 250x_4 \ge 500 ]即在满足载重的前提下用枚举法找出使总费用最低的整数组合。4.2 枚举法的边界条件为什么不需要用复杂的组合优化算法这里有一个新手最容易犯的错误一上来就想用整数规划求解器摆出线性规划的标准形式甚至想用遗传算法。实际上这类题目里伞的规格种类通常只有 3~4 种每种伞的采购数量在可运输吨位的限制下基本不会超过 10 顶。搜索空间是 (10^4 10000) 个组合完全可以用四重循环暴力枚举一秒钟出结果。枚举法还有一个好处它天然能得到全部可行解方便你观察最优解附近的次优解这在论文的“方案对比”章节非常有用。具体枚举范围怎么定由载重约束可以推出每种伞数量的上限(x_i \le \lceil M / w_i \rceil)其中 (w_i) 是第 i 种伞的最大载重。下限是 0。总量方面伞数量太多会导致总费用必然偏高单价都是正数所以可以设置一个最大伞数上限比如总物资除以最小伞载重再加 2进一步缩小搜索空间。4.3 Python 枚举实现串行计算每顶伞的落地速度import itertools import numpy as np # 伞规格: (半径, 价格, 最大载重) parachutes [ {radius: 2.0, price: 500, capacity: 100}, {radius: 2.5, price: 650, capacity: 150}, {radius: 3.0, price: 900, capacity: 200}, {radius: 3.5, price: 1100, capacity: 250}, ] # 根据伞的半径推算阻力系数 k # 这里用简化模型k 与伞面积成正比参考半径为2m的伞 k15 # 实际比赛中应使用第3章拟合出的比例关系 def k_from_radius(r): return 15.0 * (r / 2.0)**2 M_total 500 # 总物资质量 kg V_safe 20.0 # 安全落地速度 m/s H 500 # 空投高度 m best_cost float(inf) best_plan None # 计算每种伞的采购数量上界 max_counts [int(np.ceil(M_total / p[capacity])) 1 for p in parachutes] # 枚举所有数量组合 for counts in itertools.product(*[range(mc 1) for mc in max_counts]): # 载重约束总载重必须不小于物资质量 total_capacity sum(c * p[capacity] for c, p in zip(counts, parachutes)) if total_capacity M_total: continue # 费用 cost sum(c * p[price] for c, p in zip(counts, parachutes)) # 对每一顶伞计算落地速度检查是否超标 feasible True for c, p in zip(counts, parachutes): if c 0: continue k k_from_radius(p[radius]) load p[capacity] # 每个伞装满最大载重 m_total load 20 # 物资 伞自身质量20为伞自重估值 v_land landing_speed(mm_total, kk, HH) # 复用第2章函数 if v_land V_safe: feasible False break if feasible and cost best_cost: best_cost cost best_plan counts print(f最优方案: 各伞数量 {best_plan}, 总费用 {best_cost} 元)这段代码的枚举结构可以拆成三层来看。第一层是itertools.product生成所有可能的数量组合本质上是一个多重笛卡尔积总组合数约 6 × 5 × 4 × 4 480 种非常快。第二层是载重约束检查用total_capacity M_total直接过滤掉装不下的组合这一步在费用计算之前做可以节省不必要的计算。第三层是安全速度检查对每个组合里的每种伞循环计算落地速度。代码里有一个值得讨论的取值load p[capacity]即每顶伞都按最大载重装载。这是一种偏保守的策略因为它会使落地速度偏大质量越大惯性越大同样阻力系数下更难减速到安全速度。如果题目允许载荷可调你可能需要把载荷也作为变量但那样计算量会上升。更安全的做法是在枚举前先求一下每种伞的“最大允许载荷”把落地速度阈值反代回运动方程用二分法或直接复用 RK4 求解器二分查找满足 (v_{\text{land}} \le V_{\text{safe}}) 的最大载荷质量。这样每个组合下每顶伞装多少就变成一个预计算值枚举时直接查表精度与速度兼得。4.4 从枚举到整数规划的扩展枚举法能处理 3~4 种伞、数量个位数的情况。但是如果伞的规格增加到 8 种、数量上到 50枚举空间就会膨胀到难以接受此时应该转用整数规划求解器。用scipy.optimize.milp或者pulp可以把这个模型写成标准形式决策变量 (x_i) 为整数约束为[ \sum x_i w_i \ge M ][ x_i \ge 0, x_i \in \mathbb{Z} ]目标函数是 (\min \sum p_i x_i)。但需要注意整数规划只能处理线性的载重约束无法直接表达“落地速度与载荷的非线性关系”。要把它纳入整数规划框架你需要把“安全性”转为每个伞规格的载重上限即安全约束被提前消化到参数 (w_i^{\text{max}}) 中这本质上就是 4.3 节末尾提到的预计算策略。这样处理后的整数规划模型与枚举法等价只是求解规模更大时更快。对于参赛论文建议同时展示枚举法的结果和整数规划的形式化表述前者证明方案的完备性后者展示你具备模型抽象能力。5. 敏感性分析、论文图表与答辩问答让结果从“算出来”变成“站得住”5.1 对阻力系数和伞自重做敏感性分析评审和高分论文之间最常见的分水岭是后者会对模型参数做扰动分析证明最优方案不会因为参数的小幅波动而变成不可行方案。对于降落伞问题最值得分析的参数有两个阻力系数 (k) 和伞自身质量 (m_0)。阻力系数 (k) 来自试验数据拟合天然有不确定性。做法是给 (k) 施加 ±10% 的扰动重新跑第 4 章的枚举观察最优方案是否改变、落地速度是否仍低于安全阈值。通常结论会是这样(k) 扰动最优伞数量总费用最大落地速度-10%(2, 1, 1, 0)2250 元18.7 m/s0%(2, 1, 1, 0)2250 元17.9 m/s10%(2, 1, 1, 0)2250 元16.8 m/s这样的表说明两个事实第一最优方案对 (k) 的变化不敏感具备鲁棒性第二即使阻力系数被低估 10%落地速度也仍有 1.3 m/s 的安全余量。与之相对如果结算结果显示某档 (k) 扰动下落地速度超过 20 m/s那你就要在论文里明确写出“该方案仅在阻力系数不低于某值时成立”并建议实际空投前先做一次实测。伞自重的扰动同理。题目中给出的伞自重如果有 ±5kg 的不确定性那么对最优解的影响可以通过重跑枚举得出。这里我给一个高效技巧不需要对每个扰动值重写一个脚本直接把参数作为函数的输入用列表推导或者循环批量执行即可。5.2 用三线表和曲线图提升论文质量数据呈现方式决定了论文是否容易被评审读懂。数学建模竞赛阅卷时间有限图表承担的信息传递作用比文字更大。三线表是必备技能表格内容只保留表头上下两条横线和表底一条横线表头与数据之间一条横线其余竖线全部去掉。前面几章展示的表格结构已经符合这个规范。曲线图方面最重要的是“速度-时间曲线对比图”在同一张图中画出不同伞规格的下落速度曲线并标注 20 m/s 的安全阈值线。这张图能直观说明为什么某些小伞不可用、大伞虽然贵但速度余量大。其次是“总费用-伞数关系散点图”横轴是总伞数纵轴是总费用把所有可行解点出来形成一条下包络线最优解就是包络线的最左点。这两张图配合前文的枚举代码可以让你的结果展示从“自说自话”升级为“可验证”。5.3 答辩时最常见的三个追问与应对答辩环节评委最爱问的三个问题几乎每个参赛队都会碰到。第一个问题是“你的阻力系数是从哪来的如果数据不准确怎么办”应对思路是说明自己的拟合过程第 3 章的最小二乘和跨规格联合拟合并给出敏感性分析的结果区间强调“即使 k 偏差 10%方案仍安全”。第二个问题是“为什么不用连续优化而是枚举”答案是问题规模小、枚举能获得全局最优且可解释同时也能通过受约束的整数规划形式验证枚举结果。第三个问题是“如果把空投高度改为 800m你的方案还成立吗”这类参数改动的回答最考验对模型的掌握程度。直接重新跑求解器的代码很快但评委更期待看到结论高度增加后落地速度更高还是更低由于开伞时间变长、更早达到收尾速度落地速度会略增但趋近于收尾速度极限。给出这个定性判断后再报数值结果比只报数据有力得多。这一轮的思考模式和你拿到题目后先问“物理规律是什么”是同一条思路——建模的终点永远是让外界对模型的边界和可靠性一目了然。本文还有配套的精品资源点击获取