2024美赛B题参赛分享

题目翻译

一、引言

问题重述

• 我们需要构建一个用于预测潜水艇在水域中位置如何随时间演变的模型，考虑潜水艇可能位于海底或水体中间层，并且其运动会受到海流、海水的密度以及海底地形等因素的作用.
• 问题1，我们需要思考为了降低预测不确定性，潜水艇需要定期地向主船报告哪些信息，并且配备哪些设备.
• 问题2，我们应该思考一下公司应该为主船配备哪些额外的搜寻设备以应对突发情况，同时考虑到这些设备的成本、维护、准备及使用方面的需求.
• 问题3，我们需要创建一个模型，通过分析潜水器位置信息来指导初始部署点和搜寻策略的选择，从而缩短寻找失联潜水艇所需的时间.接着我们将计算出随着时间推移和搜寻结果积累，成功定位潜水艇的可能性如何变化.
• 问题4：我们要阐明如何将该模型应用于其他旅游地点，例如加勒比海.同时讨论需要对模型进行哪些调整，以便它能够处理同一地区内多潜水艇的运动.

二、分析问题

模型采用

• 对于问题1,我们收集了爱奥尼亚海的三维模型、洋流信息等相关海域参数,根据受力情况建立了潜艇动力学方程,结合卡尔曼滤波建立了模型一：动态位置预测模型,并后建立了相关模拟程序进行仿真,同时分析了与这些预测相关的不确定性因素,最后给出建议可定期向主船发送潜水艇位置、当时环境参数等,以及可以携带减少这些不确定因素的相关设备.
• 对于问题2，我们收集了市面上常用的水下搜索设备，如Klein3000侧扫声纳、MEPUS-AUV3000L（AUV）、Falcon-DR（ROV）等。本文介绍了这些搜索设备的相关，以及搜索效率，建立了模型二：成本效益比率模型，通过比较不同设备的成本效益，得出应该携带主船部署的设备，此外本文讨论了搜救船可能需要配备的额外设备.
• 对于问题3，我们首先分析常用的搜寻模式，比较它们的效率，进而推荐出最佳搜寻模式，之后将搜索区域被划分为 个网格，建立模型三：基于贝叶斯网络理论的模型，根据搜索设备搜索到新的观测信息，使用贝叶斯定理更新概率分布，后验概率最高的网格 即为最佳部署点，最后给出了找到潜水器的概率与时间和累计的搜索结果相关的函数.
• 对于问题4，我们收集了加勒比海的三维模型、洋流信息等相关海域参数,讨论了模型扩展到加勒比海地区需要的改动，并结合模型一建立了相关模拟程序进行仿真，对于存在多个潜水艇的情况，本文使用Deep-Q-Network算法，进行对多个潜水艇目标的追踪与预测，最后通过仿真模拟验证了模型的有效性.

三、预测定位

• 动力学方程
  考虑潜水艇的动力学，使用牛顿第二定律表示潜艇的动力学方程：

  $m\frac{d^2\vec{r}}{d{t}^2}=\overrightarrow{F_b}+\overrightarrow{F_g}+\overrightarrow{F_c}-\overrightarrow{F_d}$

  $F_d=\frac{1}{2}\rho v^2{C}_{d}A$

其中，$\vec{r}$ 是潜水器的位置向量，$t$ 是潜水器运动的时间； $m$是潜水艇的质量； $F_g=mg$是潜水艇的重力;
$F_c$是水流对潜水器的作用力； 浮力 $F_b$可以通过潜水器排水量和水的密度来计算；
阻力$F_{\mathrm{a}}$可以根据潜水器的形状、表面粗糙度和运动速度来估算，忽略其他方向上的受力，取竖直方向 其中,
$\rho$是水的密度，$v$是潜水器相对于水的速度，$C_d$是阻力系数，$A$是潜水器迎水面积

可建立以下微分方程组:

