关于2003D抢渡长江的研读报告

题目背景

本论文一抢渡长江为背景 需要完成在已知江宽 水流速度分布 选手游泳速度的前提下 建立模型计算出选手应采取的最佳游泳方向或路径

建模准备

建模前先推导了正弦定理和余弦定理

问题一

我对论文中的相关公式进行了推导
由论文中给出的参数得知

(10-4)

theta是游泳者的合成速度方向

论文中的

用到了正弦定理
在三角形AEF中

由于

EF = DA = v1
sin(theta) = sin(∠FAE)

因此推导出

经过变换得

(10-5)

游泳者的速度方向角

观察四边形ADEF
此处的合成速度v2（对应向量 AE）由两个分速度v0和v1​通过矢量相加的方式合成 作为矢量 速度的合成遵循平行四边形法则 即以v0和v1为邻边构成的平行四边形 其对角线表示合速度v2的大小和方向
因此∠AEF = ∠EAD(对顶角相等)
得到

在三角形AEF中 根据余弦定理

又由于

(10-6)

因此

(10-7)

根据速度公式t = s/v
在论文中其将路程分为了AP PC两段 因此公式中前半部份是AP的时间 后半部份是PC的时间 由于在PC时v0 v1方向同向 因此是v0+v1

代码复现

通过一下代码求的相对应的结果

import math

# 已知参数
a = 1000
c = 1160
v0 = 1.89  # 水的流速
v1 = 1.5  # 游泳速度

# 公式
def calculate_theta(x, c):
    """
    计算theta
    """
    # 计算theta
    theta = math.pi / 2 - math.atan(x / c)
    # 返回theta
    return theta

def calculate_AEF(theta, v0, v1):
    """
    计算AEF
    """
    AEF = math.asin(v0 * math.sin(theta) / v1)
    return AEF

def calculate_alpha(theta, AEF):
    """
    计算alpha
    """
    alpha = theta + AEF
    return alpha

def calculate_v2(v0, v1, alpha):
    """
    计算v2
    """
    cos_alpha = math.cos(alpha)
    v2 = math.sqrt(v0**2 + v1**2 + 2 * v0 * v1 * cos_alpha)
    return v2

通过一下代码可以得到当v1 = 1.5 1.89 2.11 v0 = 1.89时游泳时间与x的趋势

import math
import matplotlib.pyplot as plt

# 参数
a = 1000
c = 1160
v0 = 1.89
v1 = 1.5  # 给定游泳速度

model = SwimModel(a=a, c=c, v0=v0, v1=v1)

x_values = list(range(900, 1001, 5))
t_values = []

for x in x_values:
    result = model.forward(x)
    t_values.append(result["time"])

plt.figure(figsize=(8, 5))
plt.plot(x_values, t_values, marker='o', linestyle='-', color='blue')
plt.title("t VS x")
plt.xlabel("x")
plt.ylabel("t")
plt.grid(True)
plt.legend()
plt.show()

通过以下代码可以解得在x = 1000时游泳者的速度 方向 时间

import math
from scipy.optimize import fsolve

# 已知参数
T_target = 910
x = 1000
c = 1160
v0 = 1.89
a = 1000
# 求theta
theta = math.atan(c / x)

def time_equation(v1):

    ratio = (v0 * math.sin(theta)) / v1
    if abs(ratio) > 1:
        return 1e9

    AEF = math.asin(ratio)
    alpha = theta + AEF
    v2 = math.sqrt(v0**2 + v1**2 + 2 * v0 * v1 * math.cos(alpha))
    T_computed = math.sqrt(x**2 + c**2) / v2 + (a - x) / (v0 + v1)
    return T_computed - T_target

v1_guess = 1.0
v1_solution = fsolve(time_equation, v1_guess)[0]

print(f"所需的游泳者速度 v1 ≈ {v1_solution:.4f} m/s")

所需的游泳者速度 v1 ≈ 1.5003 m/s

import math
from scipy.optimize import fsolve

class SwimModel:
    def __init__(self, a, c, v0, v1=None):
        """
        a: 终点
        c: 距离
        v0: 水流速度
        v1: 游泳速度
        """
        self.a = a
        self.c = c
        self.v0 = v0
        self.v1 = v1

    def theta(self, x):
        return math.pi / 2 - math.atan(x / self.c)

    def AEF(self, theta, v1):
        ratio = self.v0 * math.sin(theta) / v1
        return math.asin(ratio)

    def alpha(self, theta, AEF):
        return theta + AEF

    def v2(self, v1, alpha):
        return math.sqrt(self.v0**2 + v1**2 + 2 * self.v0 * v1 * math.cos(alpha))

    def total_time(self, v2, v1, x):
        swim_time = math.sqrt(self.c**2 + x**2) / v2
        flow_time = (self.a - x) / (self.v0 + v1)
        return swim_time + flow_time

    def forward(self, x):
        """前向：已知 x 和 v1，输出 α, v2, t"""
        if self.v1 is None:
            raise ValueError("请先设置 v1")

        theta = self.theta(x)
        AEF = self.AEF(theta, self.v1)
        alpha = self.alpha(theta, AEF)
        v2 = self.v2(self.v1, alpha)
        t = self.total_time(v2, self.v1, x)
        return {
            "theta_rad": theta,
            "alpha_rad": alpha,
            "alpha_deg": math.degrees(alpha),
            "v2": v2,
            "time": t
        }

    def inverse_solve_v1(self, x, T_target):
        """已知 T_target 和 x，求 v1"""
        def equation(v1):
            theta = self.theta(x)
            AEF = self.AEF(theta, v1)
            alpha = self.alpha(theta, AEF)
            v2 = self.v2(v1, alpha)
            t = self.total_time(v2, v1, x)
            return t - T_target

        result = fsolve(equation, 1.0)[0]
        return result

