从物理动力学到策略优化：自行车运动员能量建模与MATLAB实现

1. 从一道赛题到一套方法论：自行车运动员能量特征建模的深度复盘

去年带学生备赛美赛，A题“自行车运动员的能量特征”让不少队伍直呼“物理劝退”。题目本身并不复杂，核心是建立一个数学模型，描述运动员在给定功率输出下，速度、时间、能量消耗之间的关系，并分析不同骑行策略（如恒定功率、先快后慢等）对总成绩和能量分配的影响。但难点在于，它要求你从一个经典的物理动力学模型出发，结合生理学限制，构建一个能用于策略优化的、有实际指导意义的模型。这恰恰是数学建模竞赛的精髓——不是解一道数学题，而是用数学工具解决一个简化但完整的现实问题。很多同学拿到题目，要么一头扎进微分方程求解，忽略了模型的生理学解释和策略分析部分；要么在策略优化部分想得太复杂，试图引入机器学习等“重型武器”，结果模型臃肿，求解困难，偏离了题目考察的核心。今天，我就以这道题为例，完整复盘从问题理解、模型建立、求解到论文撰写的全过程，并附上核心的MATLAB程序框架。无论你是正在备战美赛、国赛，还是对“如何用数学建模解决实际问题”感兴趣，这篇深度解析都能帮你理清思路，避开我们当年踩过的坑。

2. 问题拆解：我们到底要建一个什么样的模型？

题目描述了一个自行车运动员在一条平直赛道上骑行6英里（约9.66公里）。给出了运动员的质量、自行车质量、空气阻力系数、滚动阻力系数、传动效率等物理参数。核心要求是：建立模型，确定运动员以恒定功率输出骑行所需的时间，并探讨更复杂的骑行策略。

2.1 核心物理模型：从牛顿第二定律出发

这是整个问题的基石。自行车运动员的骑行过程，本质上是一个受多种力作用的动力学问题。我们首先需要建立运动方程。

运动员和自行车系统受到的主要力有：

1. 驱动力 (F_propulsive)：由运动员踩踏产生，通过传动系统作用于后轮。题目给出了功率P，根据物理公式，瞬时驱动力F_propulsive = (η * P) / v，其中η是传动效率，v是瞬时速度。这里有一个关键点：功率恒定不代表力恒定。在低速时，同样的功率能产生更大的力；高速时，力会变小。这是模型非线性的来源之一。
2. 空气阻力 (F_air)：与速度的平方成正比，方向与运动方向相反。公式为F_air = 0.5 * ρ * CdA * v^2，其中ρ是空气密度，CdA是空气阻力系数与迎风面积的乘积。这是模型中最重要的阻力项，尤其在高速阶段。
3. 滚动阻力 (F_rolling)：近似为常数，与正压力成正比，方向与运动方向相反。公式为F_rolling = Crr * m * g，其中Crr是滚动阻力系数，m是总质量，g是重力加速度。
4. 惯性力：在加速或减速时体现，即m * a（质量乘以加速度）。

根据牛顿第二定律，沿运动方向有：F_propulsive - F_air - F_rolling = m * a将各项代入，我们得到关于速度v(t)的微分方程：(η * P) / v - 0.5 * ρ * CdA * v^2 - Crr * m * g = m * dv/dt

为什么是这个形式？很多同学一开始会疑惑，为什么驱动力是P/v。这需要从功率的定义理解：功率P = 力F × 速度v。因此，在某一时刻，如果输出功率P是恒定的，那么此时刻的驱动力F自然等于P/v。这个关系是瞬时的，它将功率这个“原因”与力和速度这两个“结果”联系了起来。

2.2 模型的两个核心任务与边界条件

建立微分方程只是第一步。题目要求我们完成两个层面的任务：

任务一：恒定功率情景下的时间预测。这是基础。给定一个恒定的功率值P（例如225W，300W等），我们需要求解上述微分方程，得到速度随时间变化的关系v(t)，然后计算从起点（v=0）到完成6英里距离所需的总时间T。

