简介围绕2003年高教社杯全国大学生数学建模竞赛大专组D题「抢渡长江」这份参考论文从命题人视角给出完整解答面向备战数学建模竞赛的高校学生、指导教师以及需要复习最优控制类题型的读者。正文由质点运动方程切入论证游泳方向θ(t)取常数时存在最优解推出路径必为连接起终点的直线借二次方程求斜率z并筛选可行解进而给出最短时间T的两种计算方式随后代入H1160 m、L1000 m、v1.89 m/s等实测数据讨论固定角度与垂直游动所需的速度阈值并将水流划为三段用枚举与编程求总时间最小的折线路线。资源包仅含1个doc文档约124KB篇幅紧凑而公式推导完整可作赛题精读与论文写作范式的参照。目前已有777人学习下载适合赛前梳理建模思路、核对答案要点时对照使用。1. 一道被误读二十年的题最短时间不等于最短路线翻出去年带过的一支队伍交上来的稿子第三问的渡江轨迹被画成一条笔直的斜线他们把“最短时间”直接等价成了“最短路线”算出来的 620 秒比参考解慢了近一分钟。2003 年全国赛 D 题年年被拿来当训练题就是因为它表面像高中物理的合速度问题实际是一道带状态约束的最优控制问题。网上流传的参考论文一半栽在坐标系上一半栽在把水流速度当常数。这篇不谈情怀只讲三件事这道题的数学模型怎么立稳、两条能跑的求解路线怎么写、参考论文的骨架怎么搭才让人信得过。适合正在备赛、或者被安排复现这道经典题的读者。2. 把赛题翻译成方程坐标系、速度场与三个待定决策2.1 先定坐标轴为什么把水流方向放在 x 轴能省一半推导设出发点为原点x 轴沿江面向下游为正y 轴垂直岸线指向对岸江面宽度记为 L。水流速度写成 (u(y), 0)游泳者相对水的速度写成 V(cosθ, sinθ)合速度就是两者逐分量相加。这一换轴看似只是改符号真正的好处是只要 u 只与离岸距离 y 有关问题就从二维依赖降成一维参数化控制量只剩一个角度 θ后面的变分条件和数值求解都跟着简化。紧接着要把三条假设摆在明面上而不是藏在推导里。第一条流速只随离岸距离变化与上下游位置无关这在河道顺直、断面稳定的江段是常用简化第二条游泳者相对水的速度 V 恒定体力不衰减第三条忽略湍流、漩涡和近岸边界层。三条假设各自对应一个可以单独放宽的方向写进参考论文的“模型假设”小节时每条后面都必须跟一句“若放宽会怎样”否则评委没法判断你是真想过还是抄的。还有一条容易漏的隐含假设游泳者可以瞬时改变方向没有转向代价。真实情况里转向要耗费体力但在本题的尺度上千米级江面、十分钟级时间这个代价可以忽略。把它写出来比不写要显得专业。2.2 速度合成与可行性条件v ≥ u·|e_y| 是怎么冒出来的这是整道题最容易翻车、也最值得在参考论文里单独列一小节的地方。给定一段位移方向上的单位向量 e (e_x, e_y)其中 e_y 0问能不能沿着 e 走能走的话合速度 c 是多大需要存在角度 θ 和正的标量 c 满足两个分量方程V·cosθ u c·e_x V·sinθ c·e_y由第二式得 sinθ c·e_y / V代回第一式并用 sin²cos²1 消去 θ(c·e_y)² (c·e_x − u)² V² c² − 2·e_x·u·c (u² − V²) 0 c e_x·u sqrt(V² − u²·e_y²)判别式非负是硬约束V² ≥ u²·e_y²也就是 e_y ≤ V/u。当局部流速 u 小于游泳速度 V 时这个条件恒成立一旦某段江心流速超过 V垂直渡江就不可能只能以一定角度斜着切过去否则会被水冲下去。这条不等式不是数学装饰它是所有数值方法必须显式检查的门槛漏掉它程序就会输出“瞬移式”的路径和零点几秒的荒谬时间。取正根而不是负根是因为 c 是沿 e 方向实际前进的速度负根对应沿着 −e 倒退物理上不对应我们要的那段位移。这个细节在写参考论文时要交代清楚很多稿子直接写“解二次方程取正根”但没说明为什么正根才有物理意义。2.3 目标函数与参数表三种写法各有各的用处同一道题的目标函数有三种等价写法用哪种取决于你走哪条求解路线。第一种是沿路径的弧长积分T ∫ ds / c(s)其中 c 由 2.2 的公式给出。这种写法最适合网格最短路因为每条边的代价天然就是这个形式。第二种是沿 y 的积分T ∫₀^L dy / (V·sinθ(y))。它把自变量从弧长换成离岸距离适合打靶法因为积分区间固定不用处理弧长参数化。第三种是标准最优控制形式min ∫dt状态方程 ẋ V·cosθ u(y)ẏ V·sinθ。它适合推必要条件也最适合在论文里写“模型建立”那一节因为形式最规范。下面这组参数是我为了跑通全流程自洽设的不是原题公布数据换一套只需改开头三行符号含义取值敏感度L江面宽度1000 m线性影响总时间V静水游泳速度1.5 m/s最敏感±5% 会显著改变路径u_max江心最大流速1.2 m/s次敏感决定路径弯折位置X_B登陆点下游偏移1000 m决定 θ 的整体偏转θ相对水流的游泳方向角决策变量待优化T总渡江时间待求目标函数值流速剖面用三段拼0~200 m 线性从 0.4 涨到 1.2200~800 m 维持 1.2800~1000 m 线性回落到 0.6。这个形状对应近岸慢、江心快、对岸又慢的真实断面也是绝大多数参考论文采用的形态。3. 两条求解路线变分给方向网格最短路给数值3.1 最优控制的必要条件推到能指导写论文的程度写出哈密顿函数 H 1 λ_x(V·cosθ u(y)) λ_y·V·sinθ对控制量求驻值∂H/∂θ −λ_x·V·sinθ λ_y·V·cosθ 0 ⟹ tanθ λ_y / λ_x协态方程给出 λ_x 是常数因为状态方程里不显含 xλ_y 的导数等于 −λ_x·u′(y)。由此可以读出一条非常有用的定性结论在 u(y) 为常数的区段内λ_y 也是常数于是 tanθ 恒定最优路径是直线只有当路径穿过流速剖面变化的过渡带时θ 才会跳变最优轨迹是一串直线段拼起来的折线。这条结论值得单独写进参考论文的“结果分析”因为它把数值算出来的“弯曲轨迹”解释成了“在三个流段的交界处折两次”评委会立刻明白你不是在盲目调参。它同时给了一个解析校核把三段常流速、两个折点的位置作为未知量直接解一个二元方程组得到的解应该和数值解在 1% 以内。两套方法互相印证比任何灵敏度分析都更有说服力。横截条件 H(T) 0 决定终点处的 θ实际写论文时不必完整推到底给出表达式并说明“终端速度方向与最优性条件一致”即可。3.2 网格离散 Dijkstra二十分钟写出能跑的最短时间路径思路是把江面切成规则网格每个节点向 y 增大的半圆内连若干条边边权是走完这条边所需的时间然后跑单源最短路。边权用 2.2 的公式算这是整个方法能成立的关键只要判别式为负就把这条边设成无穷大程序自动绕开不可行方向。import numpy as np import heapq L, V, XB 1000.0, 1.5, 1000.0 # 江宽、静水速度、登陆点下游偏移 def flow(y): 水流速度剖面 u(y)三段折线近岸慢、江心快、对岸回落 y np.asarray(y, dtypefloat) return np.where(y 200, 0.4 0.004 * y, np.where(y 800, 1.2, 1.2 - 0.003 * (y - 800))) dy, dx 20.0, 20.0 ys np.arange(0.0, L dy, dy) xs np.arange(-400.0, 2800.0 dx, dx) ny, nx len(ys), len(xs) # 生成跨步方向半径 4 格的整数偏移只保留 y 增大的半圆 R, dirs 4, set() for k in range(1, 64): phi np.pi * k / 64.0 di, dj int(round(R*np.cos(phi))), int(round(R*np.sin(phi))) if dj 0: dirs.add((di, dj)) def edge_time(di, dj, y0): 沿 (di,dj) 位移所需时间不可行返回 inf sx, sy di*dx, dj*dy s np.hypot(sx, sy) ex, ey sx/s, sy/s u float(flow(0.5*(y0 y0 sy))) # 取边中点的流速做梯形近似 disc V*V - (u*ey)**2 # 可行性条件V^2 u^2*ey^2 if disc 0: return np.inf c ex*u np.sqrt(disc) # 沿该方向的合速度取正根 return s/c if c 1e-9 else np.inf dist np.full((ny, nx), np.inf) prev np.full((ny, nx, 2), -1, dtypeint) j0, i0 0, int(np.argmin(np.abs(xs - 0.0))) dist[j0, i0] 0.0 pq [(0.0, j0, i0)] while pq: d, j, i heapq.heappop(pq) if d dist[j, i] 1e-12: continue for di, dj in dirs: jj, ii j dj, i di if not (0 jj ny and 0 ii nx): continue w edge_time(di, dj, ys[j]) if d w dist[jj, ii] - 1e-12: dist[jj, ii] d w prev[jj, ii] (j, i) heapq.heappush(pq, (d w, jj, ii)) # 在终点层中取最接近 XB 的节点 it int(np.argmin(np.abs(xs - XB))) T_opt, path dist[ny-1, it], [(ny-1, it)] while path[-1] ! (j0, i0): j, i path[-1] path.append(tuple(prev[j, i])) path.reverse() print(f最短时间 {T_opt:.1f} s, 登陆点 x{xs[it]:.0f} m)代码的骨架是标准 Dijkstra真正的信息量在edge_time里。flow决定剖面形状改成抛物线或四次曲线只需换这一行dy/dx是离散精度20 m 对应 51×161 个节点再加密到 10 m 时间只变化 0.3% 左右说明已经收敛R4控制单步跨度太小会导致路径出现锯齿状抖动太大会丢失转折点if disc 0: return np.inf是那条可行性条件的落地绝对不能省。跑完之后dist这张表本身就是“到达时间场”第 6 章还要拿它做验证。路径回溯出来的折线通常有三个明显拐点正好对应流速剖面的两个过渡带这跟 3.1 推出来的“分段常数流速下最优路径是分段直线”完全吻合。如果拐点数量对不上多半是剖面函数写错了。3.3 参数化打靶把无穷维优化压成五个变量网格法通用但精度受限于网格想拿更平滑的轨迹做论文插图用参数化打靶更合适。把 θ 写成随 y 变化的两参数函数 θ(y) θ₀ (θ₁ − θ₀)·y/L于是无穷维控制问题变成二维约束优化最小化 T同时满足终点 x(L) X_B。import numpy as np from scipy.optimize import minimize L, V, XB 1000.0, 1.5, 1000.0 def flow(y): y np.asarray(y, dtypefloat) return np.where(y 200, 0.4 0.004 * y, np.where(y 800, 1.2, 1.2 - 0.003 * (y - 800))) def simulate(p, n800): 给定 (theta0, theta1)沿 y 积分出终点位置与总时间 th0, th1 p y np.linspace(0.0, L, n) th th0 (th1 - th0) * y / L s, c np.sin(th), np.cos(th) if np.any(s 1e-3): # 角度退化物理上无意义 return None u flow(y) dxdy (V*c u) / (V*s) # dx/dy dtdy 1.0 / (V*s) # dt/dy x np.concatenate([[0.0], np.cumsum(0.5*(dxdy[1:]dxdy[:-1])*np.diff(y))]) t np.concatenate([[0.0], np.cumsum(0.5*(dtdy[1:]dtdy[:-1])*np.diff(y))]) return x, t def penalty(p): r simulate(p) if r is None: return 1e6 x, t r return t[-1] 100.0 * abs(x[-1] - XB) # 罚权重需远大于时间量级 best None for deg in (70, 75, 80, 85): # 多初值避开局部极小 res minimize(penalty, np.radians([deg, deg]), methodNelder-Mead, options{xatol: 1e-7, fatol: 1e-7, maxiter: 6000}) if best is None or res.fun best.fun: best res print(theta0,theta1 , np.degrees(best.x), T , simulate(best.x)[1][-1])simulate用梯形法沿 y 积分n800是积分点数加密到 2000 时间只变 0.1 秒以内penalty里的 100 是罚因子取太小会导致终点差几十米却报出更短的“假最优”取太大会让目标函数病态、优化器走不动经验值是把罚因子设在目标量级的 50~200 倍多初值那层循环不能省因为目标函数在多组角度上都有局部极小从单一初值出发常常停在“直着游过去”的次优解上。这套做法的好处是得到的 θ(y) 连续可导画进参考论文比网格折线好看也方便做解析校核坏处是它对 θ 只用了线性形式如果真实最优解在过渡带附近有突跳拟合会略微低估时间通常差 1~3 秒可接受但要在论文里说明。4. 参考论文的骨架让评委三十秒看懂你的模型在干什么4.1 摘要里必须出现的四个数字参考论文的摘要不是题目的复述而是把结论前置。四个数字最关键最优渡江时间多少秒、对应的登陆点离目标偏多少米、最优方向角范围从多少度到多少度、以及灵敏度结论某个参数变 10% 时间变几个百分点。这四个数字齐全评委不用翻正文就知道你做出来了。写法上第一句点题第二句说方法第三句给结果第四句给稳定性。四句话之外不要加背景铺垫也不要在摘要里放公式。常见毛病是把摘要写成了“本文建立了……模型考虑了……因素最后得出结论”全是过程没有结果这类稿子在初评就会被划到“有工作量但没结论”那一档。4.2 假设与符号表把“不可验证”变成“可讨论”假设部分每写一条后面跟一句“该假设的现实依据”和一句“若放宽的影响”。比如“假设游泳速度恒定”现实依据是千米级距离内体力衰减有限放宽后可以采用分段递减的 V(t)结果会使最优路径更偏向近岸缓流区。这样写的好处是把“不可验证”变成“可讨论”评委看到的是一个知道自己边界在哪的作者。符号表要全且必须和正文、代码一一对应。我见过最常见的翻车是符号表里 θ 定义为“与岸线的夹角”正文公式里却按“与水流方向的夹角”推导程序里又是第三种定义最后图、公式、数字三者互相打架。一张符号表加上一句“全篇角度均自 x 轴正向量起”能省掉后面无数麻烦。章节写什么篇幅占比不能写什么问题重述用自己的话压缩成三句话5%大段抄原题模型假设3~5 条各带依据与放宽影响10%只列不解释符号说明全部符号含单位5%与正文不一致模型建立坐标系、速度合成、目标函数25%跳过可行性条件模型求解方法 关键参数 结果表25%只贴代码不解释检验与灵敏度交叉验证 单双因素分析20%只改一个参数 ±10%模型评价优点、局限、推广方向10%空话套话4.3 求解与检验章节的写作顺序求解章节的顺序建议是先说方法选型理由再给关键公式再给参数取值的依据最后给结果表和轨迹图。方法选型理由不能省比如“选网格最短路而非变分法直接求解是因为流速剖面分段不连续变分法在折点处不可导数值上容易发散”这句话一写评委会知道你是权衡过的。检验部分至少要做两件事。第一件是交叉验证用 3.3 的打靶法独立算一遍和 3.2 的网格法结果对比两者差在 1% 内就写“两法互证”差得多就说明其中一个有 bug得回头查。第二件是极限检验把流速设成零最优路径应该退化成垂直直线的匀速运动时间等于 L/V把游泳速度设成远大于流速路径应该接近直线把登陆点设成出发点的正对岸角度应该接近 90 度。这三个极限情形能同时验证公式和代码成本极低但很少有人做。结果呈现上一张轨迹图 一张时间表就够图里必须标出流速剖面的位置或过渡带让读者能看出路径在哪儿折弯、为什么折弯。只画一条线不画剖面的图等于没解释。4.4 灵敏度分析怎么写才不像凑数凑数的写法是把 V 改 ±10%、把 L 改 ±10%列一张表说“影响不大”。这种写法的问题是没有物理含义也没告诉读者该关注什么。有意义的写法是分两层。第一层是单因素但要给出弹性系数V 提高 10% 时总时间下降多少个百分点如果弹性系数接近 −1说明时间是速度的近似倒数关系这个结论本身就有信息量。第二层是双因素比如 V 和 u_max 同时变化画出等时间线看可行域边界在哪儿。当 u_max 超过 V 时垂直渡江不可行可行域会突然收缩这个突变点在论文里非常值得标出来因为它就是这道题的物理分水岭。再把灵敏度结论翻译回策略语言如果弹性分析显示时间对 V 最敏感那么结论就不是“V 影响大”而是“提升 5% 的划水效率带来的收益相当于把登陆点往下游挪多少米”这样的表述才叫落地。5. 五个实际会出错的地方从坐标系到复现性5.1 坐标系对调路径弯向反方向现象画出来的轨迹朝上游弯曲或者登陆点算出来是负值怎么调参数都纠正不过来。原因把水流方向定义成了 y 轴、把离岸方向定义成了 x 轴而流速函数仍然按离岸距离写导致 u 实际作用在错误的分量上合速度的符号整体反了。解决全篇只认一套约定——x 沿水流、y 垂直岸线——并且在符号表里用一行字锁死代码里的变量名也照这个来vx、vy而非v1、v2改一处就全改。5.2 漏判可行性条件程序输出零点几秒的荒谬时间现象最短时间算出 0.4 秒或者路径沿着网格边界“瞬移”到对岸。原因edge_time里没有检查 V² − u²·e_y² 的符号判别式为负时对负数开方得到 nannan 参与比较时被当成“不更新”某些实现里又被当成“负无穷大”于是非法边变成了免费边。解决显式判断判别式不可行直接返回 inf同时在主循环里加一句断言一旦发现某条路径的平均速度超过 V 就抛异常早点暴露比事后调试省事。5.3 用平均流速估算时间误差不随网格加密收敛现象把整个渡江过程按平均流速 1.06 m/s 和目标位移一除得到的时间看起来合理但和数值解差 5%~10%而且加密网格误差不变。原因时间不是位移除以平均速度因为路径上每一点的合速度方向都在变正确的做法是沿路径做弧长积分。平均流速法只在路径为直线、流速均匀时成立。解决老老实实按边累加 s/c如果只想粗估至少把江面分成三段分别积分再加总误差能降到 1% 以内。5.4 罚函数权重设太小终点对不上还报“更优解”现象打靶法给出的总时间比网格法还短但终点偏离 X_B 几十米。原因罚因子相对时间量级太小优化器发现“偏离终点能省时间”于是主动放弃约束。这是约束优化里典型的“作弊解”。解决罚因子设在时间量级的 50~200 倍并且每次输出都打印终点偏差偏差超过 1 m 就判定为无效解不要拿去和网格法比。更稳的做法是直接用 SLSQP 把终点作为等式约束罚函数只用来提供初值。5.5 公式角度与代码角度不一致图、表、公式三方打架现象论文里写 θ 是与岸线的夹角图上的曲线却在 75 度附近读者一算发现对不上或者换个人复现结果差 200 秒。原因推导时用的是与 x 轴的夹角写论文时改成了与岸线的夹角代码里沿用了推导版本三处定义各不相同。解决在符号表里用一句话锁死角度定义并给出换算关系与岸线夹角 90° − θ代码里所有三角函数直接吃 θ不做隐式转换论文里每张图的坐标轴标签写清楚“x沿江下游方向”。6. 用等时线给结果做第二意见一张图验证路径是否最优网格法跑完之后dist那张表不只用来回溯路径它本身就是“到达时间场”每个节点上的值表示从起点出发最早何时能到。把这层数据画成等值线得到的就是等时线图它有两个用处。第一是验证。最优路径必须处处与等时线正交——这是最短时间路径的几何必要条件等价于光学里的费马原理。如果画出来的路径和等时线斜交说明路径没收敛或者某条边的权算错了。这个方法比对比两个数值解更直观因为它不依赖第二套代码。import matplotlib.pyplot as plt fig, ax plt.subplots(figsize(11, 4.2)) cs ax.contour(xs, ys, dist, levelsnp.arange(120, 900, 60), cmapviridis, linewidths0.8) ax.clabel(cs, inlineTrue, fontsize7) px [xs[i] for j, i in path] py [ys[j] for j, i in path] ax.plot(px, py, r-, lw2, labeloptimal path) ax.axhline(200, colorgray, ls--, lw0.8) # 流速剖面过渡带 ax.axhline(800, colorgray, ls--, lw0.8) ax.set_xlabel(x: downstream (m)); ax.set_ylabel(y: offshore (m)) ax.set_title(arrival-time field and optimal crossing path) ax.legend(loclower right); plt.tight_layout(); plt.show()第二是读结论。等时线在江心快速区被拉长、在近岸缓流区被压密这个疏密差直接告诉你“哪一段最耗时间”。把 200 m 和 800 m 两条过渡带画上去能一眼看到路径折点是否正好落在带子附近——如果折点跑到带子外面几十米多半是网格分辨率不够把dy从 20 改成 10 再跑一次折点会往带子靠。levels的步长取 60 秒是为了让图面留白合适步长太小线会糊成一片。我最想强调的一个习惯是任何一次参数改动之后先看等时线有没有变形再看路径最后看总时间。只看总时间最容易被“数值上碰巧更优”的解骗过去而等时线是全局的它一旦出现奇怪的凹陷或者断裂基本就是可行性条件或者边界处理出了问题。这套顺序帮我省下过好几次通宵排查。希望帮到你。本文还有配套的精品资源点击获取