model = SwimModel(a=1000, c=1160, v0=1.89, v1=1.50)
result = model.forward(x=1000)
print(f"α: {result['alpha_deg']:.2f}")
print(f"v2: {result['v2']:.4f}")
print(f"t:  {result['time']:.2f}")

α: 121.85
v2: 1.6822
t: 910.46

选手成功到达终点的概率

论文中定义概率

其意思是 在游泳速度为v1时成功的概率/在游泳速度最大时成功的概率
检测在2002年时的最低游泳成功的速度

代码复现

import numpy as np


u = 1.89  # 流速
X = 1000.0
Y = 1160.0

def v_min(u, X, Y):
    return np.sqrt(u**2 + (X/Y)**2)

v_thresh = v_min(u, X, Y)
print("2002 阈值速度 v_min =", v_thresh, "m/s")  # ~2.08 m/s

vs = np.linspace(0.5, 4.0, 36)
results = ["可成功" if v>=v_thresh else "失败" for v in vs]

for v, res in zip(vs, results):
    print(f"速度 {v:.2f} m/s: {res}")

def v_min(u, Y, X):
    return u * Y / X

# 参数设置
u_2002 = 1.3
Y_2002 = 1160 # 江宽
X_2002 = 1000 # 起点到终点的水平距离

v_threshold_1934 = v_min(u_2002, Y_2002, X_2002)
print(f"2002 年成功游过江所需最低速度为：{v_threshold_1934:} m/s")

def v_min(u, Y, X):
    return u * Y / X

# 参数设置
u_1934 = 1.2
Y_1934 = 1400 # 江宽
X_1934 = 3800 # 起点到终点的水平距离

v_threshold_1934 = v_min(u_1934, Y_1934, X_1934)
print(f"1934 年成功游过江所需最低速度为：{v_threshold_1934:} m/s")  # 约 0.44

2002 年成功游过江所需最低速度为：1.508 m/s
1934 年成功游过江所需最低速度为：0.4421052631578947 m/s

def d(v1, a, c):
    v0 = 1.89  # 水的流速
    upper = a - c * math.tan(math.acos(v1 / v0))
    d = upper / a
    return d
p_2002 = d(1.5, 1000, 1160) / d(1.7, 1000, 1160)
print(f"2002 年成功游过江的概率 p_2002 = {p_2002}")
p_1934 = d(1.5, 3800, 1400) / d(1.7, 3800, 1400)
print(f"1934 年成功游过江的概率 p_1934 = {p_1934}")

2002 年成功游过江的概率 p_2002 = 0.2538695839819621
1934 年成功游过江的概率 p_1934 = 0.8740249877014562

问题三

目标函数

求最优的成绩 因此目标函数就是最小时间(最小半江宽时间*2)

约束

第1个约束
每段的游泳时间

第2个约束
游泳者的速度 >= 水流在垂直方向的分量
v1 > v0说明游泳者能能克服水流的影响

第3个约束
游泳者相对水流方向的夹角

第4个约束
速度与水速的夹角即抵消水流的角度

第5个约束
每一段的合成速度约束

第6个约束
自由变量

第7个约束
半江宽约束

代码复现

通过代码复现得

import numpy as np
from scipy.optimize import minimize

# 参数设置
c1, c2 = 200, 380
v01, v02 = 1.47, 2.11
v1 = 1.5
a_total = 500


# 目标函数
def objective(a):
    # min f
    a1, a2 = a
    t1 = time_segment(a1, c1, v01)
    t2 = time_segment(a2, c2, v02)
    return 2 * (t1 + t2)

# 约束条件
def time_segment(ai, ci, v0i):
    # ti(ai)
    theta = np.pi / 2 - np.arctan(ai / ci)
    alpha = theta + np.arcsin((v0i * np.sin(theta)) / v1)
    v2i = np.sqrt(v0i**2 + v1**2 + 2 * v0i * v1 * np.cos(alpha))
    s = np.sqrt(ai**2 + ci**2)
    return s / v2i

def eq_constraint(a):
    # a1 + a2 = 0.5a
    return a[0] + a[1] - a_total


def ineq_constraint1(a):
    # v0i * c_i / sqrt(ai^2 + ci^2) ≤ v1
    a1 = a[0]
    return v1 - v01 * c1 / np.sqrt(a1**2 + c1**2)

def ineq_constraint2(a):
    a2 = a[1]
    return v1 - v02 * c2 / np.sqrt(a2**2 + c2**2)

