火炮射击中的数学建模问题

代码

# 问题 1
import math
import matplotlib.pyplot as plt
# 炮弹飞行轨迹模型
def calculate_trajectory(m, theta, phi, v0, k):
 g = 9.8 # 重力加速度
 dt = 0.01 # 时间间隔
 x = [0] # x 坐标列表
 y = [0] # y 坐标列表
 vx = [v0 * math.cos(math.radians(theta)) * math.cos(math.radians(phi))] # x 轴速度列表
 vy = [v0 * math.sin(math.radians(theta))] # y 轴速度列表
 max_height = 0 # 最大高度
 max_height_x = 0 # 最大高度点的水平位置
 while y[-1] >= 0:
        # 计算下一个时间点的速度和坐标
     v = math.sqrt(vx[-1]**2 + vy[-1]**2) # 当前速度
     ax = -k * v * vx[-1] / m # x 轴加速度
     ay = -g - k * v * vy[-1] / m # y 轴加速度
     x.append(x[-1] + vx[-1] * dt) # 计算下一个时间点的 x 坐标
     y.append(y[-1] + vy[-1] * dt) # 计算下一个时间点的 y 坐标
     vx.append(vx[-1] + ax * dt) # 计算下一个时间点的 x 轴速度
     vy.append(vy[-1] + ay * dt) # 计算下一个时间点的 y 轴速度
    if y[-1] > max_height:
        max_height = y[-1]
        max_height_x = x[-1]
 return x, y, max_height_x
# 计算轨迹
m = 20
theta = 28
phi = 35
k = 0.0042
dt = 0.01
v0_1 = 509
v0_2 = 1200
x1, y1, max_height_x1 = calculate_trajectory(m, theta, phi, v0_1, k)
x2, y2, max_height_x2 = calculate_trajectory(m, theta, phi, v0_2, k)
# 计算最大高度
max_height_1 = max(y1)
max_height_2 = max(y2)
# 计算落地时间
landing_time_1 = len(x1) * dt
landing_time_2 = len(x2) * dt
# 计算落点位置
landing_point_1 = (x1[-1], y1[-1])
landing_point_2 = (x2[-1], y2[-1])
# 计算最大高度的水平距离
horizontal_distance_1 = max_height_x1 - x1[0]
horizontal_distance_2 = max_height_x2 - x2[0]
# 输出最大高度、落地时间、落点位置和最大高度的水平距离
print("炮弹 1 的最大高度:", max_height_1, "米")
print("炮弹 1 的落地时间:", landing_time_1, "秒")
print("炮弹 1 的落点位置:", landing_point_1)
print("炮弹 1 最大高度的水平距离:", horizontal_distance_1, "米")
print("炮弹 2 的最大高度:", max_height_2, "米")
print("炮弹 2 的落地时间:", landing_time_2, "秒")
print("炮弹 2 的落点位置:", landing_point_2)
print("炮弹 2 最大高度的水平距离:", horizontal_distance_2, "米")
# 绘制轨迹
plt.plot(x1, y1, label="v0 = 509 m/s")
plt.plot(x2, y2, label="v0 = 1200 m/s")
plt.xlabel("x")
plt.ylabel("y")
plt.title("Projectile Trajectory")
plt.legend()
# 标记最大高度
plt.annotate(f"Max Height: {max_height_1:.2f} m", xy=(max_height_x1, max_height_1), 
xytext=(max_height_x1 + 1000, max_height_1 + 100),
 arrowprops=dict(facecolor='black', arrowstyle='->'))
plt.annotate(f"Max Height: {max_height_2:.2f} m", xy=(max_height_x2, max_height_2), 
xytext=(max_height_x2 - 2000, max_height_2 + 100),
 arrowprops=dict(facecolor='black', arrowstyle='->'))
# 标记落地时间和落点位置
plt.annotate(f"Landing Time: {landing_time_1:.2f} s", xy=(x1[-1], y1[-1]), xytext=(x1[-1] + 3000, 
y1[-1] + 500),
 arrowprops=dict(facecolor='black', arrowstyle='->'))
