2023年全国大学生数学建模比赛B题（名次未知）
import numpy as np
import pandas as pd

data = pd.read_excel('result1.xlsx')
#初始化
D0 = 70             #测线距中心线距离为0时海水深度
a = np.deg2rad(1.5) #海底坡面与水平面夹角
l0 = D0 / np.tan(a) #中心位置与理想岸边距离
print('中心位置与理想岸边距离l0 = ',l0)
ct = np.deg2rad(120) #多波束换能器的开角
a1 = (np.pi - ct)/2 - a #多波束换能器与坡面左侧夹角
a2 = (np.pi - ct)/2 + a #多波束换能器与坡面右侧夹角
Sea_depth = []          #海水深度
Covering_width = []     #覆盖宽度
Overlap_rate = []       #与前一条测线的重叠率
k = 0
for d in range(-800,800+1,200):
    D = (l0 - d)*np.tan(a) #不同测线距离海水深度
    Sea_depth.append(D)
    print('d = %d\t, D = ' % d, D,end = '\t')
    # 覆盖宽度
    W = D*np.sin(ct/2)*(1/np.sin(a1) + 1/np.sin(a2))
    print('W = ',W,end='\t')
    Covering_width.append(W)
    if(d == -800): #此测线距离不计算重叠率
        print()
    else:
        #与前一条测线重叠率
        rate = (1 - (200*np.cos(ct/2) / (Covering_width[k]*np.sin(a1)))) * 100
        Overlap_rate.append(rate)
        k+=1
        print('rate = ', rate)
#存入各个数据到表格中
for i in range(1,10): #存入海水深度
    data.iloc[0,i] = Sea_depth[i-1]
for i in range(1,10): #存入覆盖宽度
    data.iloc[1,i] = Covering_width[i-1]
for i in range(1,9): #存入重叠率
    data.iloc[2,i+1] = Overlap_rate[i-1]
data.to_excel('result1.xlsx',index=None)
支撑材料/相关代码/Question1.py
tianbosong/CUMCM2023B
import numpy as np
import pandas as pd

data = pd.read_csv('result2.csv',encoding='gbk')
#初始化
Sea_mile = 1852         #单位转化
a0 = np.deg2rad(1.5)    #海底坡面与水平面夹角
D0 = 120                #测线距中心线距离为0时海水深度
ct = np.deg2rad(120)    #多波束换能器的开角
k = 0

for b0 in range(0,315+1,45):
    Covering_width = []        #覆盖宽度
    print('当β = %d时' % b0,end=' ')
    Beta_Angle = np.deg2rad(b0) #贝塔角
    a = np.sin(Beta_Angle) * a0 #真实投影角
    print('a= ',np.rad2deg(a))
    a1 = (np.pi - ct) / 2 - a   #多波束换能器与坡面左侧夹角
    a2 = (np.pi - ct) / 2 + a   #多波束换能器与坡面右侧夹角
    for i in np.arange(0,2.2,0.3):
        x = i * Sea_mile        #测量船距海域中心点处的距离，单位米
        #覆盖宽度计算
        if Beta_Angle==0:       #特殊情况
            W = 2*(D0+x*np.tan(a0))*np.tan(ct/2)
            Covering_width.append(W)
            print(W)
        else:
            W = (D0 + x*np.tan(a)/np.tan(Beta_Angle))* np.sin(ct / 2) * (1 / np.sin(a1) + 1 / np.sin(a2))
            Covering_width.append(W)
            print(W)
    for i in range(len(Covering_width)):
        data.iloc[k, i + 1] = Covering_width[i]
    k += 1
data.to_csv('result2.csv', index=None,encoding='gbk')


支撑材料/相关代码/Question2.py
tianbosong/CUMCM2023B
import numpy as np
import pandas as pd
data = pd.DataFrame(columns=['测量带条数','间隔距离','最小覆盖率','最大覆盖率'])
#初始化
Sea_mile = 1852             #单位转化
D0 = 110                    #测线距中心线距离为0时海水深度
Length = 2 * Sea_mile       #矩形海域长
Width = 4 * Sea_mile        #矩形海域宽
a0 = np.deg2rad(1.5)        #海底坡面
ct = np.deg2rad(120)        #多波束换能器的开角
Xmin = 0                    #测量船只移动最小距离
Xmax = Width                #测量船只移动最大距离
ans = 0                     #统计满足条件的个数
#计算最大覆盖宽度
dt = Length * np.tan(a0)    #中心海域与最深处高度差
H = D0 + dt                 #最深处深度
Wmax = 2*(H-Xmin*np.tan(a0))*np.tan(ct/2)
print('最大覆盖宽度：',Wmax)