def compute_alpha(ai, ci, v0i):
    # alpha(i)
    theta = np.pi / 2 - np.arctan(ai / ci)
    alpha = theta + np.arcsin((v0i * np.sin(theta)) / v1)
    return np.degrees(alpha)

a_init = [100, 400]
bounds = [(10, 490), (10, 490)]

constraints = [
    {'type': 'eq', 'fun': eq_constraint},
    {'type': 'ineq', 'fun': ineq_constraint1},
    {'type': 'ineq', 'fun': ineq_constraint2},
]

res = minimize(objective, a_init, bounds=bounds, constraints=constraints, method='SLSQP')

if res.success:
    a1, a2 = res.x
    print(f"最优解：a1 = {a1:.2f}, a2 = {a2:.2f}, a1 + a2 = {a1 + a2:.2f}")
    print(f"总时间 = {res.fun:.2f} 秒")
    print(f"α1 = {compute_alpha(a1, c1, v01):.2f}°, α2 = {compute_alpha(a2, c2, v02):.2f}°")
else:
    print("优化失败：", res.message)

最优解：a1 = 96.84, a2 = 403.16, a1 + a2 = 500.00
总时间 = 904.02 秒
α1 = 126.06°, α2 = 118.06°

问题四

问题四的公式和问题三的类似 主要的不同点是问题三的速度v0i是两个离散的数值 问题四的v0i是连续数值离散化 可以有n个
论文中给出的v0i的公式

其中k为水速斜率 其公式为

意思是水速在∆y内变化了∆v
论文中由于水速是分段给出的连续分布 因此水速的变化是线性的 因此其水速斜率是一个常数
水速在200米内变化了2.28 （岸边速度为0 ～水速最大为2.28）
因此

代码复现

def time_segment(a_i, c1i, v0i, debug=False):
    theta_i = np.pi / 2 - np.arctan(a_i / c1i)
    sin_theta = v0i * np.sin(theta_i) / v1
    
    if np.abs(sin_theta) > 1:
        if debug:
            print(f"[警告] sin_theta={sin_theta:.4f} 超出 [-1, 1]")
        sin_theta = np.clip(sin_theta, -1, 1)

    alpha_i = theta_i + np.arcsin(sin_theta)
    v2i = np.sqrt(v0i**2 + v1**2 + 2 * v0i * v1 * np.cos(alpha_i))
    s = np.sqrt(a_i**2 + c1i**2)
    time = s / v2i

    if debug:
        print(f"a = {a_i:.2f}, c1 = {c1i}, v0 = {v0i:.3f}")
        print(f"θ = {np.degrees(theta_i):.2f}°, α = {np.degrees(alpha_i):.2f}°, v2 = {v2i:.3f}, s = {s:.2f}, time = {time:.2f}\n")

    return time

a_list = [-30, 80, 420]
c1_list = [100, 100, 380]
v0i_list = [0.57, 1.71, 2.28]
v1 = 1.5
total_time = 0
for i in range(3):
    t = time_segment(a_list[i], c1_list[i], v0i_list[i])
    print(f"t{i+1}: {t:.2f} s")
    total_time += t

print(f"总时间：{2 * total_time:} 秒")

t1: 84.65 s
t2: 73.11 s
t3: 334.95 s
总时间：985.430963388645 秒

可以在论文中看到其路径是先往上游游在往下游游
这有点不符合直觉
我们假设若按照直线游

通过代码可以得到

alpha = np.arctan(1160 / 1000)

d = 1160/ np.sin(alpha)

v = 1.79/np.cos(alpha)

t = d/v

2 * t

1117.31843575419

这不是最短的时间 因此该路径看似不符合直觉但是确实是最短的时间

灵敏度分析

调整上述代码中的参数
将v1速度改为1.515和1.48

t1: 83.56 s
t2: 71.79 s
t3: 334.99 s
总时间：980.6787856865337 秒

t1: 85.77 s
t2: 74.54 s
t3: 334.89 s
总时间：990.4005321378754 秒

T对v1敏感

将k的值改为0.011和0.0118

t1: 84.65 s
t2: 73.11 s
t3: 334.95 s
总时间：985.430963388645 秒

t对k不敏感

模型拓展

由于问题四中的模型的速度建模思路为：将连续分布所在的江宽等分为若干区间 使得水速在每个区间近似为常数
因此在问题四中k是常数 水速是线性的
当水速不是线性的时候
其定义了v关于y的公式

在[0,200]区间中水流满足

v(0) = 0
v(200) = 2.28

为了满足此条件 其选择了一个开口向上的抛物线形式
解得 k = 2.28/200^2

由于f1(y)和f2(y)关于江心中心对称
在[960,1160]区间同理

选择抛物线函数的原因（论文中未提及 可以添加）

物理现象水流速度分布 在自然河流中 靠近岸边的地方水流较慢 越靠近江心水流越快 江心处水流速度达到最大值
这种分布不是连续且平滑变化的 抛物线是一种最简单的非线性函数 能够很好地模拟这种先慢后快再慢的对称分布趋势