• 初始条件：t=0时，v=0。运动员从静止开始启动。
• 终止条件：行驶的总距离S等于6英里时，对应的时间T即为所求。
• 求解方法：这是一个一阶常微分方程初值问题。由于方程形式为dv/dt = f(v)，不显含时间t，我们可以采用数值方法求解，如经典的龙格-库塔法（Runge-Kutta）。在MATLAB中，ode45函数非常适合处理这类问题。

任务二：非恒定功率策略的探索与优化。这是题目的升华部分。运动员的功率输出不可能全程严格恒定，更现实的策略是动态分配体能。题目暗示了两种策略：先以较高功率骑行再降低功率，或者先以较低功率骑行再提高功率。我们需要模型能够评估不同功率曲线P(t)下的总成绩（时间T）和/或能量消耗。

• 能量约束：这是关键生理限制。假设运动员的总可用能量E_total是有限的。那么对于任何策略P(t)，必须满足积分约束：∫_0^T P(t) dt ≤ E_total。这引入了优化问题的约束条件。
• 优化目标：通常是最小化总时间T。即在总能量E_total的约束下，寻找最优的功率分配曲线P(t)，使得完成固定距离S的时间最短。这是一个泛函优化问题（最优控制问题）。
• 简化路径：对于美赛而言，完全求解这个最优控制问题可能过于复杂。一个更务实、也更容易出彩的做法是进行参数化策略对比。例如，定义策略为“前一半距离以功率P1骑行，后一半距离以功率P2骑行”，那么P1和P2就是两个决策变量，在满足总能量约束P1*(T1) + P2*(T2) = E_total的条件下，寻找使总时间T1+T2最小的(P1, P2)组合。这可以通过遍历搜索或简单的优化算法实现。

3. 模型求解的实战：从微分方程到代码实现

理论清晰后，我们需要将其转化为可运行的代码。这里我分享用MATLAB求解核心部分的关键步骤和代码框架。注意，以下参数值为示例，实际做题需根据题目附件数据确定。

3.1 恒定功率模型的数值求解

首先，我们需要将微分方程改写为MATLAB ODE求解器要求的标准形式。我们的方程是：dv/dt = (1/m) * [ (η*P)/v - 0.5*ρ*CdA*v^2 - Crr*m*g ]

但是，这里有一个陷阱：当速度v为0时，等式右侧第一项(η*P)/v趋于无穷大，导致计算失败。这是模型在启动瞬间的一个奇点。在实际物理中，启动瞬间需要极大的力，但我们的模型基于“功率恒定”的假设在v=0时失效。

如何处理这个启动奇点？一个常用且合理的处理方法是设置一个很小的初始速度v0，而不是严格的0。例如，设v0 = 0.1 m/s。这样既避免了计算错误，又对总时间结果影响微乎其微（因为从0加速到0.1m/s所需的时间极短）。另一种更精细的做法是在启动阶段假设一个最大牵引力限制，但这会大大增加模型复杂度。对于美赛，采用小初始速度的方法是完全可以接受且被广泛使用的。

% 参数定义 (示例值，需根据题目更改) m = 80 + 10; % 运动员质量 + 自行车质量 (kg) P = 300; % 恒定功率 (W) eta = 0.95; % 传动效率 rho = 1.225; % 空气密度 (kg/m^3) CdA = 0.3; % 阻力面积 (m^2) Crr = 0.005; % 滚动阻力系数 g = 9.8; % 重力加速度 (m/s^2) S_target = 6 * 1609.34; % 目标距离：6英里转换为米 % 定义微分方程函数 odefun = @(t, v) (1/m) * ( (eta * P) / max(v, 0.1) - 0.5 * rho * CdA * v^2 - Crr * m * g ); % 使用 max(v, 0.1) 来避免 v=0 导致的除零错误，这是一个简单的处理技巧 % 初始条件 v0 = 0.1; % 初始速度设为一个小正值 tspan = [0, 1000]; % 时间范围，设置一个足够大的上限 % 使用事件函数来终止积分：当总行驶距离达到S_target时停止 options = odeset('Events', @distanceEvent); [t, v, te, ve, ie] = ode45(odefun, tspan, v0, options); % 计算行驶距离：速度对时间的积分就是距离 s = cumtrapz(t, v); total_time = t(end); % 完成时间 fprintf('恒定功率 %.0f W 下，完成 %.2f 英里所需时间: %.2f 秒 (约 %.2f 分钟)\n', ... P, S_target/1609.34, total_time, total_time/60); % 事件函数定义 function [value, isterminal, direction] = distanceEvent(t, v, S_target) persistent distance; if isempty(distance) distance = 0; end % 累积距离（这是一个简化处理，更精确应在主循环中计算） % 在实际编程中，更推荐在主脚本中计算距离并判断 value = distance - S_target; % 当 value >= 0 时触发事件 isterminal = 1; % 事件发生时终止积分 direction = 1; % 只关心正向过零点 end