$\{\begin{array}{l}m\frac{d^2\vec{r}}{d{t}^2}=\overrightarrow{F_b}+\overrightarrow{F_g}+\overrightarrow{F_c}-\overrightarrow{F_d}\\ F_g=m\cdot g\\ F_d=\frac{1}{2}\rho v^2{C}_{d}A\end{array}$

• 非线性Kalman滤波

虽然，我们用动力学方程预测了潜水器的位置，但是，在真实的物理世界里，潜水器在海洋里会受到一系列外在干扰因素的影响，这个问题一下子变得复杂起来.假设潜水器上有GPS，可以告诉我们在t时刻的位置，是不是解决问题了呢？但是也不能完全信任GPS，因为它也存在"精度误差".第一个误差，是根据时间,潜水器质量，水流对潜水器的作用力，浮力，阻力和位置产生的，被称为“过程误差”；第二个误差，是对传感器进行观测时产生的，被称为“观测误差”.所以我们希望用Kalman滤波模型进行计算，就是在既有“过程误差”，又有“观测误差”的真实世界里，寻找一个最优的估计和真实值更加接近.

我们可以将时间间隔 Δt 分成微小的步长：

$m\frac{v_{i+1}-v_i}{\Delta t}=F_{b_i}+F_{g_i}+F_{c_i}-F_{d_i}$

其中：
$v_i$是时间步长$i$ 时刻的速度;
$F_{bi}$ 是时间步长$i$时刻的浮力；
$F_{gi}$是时间步长$i$时刻的重力；
$F_{ci}$是时间步长$i$时刻水流对潜水器的作用力；
$F_{di}$是时间步长$i$时刻的阻力；
所以我们可得到非线性Kalman滤波状态方程:
$x_{i+1}=A{x}_i+B{u}_i$
$x_i=\left[\begin{array}{c}{r}_i\\ v_i\end{array}\right],u_i=\left[\begin{array}{c}{F}_{b_i}\\ F_{g_i}\\ F_{c_i}\\ F_{d_i}\end{array}\right],A=\left[\begin{array}{cc}1& \Delta t\\ 0& 1\end{array}\right],B=\frac{\Delta t}{m}\left[\begin{array}{cc}0& 0\\ 1& 0\end{array}\right]$

其中：

• $x_i$是时间步长$i$时刻的状态向量，包含潜水器的位置和速度
• $A$: 状态转移矩阵，用于描述状态向量$x$
• $B$:输入控制矩阵，用于描述外部控制输入对系统状态的影响
• $u_i$是时间步长时刻的输入向量，包含潜水器的浮力、重力、水流作用力和阻力；
  最终整理得到可以用于描述潜水艇的位置变化在不同时间点的位置方程:

$\begin{array}{rl}& m\frac{\vec{r}(t+\Delta t)-2\vec{r}(t)+\vec{r}(t-\Delta t)}{\Delta t^2}=(\rho v^2{C}_{d}A-F_c){\vec{e}}_z\Delta t\end{array}$

• 模型求解

import numpy as np
import matplotlib.pyplot as plt

t = np.linspace(1, 100, 100)  # 在1~100s内采样100次
a = 0.6  
v0 = 0 
s0 = 0  
imu_var = 120 ** 2  # IMU测量误差方差，Q~N(0, imu_var)
gps_var = 50 ** 2  # GPS测量误差方差，R~N(0, gps_var)
num_samples = t.shape[0]

# 根据理想模型推导出来的真实位置值
real_positions = [0] * num_samples
real_positions[0] = s0
# 模拟观测值，通过理论值加上观测噪声模拟获得
measure_positions = [0] * num_samples
measure_positions[0] = real_positions[0] + np.random.normal(0, gps_var ** 0.5)
predict_positions = [0] * num_samples
predict_positions[0] = measure_positions[0]
# 最优估计值，也就是卡尔曼滤波输出的真实值的近似逼近
optim_positions = [0] * num_samples
optim_positions[0] = measure_positions[0]
pos_k_1 = optim_positions[0]

predict_var = 0
for i in range(1, t.shape[0]):
    real_v = v0 + a * i
    real_pos = s0 + (v0 + real_v) * i / 2
    real_positions[i] = real_pos
    # 模拟输入数据，实际应用中从传感器测量获得
    v = real_v + np.random.normal(0, imu_var ** 0.5)
    measure_positions[i] = real_pos + np.random.normal(0, gps_var ** 0.5)
    predict_positions[i] = predict_positions[i - 1] + (v + v + a) * (i - (i - 1)) / 2
    # 利用上个时刻的位置（上一时刻的最优估计值）和速度预测当前位置
    pos_k_pred = pos_k_1 + v + a / 2
    # 更新预测数据的方差
    predict_var += gps_var
    # 求得最优估计值
    pos_k = pos_k_pred * imu_var / (predict_var + imu_var) + measure_positions[i] * predict_var / (
                predict_var + imu_var)
    # 更新
    predict_var = (predict_var * imu_var) / (predict_var + imu_var)
    pos_k_1 = pos_k
    optim_positions[i] = pos_k

plt.plot(t, real_positions, label='hypothetical model positions')
plt.plot(t, optim_positions, label='kalman filtered positions')
plt.plot(t, predict_positions, label='predicted positions')
plt.legend()
plt.show()

● 在建立模型的过程中，我们必须要考虑很多不确定因素，比如地质环境，自身姿态等等，我们尽可能地列出这些不确定因素.
●探险过程中，由于海底地形具有特殊的构成和地貌特征，使得某些区域水流速度较快，对潜水器的运动产生影响.此外，在潜水器运动的过程中，若遇到了强劲的洋流，可能会对潜水器的运动方向和速度造成影响.
●首先从GEBCO网站获取海底地形数据，使用Global Mapper pro v23.1分析海底地形局部特征以及部分参数，再使用python对数据进行可视化处理.

import numpy as np

def read_asc_file(file_path):
    """读取ASC文件，并返回地形数据"""
    with open(file_path, 'r') as file:
        lines = file.readlines()
        
        # 找到NODATA_value所在的行
        nodata_index = [i for i, line in enumerate(lines) if 'NODATA_value' in line][0]
                
        # 从NODATA_value所在的下一行开始读取数据
        data = [list(map(float, line.split())) for line in lines[nodata_index+1:]]
        data = np.array(data)        
        # 获取行数和列数
        nrows, ncols = data.shape        
        # 从ASC文件中获取X、Y、Z坐标
        x = np.arange(0, ncols, 1)
        y = np.arange(0, nrows, 1)
        x, y = np.meshgrid(x, y)
        z = data[::-1]  # 将数据反转，使得坐标原点位于左下角
        return x, y, z
def plot_3d_terrain(x, y, z):
    """绘制3D地形图"""
    fig = plt.figure(figsize=(10, 6))
    ax = fig.add_subplot(111, projection='3d')
    ax.plot_surface(x, y, z, cmap='terrain')
    ax.set_xlabel('X')
    ax.set_ylabel('Y')
    ax.set_zlabel('')
    ax.set_title('3D Terrain Map')
    plt.show()

# ASC文件路径
asc_file_path = 'gebco_2023.asc'
# 读取ASC文件
x, y, z = read_asc_file(asc_file_path)
# 绘制3D地形图
plot_3d_terrain(x, y, z)

●分析：爱奥尼亚海是地中海最深的子盆地，由几个深度超过4000 m的海沟和高原组成.由3D Terrian Map可以看出爱奥尼亚海海底跨度在1000-2000以内，探险人员需要注意自己的人身安全，以防失去与主船的通信.

四、设备选择

a)搜索设备
我们准备了以下比较常用的搜寻设备
●Klein3000侧扫声纳：侧扫声纳用于扫描水下地形和目标，通常安装在船只上，可用于海床勘测和水下搜寻.
●MEPUS-AUV3000L（AUV）：AUV通常配备有高清摄像头、声纳系统、机械臂、水下推进器等设备，能够在水下自主航行、收集数据、执行任务，并将所获取的信息传输至地面控制中心.
●Falcon-DR（ROV）：ROV可以携带各种传感器和工具，如摄像头、操纵臂、声纳和取样设备，以便执行不同的任务，通常由操作员通过遥控设备来控制，可以在水面或岸上的控制站远程操纵.