plt.annotate(f"Landing Time: {landing_time_2:.2f} s", xy=(x2[-1], y2[-1]), xytext=(x2[-1] - 2000, 
y2[-1] + 1000),
 arrowprops=dict(facecolor='black', arrowstyle='->'))
plt.annotate(f"Landing Point", xy=(x1[-1], y1[-1]), xytext=(x1[-1] + 1000, y1[-1] + 3000),
 arrowprops=dict(facecolor='black', arrowstyle='->'))
plt.annotate(f"Landing Point", xy=(x2[-1], y2[-1]), xytext=(x2[-1] - 2000, y2[-1] + 2000),
 arrowprops=dict(facecolor='black', arrowstyle='->'))
# 标记最大高度的水平距离
plt.annotate(f"Horizontal Distance: {horizontal_distance_1:.2f} m", xy=(max_height_x1, 
max_height_1), xytext=(max_height_x1 - 2000, max_height_1 - 2000),
 arrowprops=dict(facecolor='black', arrowstyle='->'))
plt.annotate(f"Horizontal Distance: {horizontal_distance_2:.2f} m", xy=(max_height_x2, 
max_height_2), xytext=(max_height_x2 + 2000, max_height_2 - 2000),
 arrowprops=dict(facecolor='black', arrowstyle='->'))
plt.grid(True)
plt.savefig('问题 1.png')
plt.show()
# 问题 2
import pandas as pd
import numpy as np
from scipy.optimize import least_squares
# 读取雷达数据
data = pd.read_excel("附件 1.xlsx")
# 提取距离、仰角和方位角，并转化为弧度制
distance = data["目标距离（米）"].values
height_angle = np.deg2rad(data["目标仰角（度）"].values)
azimuth_angle = np.deg2rad(data["目标方位角（度）"].values)
# 计算各时刻炮弹的空间坐标
x = distance * np.cos(height_angle) * np.sin(azimuth_angle)
y = distance * np.cos(height_angle) * np.cos(azimuth_angle)
z = distance * np.sin(height_angle)
# 用抛物线方程对炮弹轨迹进行拟合
# 设抛物线方程为 z = a * (x - x0)**2 + b * (y - y0)**2 + z0，通过最小二乘法求解参数 a, b, x0, 

def func(params):
 a, b, x0, y0, z0 = params
 return a * (x - x0)**2 + b * (y - y0)**2 + z0 - z