注意：上面的distanceEvent函数是一个概念示意。在实际编码中，更清晰的做法是在调用ode45后，先积分得到距离s = cumtrapz(t, v)，然后通过插值找到s首次超过S_target时对应的精确时间。或者，可以使用更精确的事件函数，在积分过程中动态计算并判断距离。

3.2 策略对比模型的实现

我们来实现一个简单的两段式功率策略分析。假设运动员总能量E_total固定，策略是：前一段距离S1（或时间T1）内以功率P1骑行，剩余距离内以功率P2骑行。

% 假设总能量 E_total = P_avg * T_estimated % 先估算一个平均功率下的总时间，作为能量估算基准 P_avg = 250; % 平均功率估计 (W) % ... 使用上面代码计算在P_avg下的总时间 T_est % 假设 T_est 已计算出 E_total = P_avg * T_est; % 总可用能量 (J) % 定义策略：前一半距离用高功率，后一半用低功率 S1 = S_target / 2; P1 = 350; % 高功率 (W) P2 = (E_total - P1 * T1) / T2; % 根据总能量约束计算P2，但T1, T2未知，形成耦合 % 这形成了一个循环依赖问题。我们需要一个迭代或优化的框架。 % 更直接的方法是：遍历不同的P1，对于每个P1： % 1. 用P1积分方程，直到行驶距离达到S1，记录所用时间T1和此时的速度v_mid。 % 2. 计算剩余能量: E_remaining = E_total - P1 * T1。 % 3. 用P2（作为变量）积分方程，从速度v_mid开始，初始时间从T1开始，直到行驶距离达到S_target。 % 4. 在积分过程中，需要保证消耗的能量不超过E_remaining。这要求P2是随时间变化的？ % 5. 实际上，如果第二段也采用恒定功率，那么P2必须等于 E_remaining / T2。 % 这揭示了问题的复杂性。一个更可行的简化方案是：固定两段的总时间分配比例或距离比例，然后优化P1和P2。 % 例如，我们固定S1 = S2 = S_target/2。 % 步骤： % Step 1: 对于给定的P1，计算从起点到S1所需时间T1及末端速度v_mid。 % Step 2: 剩余距离 S2 = S_target - S1。 % Step 3: 问题转化为：从速度v_mid开始，在剩余能量 E_rem = E_total - P1*T1 的约束下，以恒定功率P2骑行完S2，所需时间T2是多少？ % 这需要求解一个方程：找到P2和T2，使得同时满足： % a) 动力学方程从v_mid开始，以功率P2积分T2时间后，行驶距离为S2。 % b) P2 * T2 = E_rem。 % 这是一个关于P2（或T2）的非线性方程，可以用fzero求解。 % 由于代码较长，这里给出核心的求解逻辑框架： % 1. 定义函数：给定第一段功率P1，计算总时间T_total function T_total = two_stage_time(P1, S_target, E_total, other_params) % 解算第一段 [T1, v_mid, S1] = simulate_stage(P_start, v_init, S_stage, P_const) % 这里simulate_stage是一个模拟单段恒定功率骑行的函数，返回时间、末速度、实际距离（应等于S_stage） E_used1 = P1 * T1; E_rem = E_total - E_used1; if E_rem <= 0 T_total = Inf; % 能量不足，返回无穷大时间 return; end % 2. 第二段：需要求解满足能量和距离约束的P2 % 定义方程：f(P2) = 实际行驶距离(S2) - 目标剩余距离(S_target - S1) % 其中，实际行驶距离需要通过模拟以P2从v_mid开始骑行得到，并且模拟停止的条件是消耗能量达到E_rem或距离达到目标。 % 这是一个更复杂的边界值问题。一个实用的近似方法是：假设第二段功率P2恒定，则T2 = E_rem / P2。 % 然后模拟从v_mid开始，以功率P2骑行T2时间，看行驶了多少距离S_simulated。 % 我们的目标是让S_simulated 等于 S_target - S1。 % 因此，可以定义误差函数：error(P2) = S_simulated(P2) - (S_target - S1) % 用fzero寻找使error(P2)=0的P2。 % ... fzero求解过程 ... P2_opt = fzero(@(P2) calc_distance_error(P2, v_mid, E_rem, S_target-S1, other_params), [1, 1000]); % 3. 计算第二段最优时间 T2_opt = E_rem / P2_opt; % 4. 总时间 T_total = T1 + T2_opt; end % 然后，我们可以遍历一个合理的P1范围（例如200W到400W），调用two_stage_time函数，找到使T_total最小的P1_opt。 [P1_opt, min_time] = fminbnd(@(P1) two_stage_time(P1, S_target, E_total, params), 200, 400);