b)如果要使用这些设备，在准备和使用这些搜索设备时，需要考虑成本：
●设备成本：这些设备的购买和维护费用需要纳入预算考虑
●人员培训成本：特定的设备可能需要操作员接受培训才可以熟练使用
●储存和保养成本：为了储存和保养这些设备，可能需要专门的设备储藏室以及维护计划

设备名称	设备成本	人员培训成本	储存和保养成本
Klein4000侧扫声纳	＄13903.4	＄500	＄500
MEPUS-AUV3000L（AUV）	＄20855.1	$ 1500	＄1000
Falcon-DR（ROV）	＄10427.6	＄1200	＄1000

c)搜索效率
考虑成本的同时，使用效果也是不可或缺的,以下为搜寻设备的搜索效率

设备名称	搜索深度	搜索范围	搜索速度
Klein4000侧扫声纳	2000m	600m	170m/s
MEPUS-AUV3000L（AUV）	3000m	1852m	1m/s
Falcon-DR（ROV）	1100m	1100m	2m/s

3、模型的建立

• 假设主船和救援船总共允许携带设备个数小于等于5，根据搜集到的设备成本和设备搜救参数，我们在考虑设备可用性的同时需要尽可能节约相关的开支成本.

• 我们先建立成本效益比率模型：$CER(e)=E(e)/C(e)$

  其中，$e$表示不同的搜索设备

  $E(e)$表示设备的效益值，可以进一步细化为：

  $E(e)=a{E}_{profundity}(e)+b{E}_{radius}(e)+c{E}_{speed}(e)$

