2003年数学建模D题“抢渡长江”：最优控制建模与Python求解

简介数学建模竞赛常以现实问题为背景考察模型构建与优化能力。最优控制作为运筹学与自动控制的核心分支其建模思路在路径规划、资源调度等领域有广泛应用。本文从运动学与速度合成原理出发拆解2003年全国大学生数学建模竞赛D题“抢渡长江”的完整建模过程先建立质点运动方程推导直线路径与可行性条件再针对分段流速构建目标函数通过Python二分法与网格枚举求解最短时间。聚焦游泳速度、水流速度与横向位移三者间的约束关系并指出垂直游陷阱、vu约束等常见误区最终给出可复现的代码与优化策略。旨在帮助竞赛选手掌握“从物理概念到数值优化”的解题框架提升工程实践与模型分析能力。1. 数学建模2003年全国赛D题参考论文为什么这道“抢渡长江”值得反复拆2003年高教社杯全国大学生数学建模竞赛D题“抢渡长江”表面是一个运动学问题实际上是最优控制里最经典的“最短时间渡河”变体。命题人给出的参考答案只有四页纸却把解析推导、可行性判断、分段建模和数值枚举全走了一遍。这份数学建模2003年全国赛D题参考论文直到今天依然是练手的好素材。它适合三类人准备国赛的建模队员、想补最优控制入门的工科生、带竞赛的指导老师。我拆完之后的最大感受是这道题真正的门槛不在公式而在从“能到终点”到“最短时间”之间的那条可行性边界。2. 最优控制建模先把运动方程和直线路径的推导吃透2.1 把游泳者当成质点变量与约束原题把竞渡放在平面区域内参赛者看成质点沿路线 (x(t), y(t)) 前进。游泳者相对水的速度大小为 u方向与 x 轴正向夹角为 θ(t)水流速度沿 x 轴方向记为 v可能是常数也可能随位置变化。于是绝对速度在两个方向的分量是dx/dt u cosθ(t) v dy/dt u sinθ(t)起点是 (0,0)终点是 (L,H)。给定 u、v 和终点位置要选一个 θ(t) 使到达时间 T 最小。这是一个终端固定、控制受约的最优控制问题。命题人第一句话就点出一条关键结论若 θ(t) 是连续函数则 θ(t) 等于常数时上述问题有最优解证明可参考 Leitmann 的《The Calculus of Variations and Optimal Control》。这个结论初看反直觉水流又不是处处相同为什么最优方向角是常数其实在这个模型里游速恒定、终点固定目标函数是时间积分变分法得到的极值条件会让控制量退化为常值。对参赛者来说这条结论的实际价值是后面所有分段计算里每一段都可以放心地把 θ 当成常数来处理不需要处理时变控制。2.2 直线路径的必然性从运动方程到二次方程当 u、v 为常数且 θ 恒为常数时速度分量恒定轨迹一定是直线。设到达时间为 T则有L (u cosθ v)T H u sinθ T由第二式得 T H/(u sinθ)代入第一式L u sinθ H(u cosθ v)移项得 u(L sinθ - H cosθ) H v。这就是命题人所说的“路径一定是连接起终点的直线”的代数表达式。为了方便解算令 c cotθ cosθ/sinθ。把两个运动方程相除L/H (u cosθ v)/(u sinθ) c v/(u sinθ)再用 H u sinθ T 消去 u sinθ得到L/H c vT/H即 c (L - vT)/H这个式子非常实用。已知 L、H、v、T一步算出 c进而得到 θ。比如 2002 年数据 H1160mL1000mv1.89m/s第一名成绩 T848sc(1000-1.89×848)/1160≈-0.5196对应 θ≈117.5°。再由 H u sinθ T 得 u≈1.54m/s。如果反过来已知 L,H,u,v 要求 T则需要解一个二次方程。我习惯从速度合成关系出发绝对速度的水平分量是 L/T垂直分量是 H/T于是(L/T - v)^2 (H/T)^2 u^2这个方程展开后是关于 1/T 的二次方程。给定 L,H,v,u判别式是否非负决定了是否存在可行路径。命题人的式(6)-(8)本质就是这个判别式条件只是他用了一个中间变量 z原文档扫描件里下标比较糊照着抄容易抄错。我这里直接用速度合成式形式上更不容易错。二次方程的两个根分别对应两种游法一个根对应“朝上游斜着游、但绝对速度还是指向终点”的路线另一个根对应“朝下游斜着游、绕大圈”的路线。前者时间短后者时间长。命题人说要取较大的根翻译成我的写法就是取 T 较小的那个正根。2.3 可行性条件u 到底得多快才能不被冲走把判别式大于等于零展开就得到“能游到终点”的临界游速。以 2002 年数据为例命题人算出的临界值是 u 1.43 m/s。这个数意味着什么男子1500米自由泳世界纪录平均速度约1.7m/s普通人根本达不到。但1934年那场竞渡比赛直线距离约5000m横向距离 L4864m垂直距离 H 仍为1160m流速还是1.89m/s这时临界游速只需要 u0.44m/s。所以历史上很多人能完成不是因为游得快而是因为赛道长给了水流足够的时间把人带过去。我建议你把这个临界值单独算一遍并写进论文。它有两个用途一是确定后续搜索区间的下界二是评阅时一眼就能看出你有没有理解“可行性”这个隐含约束。很多队伍会算出角度、时间却漏写临界游速评阅时很吃亏。3. 两个必算的经典算例反推游速和“垂直游”的陷阱3.1 从第一名成绩反推848秒对应什么水平2002年竞渡实际数据H1160mL1000mv1.89m/s第一名成绩848s。用 c(L-vT)/H 得 c≈-0.5196θ≈117.5°。注意这个 θ 是游泳者相对水的方向不是绝对运动方向。绝对运动方向是从起点指向终点的方向约与 x 轴正方向成49°角而相对水的方向是“朝左上方游”靠这个角度抵消水流把它往右冲的效应。再由 H u sinθ T 算出 u H/(T sin117.5°) ≈ 1.54m/s。这个速度比世界纪录平均速度低一些说明第一名选手并没有“拼命逆流”而是顺着水流斜切用更经济的速度完成渡江。这个反推过程在参赛论文里是很好的验证段落先用已知成绩检验模型再做预测。如果你把 u1.5m/s 代入速度合成式求 T会得到 T≈910s也就是15分10秒这和命题人算的结果完全一致说明模型自洽。3.2 垂直游为什么 z0 是必死路线如果游泳者始终垂直岸边游即 θ90°那么水平绝对速度只剩水流速度 v横向位移只能靠水流带。由 TL/v1000/1.89≈529s垂直方向需要 uH/T≈2.19m/s。这个速度已经超过1500m自由泳世界纪录平均速度1.7m/s职业运动员短时间能冲但要在529秒内维持这个速度横渡绝无可能。命题人给这个反例目的是戳破一个直觉误区很多人觉得“起点到终点的连线就是最优路线”但那是没有水流时的结论。有水流时游泳者必须朝上游偏一个角度。垂直游看起来方向最正实际会被冲到下游根本到不了终点。这个算例几乎每篇优秀论文都会提到用来做模型验证很合适。我在复现时会把“垂直游需要的 u”和“临界游速 u_min”放在同一张表里评阅人一眼就能看出你理解了两个约束的区别。3.3 用速度合成式快速估算最短时间一个可抄的模板除了手算我更喜欢用二分法直接求根。把 (L/T - v)^2 (H/T)^2 u^2 看成关于 T 的方程在可行区间内二分找最小正根。对 2002 年数据给定 u1.5m/s解得 T≈910s。模板如下import math def shortest_time(L, H, v, u): def f(T): return (L/T - v)**2 (H/T)**2 - u*u lo, hi 1e-6, 5000 # f(T) 随 T 增大单调递减所以二分一定能找到根 for _ in range(200): mid (lo hi) / 2 if f(mid) 0: lo mid else: hi mid return (lo hi) / 2 print(shortest_time(1000, 1160, 1.89, 1.5))代码逻辑f(T)是速度合成残差可行时存在唯一 T 使 f0。由于 f 随 T 增大而减小二分区间 [1e-6, 5000] 覆盖所有合理情况。参数说明L1000、H1160是终点坐标v1.89是流速u1.5是游泳者速度。这个函数在后面分段模型中也可以复用只要把每一段的 L、H、v、u 传进去就能算单段时间。4. 从三段等流速到连续流速用 Python 枚举最优 L24.1 问题3的模型流速分三段怎么拆命题人把河分成三段H1H3200mH2760m对应流速 v1v31.47m/sv22.11m/s游泳者速度恒为 u1.5m/s。注意 v2 u这是第二个关键约束。当一段流速大于游速时游泳者相对水的速度不足以完全抵消水流必须借助水流的横向搬运所以该段的横向位移 L2 不能太小。如果 L2 太小合成绝对速度的横向分量不足以在垂直位移 H2 内完成横向推进方程无解。三段路径为 ABCDAB 平行于 CD也就是第一段和第三段方向角相同。这是因为第一段和第三段流速相同、垂直高度相同且总横向位移对称。设 L1L3L2 待定满足 L1L2L31000m。这样优化变量只有一个 L2问题降维到一维搜索。4.2 枚举脚本扫描 L2 与 L1 网格参考论文用枚举法L2 从760m到1000m每20m一段L1、L3从120m到0m每10m一段。我复现了一个 Python 脚本完全按同样的网格扫描import math L_total 1000.0 H1 H3 200.0 H2 760.0 v1 v3 1.47 v2 2.11 u 1.5 def seg_time(H, L, v, u): # 给定垂直距离H、横向位移L、流速v、游速u # 返回 (时间, 方向角)若不可达返回 None a L*L H*H b -2.0 * L * v c v*v - u*u disc b*b - 4.0*a*c if disc 0: return None T1 (-b math.sqrt(disc)) / (2*a) T2 (-b - math.sqrt(disc)) / (2*a) Ts [t for t in (T1, T2) if t 0] if not Ts: return None T min(Ts) vx L / T # 绝对速度水平分量 vy H / T # 绝对速度垂直分量 ux vx - v # 相对水速度水平分量 theta math.degrees(math.atan2(vy, ux)) return T, theta best None for L1_int in range(0, 121, 10): L1 float(L1_int) L3 L1 L2 L_total - L1 - L3 if not (760 L2 1000 1e-9): continue r1 seg_time(H1, L1, v1, u) r2 seg_time(H2, L2, v2, u) r3 seg_time(H3, L3, v3, u) if r1 and r2 and r3: T r1[0] r2[0] r3[0] if best is None or T best[0]: best (T, L1, L2, r1, r2, r3) if best: T, L1, L2, r1, r2, r3 best print(f最优: L1L3{L1:.0f}m, L2{L2:.0f}m, T{T:.2f}s) print(fθ1θ3{r1[1]:.2f}°, θ2{r2[1]:.2f}°)代码逻辑说明seg_time解速度合成方程乘开得到关于 T 的二次方程其中 aL^2H^2b-2Lvcv^2-u^2。判别式小于0说明该段在给定 L 下不可达返回 None。两个正根中取较小时间对应“合理的、朝上游的方向”。atan2(vy, ux)求相对水的速度方向角注意ux可能是负的此时 θ 会落在第二象限对应“朝上游偏”。参数说明枚举步长 L1 从0到120m每10m一格L2 自动由 1000-2*L1 得到覆盖 L2 从1000m到760m。原论文表格里 L2800m 对应 L1100m是网格最优点。如果你想复现那张表可以把循环改成 L2 从760到1000步长20再换算 L1。运行这个脚本最优结果L1L3100mL2800m总时间 904.58sθ1θ3≈124.66°θ2≈119.19°和参考论文表格完全吻合。这个脚本把“模型公式”和“数值搜索”分开了后续换任何参数只要改 H、v、u 即可。4.3 结果分析为什么最优 L2 不是 760 而是 800从枚举表看L2760 时第二段时间约652.68s方向角约129°说明要非常朝上游顶流游。L2800 时第二段时间降到约580.37s方向角约119°虽然第二段横向距离变长了但方向更顺流绝对速度反而提高。继续增大 L2第一、三段缩短它们的时间减少但第二段因为需要更多横向位移方向角又要增大时间回升。T2 先降后升T1T3 一直在降两者之和形成一个单峰函数最优出现在 L2800m 附近。这就是典型的“分段权衡”每段单独最优加在一起不是全局最优必须放在一起算。4.4 问题4升级第一段流速连续变化路径变成抛物线问题4把第一段流速改成从岸边到深处连续变化命题人取 v2.28m/s这个值是为了保证三段平均流速仍为1.89m/s。第一段内水平绝对速度不再恒定路径不再是直线而是抛物线。推导思路第一段 H1200m流速随 y 变化运动方程写成 dx/dy (u cosθ v(y))/(u sinθ)。如果 θ 仍取常数且 v(y) 是线性函数那么 dx/dy 是 y 的线性函数积分后 x 是 y 的二次函数路径就是抛物线。命题人给出对应的二次方程并强调此时应取绝对值较小的根。与等流速情况不同抛物线解有两个分支只有靠近岸边的那一支有物理意义。我用同样的枚举思路把第一段时间替换成抛物线公式第二段仍用等流速公式第三段与第一段对称。枚举得到 L1L340mL2920mT892.56sθ1θ3≈126.87°θ2≈115.04°。这个结果比三段等流速模型的 904.58s 快了约12秒原因正是第一段更贴近岸边避开了深水区的大流速。如果想复现连续模型的精确最优可以把 L2 的枚举步长改成 1m得到 T892.5s与命题人的“编程计算可得 892.5s”一致。提示原论文表格里 L140、L2920 是 10m 网格下的最优点不是连续最优。写论文时要说明网格精度别把“网格最优”说成“全局最优”。5. 避坑指南做“抢渡长江”最容易翻车的五个点5.1 现象直接用直线距离除以游速当结果有些队伍会算出起点到终点距离约1532m然后除以1.5m/s得到1021s还觉得挺合理。实际这个做法把水流速度完全丢弃了。原因游泳者的绝对速度是相对水的速度与水流速度的矢量合成不是简单相加。直线距离除以游速只在静水中成立。解决至少要先做速度合成把问题分解为水平、垂直两个方向。更准确的做法是先算可行性再求最短时间。我在复现时第一步就是写(L/T - v)^2 (H/T)^2 u^2把这个方程放在论文开头评阅人一看就知道你懂矢量合成。5.2 现象把“垂直游”当成可行方案垂直游即 θ90°水平方向完全靠水流带。我算过需要 u≈2.19m/s超过世界纪录所以不可行。但如果有人只算时间不验证就会得出“垂直游最快”的错误结论。原因垂直游在无流时确实是“最短路径”但有流时你的绝对水平速度就是流速垂直到不了对岸的横向位置。解决把 c0 代入公式算出理论需要的游速再与实际成绩或常人体能对比。这道题里“垂直游”是一个绝佳的反例写进论文可以作为模型检验。5.3 现象枚举网格步长太大错过最优解参考论文用20m网格演示但那是为了手算方便。如果你照抄20m网格去写自己的参数可能正好错过最优点。比如把 v2 改成2.0最优 L2 可能在811m20m网格只能给你800或820。原因网格分辨率不足时目标函数在最优附近变化平缓取整误差会改变最优位置。解决先用粗网格定位大致区间再在最优邻域用1m或0.5m步长细扫。也可以直接用 scipy.optimize.minimize_scalar 做一维连续优化得到全局最优后再用网格验证论文里的表格。5.4 现象分段模型里没有检查“vu”的约束问题3中 v22.11m/su1.5m/s第二段流速大于游速。这意味着第二段不能逆流游必须让水流把自己往横方向带所以 L2 有一个最小允许值。如果代码里不检查判别式某些 L2 会算出虚数时间或者你会忽略这个约束。原因速度合成方程在判别式小于零时无实根物理上就是“到不了”。解决在seg_time里先算判别式小于0直接返回 None这样枚举时自动跳过不可达段。我第一次写就漏了导致 L2760 时第二段时间算成负数后来加了一行disc 0的判断就正常了。5.5 现象把参考论文的“最优”当成唯一标准命题人在文末特意加了一句“学生可以合理地改变数据”。如果你把参考答案里的 H、v、u 当成定死的论文就失去了建模比赛最重要的“自圆其说”环节。原因实际比赛时流速数据是估算的不同队伍对“平均流速1.89m/s”的理解不同分段方式也不同。解决做参数敏感性分析。比如把 v 从1.89改成1.7、1.8、2.0看最优 L2 怎么变。如果结论变化很大说明模型对流速估计敏感反而是一个值得写进论文的发现。我通常会在正文里放一张“v-最优L2”的表格评阅人很吃这一套。6. 验证与写作技巧把参考答案变成自己的解题框架6.1 用一维优化复核枚举结果枚举法能给出网格最优但网格不是连续最优。我在复现时用scipy.optimize.minimize_scalar对 L2 做连续优化目标函数就是分段总时间。对问题3连续最优约 L2806.33mT904.02s问题4连续最优约 L2922.9mT892.5s。这两个数字才是命题人说的“编程计算可得”。写论文时我会先给连续最优再说明枚举表只是为了展示“角度-时间”的变化趋势。这样既严谨又能体现你比网格法多做了一步。6.2 把推导过程和代码分开呈现很多参赛论文喜欢贴一堆公式然后直接给结果。更好的做法是公式只保留运动方程、速度合成方程和可行性条件其余中间代数变换放到附录。正文里用一张表给出“分段最优参数”比如 L1、L2、θ1、θ2、T。代码只作为说明“如何枚举”不要大段贴评阅人不关心你的缩进。我自己的习惯是正文给结论表和一张速度合成图附录放完整代码这样页数可控信息量不减。6.3 我的一个固定习惯从那以后我每次做这道题的复现都强制走一遍“先算临界游速、再解二次方程、再枚举、再用连续优化复核”的流程。这道题最迷惑人的地方不是模型推导而是那些看起来不起眼的根号正负号、网格步长、可行性判断。把这套流程固定下来换任何比赛题都能快速找到“可行解在哪、最优解在哪”。希望你也能把这份四页纸的参考答案拆成自己的解题框架而不是只记住几个结论——希望帮到你。本文还有配套的精品资源点击获取