initial_params = [1, 1, 0, 0, 0]
result = least_squares(func, initial_params)
a, b, x0, y0, z0 = result.x
print(f"Fire cannon might be located at ({x0:.2f}, {y0:.2f}).")
# 炮弹落点的空间坐标为抛物线的顶点，即 (x0, y0, z0)
print(f"Shell might fall at ({x0:.2f}, {y0:.2f}, {z0:.2f}).")
import matplotlib.pyplot as plt
# 计算各时刻炮弹的速度
velocity = np.sqrt(np.diff(x)**2 + np.diff(y)**2 + np.diff(z)**2)
# 假设每次观测的时间间隔为 delta_t
delta_t = 1.0 # 换算成合适的单位
# 计算速度
velocity = velocity / delta_t
# 使用第一个速度数据填充，使得速度数据和位置数据一样多
velocity = np.insert(velocity, 0, velocity[0])
# 绘制速度随时间的变化曲线
plt.plot(np.arange(len(velocity)), velocity)
plt.xlabel('Time')
plt.ylabel('Velocity')
plt.title('Velocity of the shell over time')
plt.show()
# 问题 3
import pandas as pd
import numpy as np
# 读取附件 2 数据
df_attach2 = pd.read_excel('附件 2.xlsx')
x_attach2 = df_attach2['X 坐标（米）']
y_attach2 = df_attach2['Y 坐标（米）']
# 火炮目标点
target_point = (40, 60)
# 计算诸元误差的分布参数
mean_xe = np.mean(x_attach2 - target_point[0])
mean_ye = np.mean(y_attach2 - target_point[1])
std_xe = np.std(x_attach2 - target_point[0])
std_ye = np.std(y_attach2 - target_point[1])
print("诸元误差的分布参数：")
print("mean_xe =", mean_xe)
print("mean_ye =", mean_ye)
print("std_xe =", std_xe)
print("std_ye =", std_ye)
# 读取附件 3 数据
df_attach3 = pd.read_excel('附件 3.xlsx')
x_attach3 = df_attach3['X 坐标（米）']
y_attach3 = df_attach3['Y 坐标（米）']
# 火炮目标点
target_point_1 = (40, 60)
target_point_2 = (43, 70)
# 计算散布误差参数和两次任务的诸元误差
std_xb = np.std(x_attach3 - target_point_1[0])
std_yb = np.std(y_attach3 - target_point_1[1])
mean_xe_task2 = np.mean(x_attach3 - target_point_2[0])
mean_ye_task2 = np.mean(y_attach3 - target_point_2[1])
mean_xe_task1 = np.mean(x_attach3 - target_point_1[0])
mean_ye_task1 = np.mean(y_attach3 - target_point_1[1])
print("\n 散布误差参数和两次任务的诸元误差：")
print("std_xb =", std_xb)
print("std_yb =", std_yb)
print("mean_xe_task2 =", mean_xe_task2)
print("mean_ye_task2 =", mean_ye_task2)
print("mean_xe_task1 =", mean_xe_task1)
print("mean_ye_task1 =", mean_ye_task1)
# 问题 4
import pandas as pd
import numpy as np
from sklearn.model_selection import train_test_split
from sklearn.neighbors import KNeighborsClassifier
# 从 Excel 读取数据
target_points_data = pd.read_excel('附件 4.xlsx', sheet_name='目标点')
error_params_data = pd.read_excel('附件 4.xlsx', sheet_name='射击误差参数')
landing_positions_data = pd.read_excel('附件 4.xlsx', sheet_name='落点位置')
# 提取目标点数据
target_points = target_points_data[['X 坐标（米）', 'Y 坐标（米）']].values
# 提取落点位置数据
landing_positions = landing_positions_data[['X 坐标（米）', 'Y 坐标（米）']].values
# 建立目标点和落点位置之间的映射
target_to_landing_map = {tuple(target): landing for target, landing in zip(target_points, 
landing_positions)}
# 提取火炮编号和设计误差参数数据
cannon_ids = error_params_data['火炮编号'].values
error_params = error_params_data[['诸元误差 X 坐标（米）', '诸元误差 Y 坐标（米）','散布误
差方位向（X）参数（米)','散布误差距离向（Y）参数（米）'].values
# 创建训练数据
X_train, X_test, y_train, y_test = train_test_split(target_points, cannon_ids, test_size=0.2, 
random_state=42)
# 建立 K 最近邻分类模型
knn = KNeighborsClassifier(n_neighbors=3)
knn.fit(X_train, y_train)
# 使用模型进行预测
predicted_cannon_ids = knn.predict(X_test)
# 打印预测结果
print("预测结果:")
for i, pred in enumerate(predicted_cannon_ids):
 print("目标点 ({:.3f}, {:.3f}) 的预测火炮编号为 {}".format(X_test[i, 0], X_test[i, 1], pred))
# 自行生成数据进行验证
new_data = np.array([
 [-110.0, -50.0], # 示例新的目标点数据
 [90.0, 90.0] # 示例新的目标点数据
])
# 打印验证结果
print("\n 验证结果:")
for i, data in enumerate(new_data):
 if tuple(data) in target_to_landing_map:
 pred = knn.predict([data])[0]
 landing = target_to_landing_map[tuple(data)]
 print(" 目 标 点 ({:.3f}, {:.3f}) 的 预 测 火 炮 编 号 为 {} ， 落 点 位 置 为 ({:.3f}, 
{:.3f})".format(data[0], data[1], pred, landing[0], landing[1]))
 else:
 print("目标点 ({:.3f}, {:.3f}) 没有找到对应的落点位置".format(data[0], data[1]))