• $E_{profundity}$其中是搜索深度，$E_{radius}$是搜索范围，$E_{speed}$是搜索速度.

  a、b、c为这些因素的权重

  C(e)表示设备的总成本，可以进一步细化为：
  $C(e)=C_{equipment}(e)+C_{training}(e)+C_{storage}(e)$

4、模型的求解

• 根据数据我们可以列出:
  $s.t.\{\begin{array}{l}f(C_{\text{equipment }})=13903.4X_1+20855.1X_2+10427.6X_3\\ f(C_{\text{training }})=500X_1+1500X_2+1200X_3\\ f(C_{\text{storage }})=500X_1+1000X_2+1000X_3\end{array}$
  $s.t.\{\begin{array}{ll}f(E_{\mathrm{profundity}})=2000X_1+3000X_2+1100X_3& \\ f(E_{\mathrm{radius}})=600X_1+1852X_2+1100X_3& \\ f(E_{\mathrm{speed}})=170X_1+X_2+2X_3& \end{array}$
  $\min f_{C(e)}=C_{\text{equipment }}(e)+C_{\text{training }}(e)+C_{\text{storage }}(e)$
  $=14903.4X_1+23355.1X_2+12627.6X_3$
  $\min f_{E(e)}=a{E}_{\text{profundity }}(e)+b{E}_{\text{radius }}(e)+c{E}_{\text{speed }}(e)$
  $=(2000a+600b+170c)X_1+(3000a+1852b+c)X_2+(1100a+1100b+c)X_3$
  $s.t.\{\begin{array}{ll}C(e)=14903.4X_1+23355.1X_2+12627.6X_3& \\ X_1+X_2+X_3\le 5& \\ X_1,X_2,X_3>0& \end{array}$
  $s.t.\{\begin{array}{ll}E(e)=(2000a+600b+170c)X_1+(3000a+1852b+c)X_2+(1100a+1100b+c)X_3& \\ X_1+X_2+X_3\le 5& \\ X_1,X_2,X_3>0& \end{array}$
  $CER(e)=E(e)/C(e)=\frac{(2000a+600b+170c){\mathrm{X}}_1+(3000a+1852b+\mathrm{c}){\mathrm{X}}_2+(1100a+1100b+\mathrm{c}){\mathrm{X}}_3}{14903.4{\mathrm{X}}_1+23355.1{\mathrm{X}}_2+12627.6{\mathrm{X}}_3}$

• 对于线性整数规划问题，我们可以采用蒙特卡洛法进行估算，我们假设分别为0.4，0.3，0.3 ，然后我们继续进行分析：

• $s.t.\{\begin{array}{ll}E(e)=1031X_1+1755.9X_2+770.3X_3& \\ X_1+X_2+X_3\le 5& \\ X_1,X_2,X_3>0& \end{array}$

• $CER(e)=E(e)/C(e)=\frac{1031{\mathrm{X}}_1+1755.9{\mathrm{X}}_2+770.3{\mathrm{X}}_3}{14903.4{\mathrm{X}}_1+23355.1{\mathrm{X}}_2+12627.6{\mathrm{X}}_3}$

import numpy as np
import matplotlib.pyplot as plt

# 目标函数
# def objective_function(x):
#     return x[0]*14903.4 + x[1]*23355.1 + x[2]*12627.6 

#目标函数2
def objective_function(x):
    return x[0]*1031 + x[1]*1755.9 +x[2]*770.3

# 约束条件
def constraint1(x):
    return x[0] + x[1] + x[2] - 5


def constraint2(x):
    return x[0] 

def constraint3(x):
    return  x[1] 

def constraint4(x):
    return x[2] 

# 定义问题的上下界
bounds = [(0, 5)] * 5

# 定义蒙特卡洛模拟的次数
num_samples = 1000

# 初始化目标函数值列表
values = []

# 执行蒙特卡洛模拟
for _ in range(num_samples):
    # 生成随机解
    x = np.random.uniform(0, 5, size=5)
    
    # 检查是否满足约束条件
    if (constraint1(x) <= 0 and constraint2(x) > 0 and constraint3(x) > 0 and constraint4(x) > 0):
        # 计算目标函数值
        value = objective_function(x)
        
        # 添加到目标函数值列表中
        values.append(value)

# 绘制散点图
# 这段代码将每次随机得到的目标函数值绘制成散点图，横坐标为迭代次数，纵坐标为目标函数值。通过观察散点图，可以了解到目标函数值的分布情况。
plt.scatter(range(len(values)), values, s=5, color='#845EC2')
plt.xlabel('Iteration')
plt.ylabel('Objective Function Value')
plt.title('Monte Carlo Simulation Results')
plt.grid(True)
plt.show()

• 对设备总成本使用蒙特卡洛算法模拟1000次后得到下图，不难看出大多分布在60000-80000区间范围内

• 对设备效益值使用蒙特卡洛算法模拟1000次后得到下图，不难看出大多分布在
  3000-6000区间范围内
  

对成本效益比率，我们首先要确定每个变量的概率分布，根据概率分布生成大量随机样本，再计算相对应的目标函数值，收集所有的目标函数值，进行统计分析，计算出均值和方差，经过处理得到：

Mean value of CER(e): 0.07217944998576989
Standard deviation of CER(e): 0.0036991256889698423

• 5、额外的设备
  a) 急求设备
  ●紧急呼吸设备（EBAs）：用于紧急救治伤者，挽救伤者生命
  b）深海设备
  ●Deep-sea Camera:使用可以承受高水压的摄像机，可以在深海中捕捉图像，帮助定位潜水器
  ●惯性导航系统（INS）：惯性导航系统利用加速度计和陀螺仪来测量和计算潜水器的位置、速度、方向等动态参数;它不依赖于外部信号，因此在水下环境中非常适用
  c) 浅水区域或水面设备
  ●卫星定位系统（如GPS接收机）：虽然GPS信号无法直接在水下使用，但潜水器在浮出水面时可以利用GPS进行精确定位,这些数据可以用来校准和验证其他导航系统的准确性

五、模拟搜索

1、问题分析
我们需要建立一个合适的模型来推荐设备的初始部署点和搜索模式，以尽量减少丢失潜水器的定位时间.
2、模型的准备
a)模型假设
●假设搜索区域为预测路径周围
●假设10min未收到潜艇发送的消息，开始搜索
●考虑到潜水器定位是一个动态过程，模型可能需要定期更新以适应不断变化的环境和情况
b)处理步骤
●数据收集与预处理
我们将随机生成潜水器与主船的坐标位置保存为数据集，样本含有x,y,z坐标值
●提取特征值
根据随机生成的坐标值模拟出主船和潜水器之间的航行路线，先假设主船和潜水器的航行 路线为短程，再逐次增加坐标点数，观察模拟出的路线，以主船和潜水器之间的距离和时间为特征值
●模型的选择和训练
我们将模拟数据集使用贝叶斯网络分类器进行训练。然后，通过可视化展示了主船和潜水器路径以及分类器的决策边界
●模型的评估和优化
需要综合考虑潜水器的搜索效率和定位时间以及变化的环境因素等指标对模型的准确度造成的影响，找到一个合适的平衡点

c)搜索模式
●圆形搜索（Circle search）：这是一种相对简单直接的方法，通常围绕一个固定的参考点进行一系列圆形搜索，逐步增加圆的直径.
●螺旋搜索（Spiral search）：这是一种在二维矩阵中按照螺旋路径遍历元素的方法;在螺旋搜索中，我们从矩阵的左上角元素开始，按照顺时针的方向依次访问矩阵中的元素，直到遍历完所有元素为止.
●网格搜索（Grid search）：这是一种穷举搜索方法，它通过遍历超参数的所有可能组合来寻找最优超参数.
下面我们采用网格搜索模式
3、模型的建立
●贝叶斯网络概述
贝叶斯网络(Bayesian network,BN)被称为因果概率网络，基础是贝叶斯公式;贝叶斯网络各节点变量的条件概率是在贝叶斯公式与统计学的基础上发展而来.
设$x_1,x_2,x_3,...,x_n$是变量X的1个互不相容的完备事件组，且$p(x_i)>0$,$Y$是$X$引发的事件，贝叶斯公式如下：
$P\left(x_i\mid Y\right)=\frac{P\left(Y\mid x_i\right)P\left(x_i\right)}{\sum _{i=1}^{n}P\left(x_i\right)P\left(Y\mid x_i\right)}$