这段代码框架展示了解决策略优化问题的核心思路：将连续的能量分配问题，离散化为对有限个参数（如P1, P2）的优化问题，并通过数值模拟和方程求解来评估每个参数组合的性能。在实际比赛中，你需要完善simulate_stage和calc_distance_error这两个函数，并处理好各种边界情况（如能量不足、方程无解等）。

4. 结果分析与论文写作的“加分项”

算出结果只是完成了技术部分。如何将你的工作清晰、有深度地展现在论文中，才是决定奖项高低的关键。

4.1 如何呈现你的结果？

1. 恒定功率分析：制作一个表格，展示不同恒定功率（如200W, 250W, 300W, 350W）下的预测完成时间。同时，绘制速度-时间曲线图和速度-距离曲线图。在图中可以清晰看到，初始加速阶段速度上升较快，随后因空气阻力增大而趋于一个稳定值（平衡速度）。指出这个平衡速度：当加速度为0时，由(η*P)/v = 0.5*ρ*CdA*v^2 + Crr*m*g解出v。这个解析解可以作为验证你数值解正确性的一个基准。
2. 策略对比分析：这是亮点。用图表对比“先快后慢”（High-Low）、“先慢后快”（Low-High）、“匀速功率”（Constant）三种策略。
   • 图表1：功率-距离曲线。直观展示三种策略的功率分配。
   • 图表2：速度-距离曲线。可以看到High-Low策略前期速度高，但后期因疲劳（模型中体现为能量约束）速度下降明显；Low-High策略则相反。
   • 图表3：总时间对比。用柱状图清晰显示哪种策略在相同总能量下用时最短。通常，由于空气阻力与速度平方成正比，前期用较高功率达到较高速度的收益，会被后期因能量不足而大幅降速的损失所抵消，甚至可能不如匀速策略。你的模型结果很可能显示匀速策略或接近匀速的策略是最优的，这与很多长距离耐力运动的实际策略是相符的。这个结论本身就是一个重要的分析点。
3. 敏感性分析：美赛论文非常看重这个。探讨关键参数（如空气阻力系数CdA、总质量m、总能量E_total）的微小变化对最终成绩（时间T）的影响程度。
   • 方法：改变某个参数（如CdA增加10%），重新计算最优策略下的时间，计算相对变化率。
   • 呈现：用表格或柱状图展示。结论可能是：“运动员的总成绩对空气阻力系数最为敏感，降低CdA（通过改进姿势或装备）是提高成绩最有效的途径。” 这体现了模型的实际指导意义。

4.2 论文写作必须包含的要点

1. 模型假设的清晰阐述：必须明确列出所有假设，例如：赛道绝对平坦无风、运动员总能量恒定、传动效率恒定、忽略体温变化对代谢的影响等。并简要说明这些假设的合理性及局限性。
2. 模型的优缺点分析：在结论部分，务必客观评价你的模型。
   • 优点：基于物理原理，机理清晰；能够量化分析功率分配策略；可以方便地进行参数敏感性分析，为训练和装备选择提供参考。
   • 缺点/局限性：忽略了生理上的无氧阈、乳酸堆积导致的功率衰减（我们只用总能量约束简化了）；假设功率可以瞬时切换，现实中人体有惯性；模型未考虑坡度、风向等环境因素。提出改进方向：可以引入一个关于疲劳的微分方程，将功率P表示为剩余体能的函数，从而构建更复杂的生理学-动力学耦合模型。