#计算最小覆盖宽度
l0 = D0 / np.tan(a0)            #理想岸边距离
h = (l0 - Length)*np.tan(a0)    #最浅处海水深度
Wmin = 2*(H-Xmax*np.tan(a0))*np.tan(ct/2)
print('最小覆盖宽度：',Wmin)

#两侧测线距边界距离
x0 = Wmin / 2
print('北侧第一条测线距离：',x0)

for n in range(2,300):
    d = (Length - 2*x0) / n     #排除边界的测线间隔
    #print('测量带条数n = %d, 间隔d = %d' % (n+2,d),end='\t')
    rate_max = 1-d/Wmax
    rate_min = 1-d/Wmin
    #print('覆盖率范围：(%f,%f)' % (rate_min,rate_max))
    data = data.append({'测量带条数': n+2, '间隔距离': d,
                        '最小覆盖率': rate_min, '最大覆盖率': rate_max}, ignore_index=True)
    if rate_min > 0.1 and rate_max < 0.2:
        ans += 1
        #print('符合要求')
data.to_excel('result3_1.xlsx', index=None)
print('共有%d种情况符合要求' % ans)
支撑材料/相关代码/Question3_1.py
tianbosong/CUMCM2023B
import numpy as np
import pandas as pd
data = pd.DataFrame(columns=['当前测线位置序号','海水深度','距上一条测线距离','覆盖率'])
#初始化
Sea_mile = 1852             #单位转化
D0 = 110                    #测线距中心线距离为0时海水深度
Length = 2 * Sea_mile       #矩形海域长
Width = 4 * Sea_mile        #矩形海域宽
a0 = np.deg2rad(1.5)        #海底坡面
ct = np.deg2rad(120)        #多波束换能器的开角
a1 = (np.pi - ct) / 2 - a0  #多波束换能器与坡面左侧夹角
a2 = (np.pi - ct) / 2 + a0  #多波束换能器与坡面右侧夹角
target_rate = 0.1           #最优覆盖率

#计算起始测线位置
dt = Length * np.tan(a0) #中心海域与最深处高度差
H = D0 + dt              #最深处深度
s0 = H / np.tan(a0+a1)   #起始测线距西边界距离

#计算末位置测线
l0 = D0 / np.tan(a0)            #理想岸边距离
h = (l0 - Length)*np.tan(a0)    #最浅处海水深度
sn = Width - h * np.tan(ct/2)   #末位置测线距西边界距离

s = []          #各个测线间距
s.append(s0)
i = 2
cnt = 0
data = data.append({'当前测线位置序号': 1, '海水深度': np.tan(a0+a1)*s0,
                                    '距上一条测线距离': s0,'覆盖率': '——'}, ignore_index=True)
while sum(s) < sn:
    rate_temp = 0
    sx_temp = 0
    D_temp = 0
    cnt = 0
    for sx in range(20,600):
        #print('第%d个测线，迭代%d次' % (i,cnt))
        cnt += 1
        #print('距离上一条航线 %d 米' % sx)
        D = H - (sx + sum(s)) * np.tan(a0)  # 当前海水深度
        #print('\t海水深度:', D)
        W = D * np.sin(ct / 2) * (1 / np.sin(a1) + 1 / np.sin(a2))
        #print('\t测量宽度:', W)
        # ans = ((1 - rate) * W * np.sin(a1)) / np.cos(ct / 2)
        rate = 1 - (sx * np.cos(ct / 2) / (W * np.sin(a1)))
        #print('\t覆盖率：', rate)
        #print('\t当前测量到：', sum(s))
        if rate < 0.12:
            if rate < target_rate:
                data = data.append({'当前测线位置序号': i, '海水深度': D_temp,
                                    '距上一条测线距离': sx_temp,'覆盖率': rate_temp}, ignore_index=True)
                s.append(sx_temp)  #最优
                break
            else:
                D_temp = D
                sx_temp = sx
                rate_temp = rate
                ...
            ...
        ...
    i += 1
    ...