• 对应到我们的模型中的话，搜索区域被划分为 $n$个网格，每个网格$i$中潜水器存在的先验概率为$P_i$， $x_i$表示潜水器位于网格 $i$， $Y$表示观测信息，$P\left(x_i\mid Y\right)$ 是给定观测信息$Y$ 后潜水器位于网格$i$的后验概率.

• ●初始贝叶斯网络航迹模型
  贝叶斯网络将每个变量视为节点，用有向边表示它们之间的条件概率依赖关系，如图
  
  
  由上图可知，给定当前时刻深度、加速度、升降率、 转弯率 ４个变量值，根据动态贝叶斯网络及各个变量的条件概率矩阵，即可实现下一时刻航迹预测和生成.

• 基于前文航迹变量介绍，深度变量可以直接获取， 速度、加速度、转弯率和升降率则通过推导获取：
  1）速度和加速度变量
  速度通过已知速度分量整合获得：
  $V(t)=\sqrt{V_{\mathrm{E}}^1\left(t\right)+V_{\mathrm{N}}^2\left(t\right)+V_{\mathrm{G}}^3\left(t\right)}$

对前后两个时刻做差可得加速度数值：
$A(t)=V(t+1)-V(t)$

2)升降率变量
通过对前后2个时刻深度数值做差可得：
$Ah(t)=L(t+1)-L(t)$