3. 摘要的写法：摘要必须独立成文，概括问题、方法、模型、主要结果和结论。采用“我们建立了...模型，该模型基于...动力学方程，考虑了...约束。通过...方法求解，我们发现...。敏感性分析表明...。最后，我们讨论了模型的局限性与扩展方向。”这样的结构。避免在摘要中出现公式和图表引用。

5. 常见踩坑点与高阶拓展思路

回顾我们和众多参赛队的经历，以下几个坑几乎每年都有人掉进去：

1. 忽视单位换算：美赛题目常用英制单位（英里、磅），而物理公式国际标准单位是米、千克、秒。任何涉及长度的参数（如距离、高度）、质量的参数，在代入计算前，必须统一转换为国际单位制（SI）。这是最低级也最致命的错误。在程序开头就做好所有单位的转换和注释。
2. 混淆平均功率与瞬时功率：在策略分析中，最容易犯的错误是直接用“平均功率”去分配能量。例如，认为总能量E_total = P_average * T，然后就随意分配P1和P2，忽略了动力学方程的非线性。正确的做法是，任何功率策略P(t)都必须代入微分方程进行积分，得到的速度和距离才是真实的。平均功率相等，不代表成绩相同，因为空气阻力的非线性效应。
3. 能量约束处理不当：有的队伍只做了不同恒定功率的对比，完全没有考虑总能量约束，这相当于让运动员拥有无限体能，失去了策略优化的意义。另一些队伍虽然引入了能量约束，但只是在策略对比后简单检查一下总能耗是否超标，而没有将能量约束作为优化问题的前提条件。必须将能量约束融入到策略生成或评估的过程中，如上文代码框架所示。
4. 模型求解方法单一或错误：对于微分方程，只知道ode45是好的，但必须理解其原理和适用性。此外，对于策略优化部分，如果采用离散化方法（如将全程分为N小段，每段功率恒定），那么问题就转化为一个非线性规划问题，可以使用MATLAB的fmincon等优化工具箱求解。这比手动遍历两段策略更通用，能探索更复杂的功率曲线。在论文中，可以简要提及这种更通用的方法，作为模型的扩展。

高阶拓展思路（用于冲击更高奖项）：

• 引入坡度：将赛道建模为有坡度的，重力分量会成为阻力或动力的一部分。运动方程需修改为F_propulsive - F_air - F_rolling - m*g*sin(θ) = m*a，其中θ是坡度角。这会使问题变成距离s的函数（因为坡度θ(s)随位置变化），方程变为dv/dt = f(v, s)，求解更复杂，但模型实用性大增。
• 动态生理模型：将总能量约束替换为一个“体能池”动态模型。例如，设最大可持续功率（FTP）为一个基准，超过FTP的输出会快速消耗体能储备，低于FTP的输出可以缓慢恢复体能。这需要建立另一个关于体能的微分方程，与运动方程耦合。这更贴近真实运动生理学，但参数估计和求解难度也大大增加。
• 蒙特卡洛模拟：考虑参数（如CdA, Crr, 总能量）的不确定性，将其视为符合某种分布的随机变量，进行成千上万次模拟，得到完成时间的概率分布。这可以回答“在95%的置信水平下，运动员完成时间不超过多少？”这类更丰富的问题。

这道A题是一个经典的力学建模问题，它完美地展示了如何将物理原理、数学工具和实际问题结合起来。其核心价值不在于复杂的算法，而在于清晰的建模思维、合理的简化能力以及对结果深刻的物理解释。当你拿到赛题时，不妨先问自己：这个问题的本质是什么？有哪些核心变量和关系？我可以做出哪些合理假设来简化它？我的模型结果能告诉人们什么以前不知道的事情？想清楚这些，你就已经走在正确的路上了。