data.to_excel('result3_2.xlsx', index=None)
print('总测线长度为：',Length*len(s))



支撑材料/相关代码/Question3_2.py
tianbosong/CUMCM2023B
import numpy as np
import matplotlib.pyplot as plt

# 定义常量
Sea_mile = 1852         # 海里（m）
Length = 2 * Sea_mile   # 海域长度（m）
Width = 4 * Sea_mile    # 海域宽度（m）

# 问题2得到的各间距距离
s = [359, 546, 507, 470, 436, 404,
     375, 348, 323, 299, 278, 257,
     239, 221, 205, 191, 177, 164,
     152, 141, 131, 121, 113, 104,
     97, 90, 83, 77, 72, 67, 62,
     57, 53, 49, 46, 42, 39]

# 初始化测线列表
x = [s[0]]      # 测线的横坐标列表（m）
y1 = [0]        # 测线的起点纵坐标列表（m）
y2 = [Length]   # 测线的终点纵坐标列表（m）

for i in range(len(s)):
    x.append(x[-1] + s[i])                      # 添加下一条测线的横坐标
    y1.append(Length if len(x) % 2 == 0 else 0)  # 添加下一条测线的起点纵坐标
    y2.append(0 if len(x) % 2 == 0 else Length)  # 添加下一条测线的终点纵坐标

# 绘制示意图
plt.figure(figsize=(8,4))               # 创建画布
plt.plot(x, y1, 'bo', label='start')    # 绘制起点
plt.plot(x, y2, 'ro', label='end')      # 绘制终点
plt.plot([x, x], [y1, y2], 'k--')       # 绘制测线
plt.xlabel('X (m)')                     # 设置横轴标签
plt.ylabel('Y (m)')                     # 设置纵轴标签
plt.title('Line diagram')               # 设置标题
plt.legend()                            # 显示图例
plt.show()                              # 显示图像

支撑材料/相关代码/Question3_plot.py
tianbosong/CUMCM2023B
import numpy as np  
  
class Particle:  
    def __init__(self, dim, minx, maxx):  
        self.position = np.random.uniform(low=minx, high=maxx, size=dim)  
        self.velocity = np.random.uniform(low=-0.1, high=0.1, size=dim)  
        self.best_position = np.copy(self.position)  
        self.best_score = -np.inf  
  
class PSO:  
    def __init__(self, dim, n_particles, bounds):  
        self.dim = dim  
        self.n_particles = n_particles  
        self.bounds = bounds  
        self.particles = [Particle(dim, *bounds) for _ in range(n_particles)]  
          
    def optimize(self, func, n_iter):  
        global_best_score = -np.inf  
        global_best_position = None  
  
        for i in range(n_iter):  
            for particle in self.particles:  
                score = func(particle.position)  
  
                if score > particle.best_score:  
                    particle.best_score = score  
                    particle.best_position = particle.position  
  
                if score > global_best_score:  
                    global_best_score = score  
                    global_best_position = particle.position  
  
            for particle in self.particles:  
                particle.velocity += 2 * np.random.rand(self.dim) * (particle.best_position - particle.position) +  2 * np.random.rand(self.dim) * (global_best_position - particle.position)  
                particle.position += particle.velocity  
  
                for i in range(self.dim):  
                    if particle.position[i] < self.bounds[0]:  
                        particle.position[i] = self.bounds[0]  
                    elif particle.position[i] > self.bounds[1]:  
                        particle.position[i] = self.bounds[1]  
  
        return global_best_position  
def fitness(position):  
    total_length = np.sum(position)  
    missed_area_percentage = np.sum(data[position == 0]) / np.sum(data)  
    overlap_length = np.sum(position[position > 1] - 1)  
      
    # 希望最小化这三个指标  
    return -(total_length + missed_area_percentage + overlap_length)  
  
pso = PSO(dim=20, n_particles=30, bounds=(0, 1))  
best_position = pso.optimize(fitness, n_iter=1000)  

支撑材料/相关代码/Question4.py
tianbosong/CUMCM2023B