3)航向角和转弯率变量
将东向速度和北向速度做比，得到航向参考数据，设定北为轴正方向，东为轴正方向，将航向参考数据转换成坐标系下的航向角
$\begin{array}{rcl}\theta (t)& =& \frac{\arctan \left|\left(\frac{V_{{}_{\mathrm{E}}}\left(t\right)}{V_{{}_{\mathrm{N}}}\left(t\right)}\right)\right|\times 180^{\circ }}{\pi }\end{array}$

$\theta (t)$为潜水器的坐标系航向角

●针对潜水器升降率、转弯率和加速度，均可建立类似条件概率矩阵。在之后迭代过程中，通过调用变量条件概率矩阵，依照子节点取值概率分布生成预测值.
● 针对预测得到的潜水器航迹数据，进行可视化处理，建立三维坐标系。其中，Ｚ轴取值地面垂直向上为 正，Ｘ轴取值向东为正，Ｙ轴取值向北为正。在 Ｘ-Ｙ平面 上，设定航向初始值，将预测出来的转弯率变化量不断 叠加，迭代生成航向新值。基于运动学方程，生成潜水器在 Ｘ轴、Ｙ轴、Ｚ轴的实时位置:

$$\begin{gathered} \begin{array}{rcl}x(t+1)& =& x(t)+V(t)\cos \theta (t)\\ y(t+1)& =& y(t)+V(t)\sin \theta (t)\end{array} \\ \theta (t)={\sum }_{t=1}^{t}RATE(t)+{\theta }_0 \\ z(t)=L(t) \end{gathered}$$

作者：爱建模的鹿鹿
链接：https://www.zhihu.com/question/642461787/answer/3384560591
来源：知乎
著作权归作者所有。商业转载请联系作者获得授权，非商业转载请注明出处。

import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D

# 模拟潜水器在三维空间中的运动轨迹
t = np.linspace(0, 20, 1000)
x = t * np.sin(t)
y = t * np.cos(t)
z = np.log(t + 1)

# 模拟另一潜水器为对比
x2 = t * np.sin(t + 1)
y2 = t * np.cos(t + 1)
z2 = np.log(t + 0.5)

fig = plt.figure(figsize=(14, 10))
ax = fig.add_subplot(111, projection='3d')

# 绘制潜水器的运动轨迹
ax.plot(x, y, z, label='Submersible 1 Trajectory', color='blue', linewidth=2)
ax.plot(x2, y2, z2, label='Submersible 2 Trajectory', color='red', linestyle='--', linewidth=2)

# 标记起点和终点
ax.scatter(x[0], y[0], z[0], color='green', s=100, label='Start Point', edgecolor='black')
ax.scatter(x[-1], y[-1], z[-1], color='orange', s=100, label='End Point', edgecolor='black')

# 设置图形属性
ax.set_title('3D Trajectories of Submersibles in the Sea', fontsize=20)
ax.set_xlabel('X Axis', fontsize=14)
ax.set_ylabel('Y Axis', fontsize=14)
ax.set_zlabel('Z Axis (Depth)', fontsize=14)
ax.legend(fontsize=12)

plt.show()

●假设模型
我们生成简单的模拟数据来观察主船和潜水器的路径，再构建特征和标签，0表示主船，1表示潜水器，然后我们使用贝叶斯网络分类进行训练，构建网格来绘制决策边界，得到主船和潜水器的路径图，这个只是我们简单的模拟
ship_path： [[1, 2], [2, 3], [3, 4], [4, 5], [5, 6]]
submarine_path： [[1, 1], [2, 2], [3, 3], [4, 4], [5, 5]]

接下来我们采用随机生成它将生成 100 个随机的点数对来模拟主船和潜水的坐标，因为是模拟数据，所以和实际会有较大偏差，如下图

我们接下来根据坐标点随机出现的位置进行概率模拟，继续使用贝叶斯网络分类进行训练，生成网格预测数据，如下图

4、模型的误差分析与改进
a)模型分析
可以看出主船和潜水器出现的概率分布，因为是模拟数据而且是随机分布，所以造成的偏差是非常大的，不过我们这里讨论的是建立一个合适的模型来推荐设备的初始部署点和搜索模式，以尽量减少丢失潜水器的定位时间，如果我们有历来探险者的航线路线，就可以根据航行路线过程中的坐标信息进行多次训练，来达到较快搜索到丢失潜水器.
b)模型改进与建议
数据集需要采集历年探险者的航行坐标以便于增加置信度，不断提取新的特征值，比如潜水器与主船之间的相对速度，路径曲率等，使用交叉验证找到最优的参数组合，比如网格搜索，随即搜索，贝叶斯优化，训练多组模型，对模型进行集成，降低方差，提升准确度.

六、拓展模型

1、问题分析
在这个问题上，我们需要将模型拓展到其他旅游目的地，比如加勒比海等，也就是说需要将模型进行环境迁移，同时也要考虑潜水器运动模型;对于存在多个潜水艇的情况，则需要使用MARL算法来对多个潜水艇实现追踪与预测，救援首先要从概率最大的救援区域开始寻找.
2、模型的准备
a)模型假设
●假设潜水器与主船所在区域天气等不确定因素忽略
●假设多个潜水器的航行路线不互相冲突
●假设潜水器的燃料充足,通信设备正常
b)处理步骤
●数据收集与预处理
先收集加勒比海地区的主船和潜水器路径数据,并进行预处理和清洗.
●提取特征值
由于需要考虑到海域的变更,我们需要加入更多的环境因素,比如海底地形的复杂性,海流和风向的影响,可能存在的潜水点等
●模型的选择和训练
针对多个潜水器在同一区域移动的情况,一个常用的算法是多智能体强化训练学习中的Deep -Q -Network算法(MADQN),MADQN常用于解决多个智能体协作和竞争问题
●模型的评估和优化
针对加勒比海地区的特殊情况,可能要调整模型的超参数或结构以达到更好的性能

3、模型的建立
●面对上述场景和需求,我们决定分别解决地点迁移和多智能体强化学习这两大方向的问题

●地点迁移
地点迁移会导致主船和潜水器航行过程中的不确定因素增加,其中影响最大的因素便是潜水器下潜过程中所到的洋流干扰,因此我们在第一问中建立的动力学预测模型将会存在不确定性,针对潜水器存在的洋流干扰,我们需要建立洋流影响下的潜水器模型以此来提高潜水器的控制性能.

a)洋流影响下的潜水器模型
假设洋流速度恒定且无旋,那我们现在给出以下模型:

$\dot{P}=JV$
$M_{RB}\dot{V}+M_A{\dot{V}}_r+C_{RB}(V)V+C_A(V_r)V_r+D(V_r)V_r+G(P)=F+d$
将受洋流影响的潜水器模型整理得到如下潜水器动力学模型：

$M_{RB}\dot{V}+M_A{\dot{V}}_r+C_{RB}(V)V+C_A(V_r)V_r+D(V_r)V_r+G(P)=F+d+d_c$

现在我们可以确定潜水器所受的洋流干扰:
$d_c=M_{RB}\dot{V}+M_A{\dot{V}}_r+C_{RB}(V)V+C_A(V_r)V_r+D(V_r)V_r+G(P)-F-d$

●多智能体强化学习模型
我们首先将加勒比海地区的环境建模为一个多智能体环境,把每一个潜水器当作一个智能体,根据环境状态选择动作,并根据奖励信号来调整策略初步确定使用 Deep-Q-Network算法进行训练 ,通过模拟多个潜水器在环境中移动，并根据环境返回的奖励信号来训练智能体的行为.每个智能体都有自己的 DQN 模型，并独立地选择动作和学习策略。在每个 episode 中，智能体执行动作并观察环境的反馈，然后根据反馈更新自己的 Q 函数

在强化学习过程中，一个“episode”通常指的是从环境的初始状态开始，直到达到终止状态为止 的一系列动作和状态转换过程.在每个episode中，智能体与环境进行交互，采取动作，观察环境的状态和奖励，并尝试通过学习来提高其行为策略;在每个episode结束时，输出该episode的序号以及智能体在该episode中获得的奖励.

1.绘制累积奖励曲线：在每个训练episode结束时，记录智能体获得的累积奖励，并绘制累积奖励随训练轮次变化的曲线图.帮助我们观察智能体在训练过程中奖励的变化趋势.
2.绘制奖励分布直方图：将智能体获得的奖励进行统计，并绘制奖励分布的直方图,帮助我们了解智能体在训练过程中不同奖励值的分布情况.

4、模型的分析与改进
●分析:
○根据累积奖励曲线和奖励分布直方图,我们可以知道,随着训练次数的增加,对每个潜水器的运动姿态的预测性就会更加准确,观测上图会发现整体训练效果比较平稳,即使面对不同的环境也能保持良好的一致性.
○建议选取更多的数据集进行训练,增加迭代次数来达到更好的效果.
●改进:
○1在环境模型、智能体设计和学习算法确定后，使用历史数据对智能体进行训练,并不断优化智能体的策略,以最大化累积奖励.
○ 通过在仿真环境中进行测试和验证,评估训练好的智能体在实际场景中的性能，并根据需要进行调整和优化.

八、敏感性分析

分析:
●基于Kalman滤波算法改进后的动态预测模型模拟潜水器位置在实际应用中需要实时更新,总体来说,模型预测与事先给定路径基本吻合,波动性较小

●使用 Deep-Q-Learning算法来处理多个潜水器的移动问题，考虑到了潜水器之间的协作和竞争关系，从而更好地优化整体的搜寻效率和定位准确度,如果继续优化参数,迭代一定次数之后,模型的准确性会更高

小结

以上是这次参加美赛B题的一些分析,尚未完全贴出,整体来说还是较为局限的,毕竟时间也有限,下一篇再分享一下代码部分.