【数学建模2025】2025年数学建模国赛A题比赛记录
·
2025年全国大学生数学建模竞赛A题比赛记录
欢迎大家提出批评，提出建议，分享见解，交流经验
本人首次参赛，独狼，兼顾建模，编程，论文
本人比赛大体流程安排如下：
9月4日 18:00 发题后，读题选题，
9月5日 第一问第二问求解
9月6日第三问第四问求解，论文完成初稿
9月7日第五问求解，论文完成,晚20:00提交
首先分析题目
全文要解决的核心问题是寻找无人机参数的优化策略方案，使无人机抛掷的烟幕弹干扰导弹的时间尽可能长。
本人前五问答案如下，大家可以根据我的答案判断这篇文章有没有你参考的价值：
第一问1.455s
第二问4.479s
第三问6.042s
第四问25.9s
第五问17.8s(一定是错误答案，实际答案我没算出来)
本人主要解题参考了最优控制理论，遗传算法与约束规划
接下来是我的解题思路与源代码
导弹的飞行方向指向假目标（即原点），假目标旁边有个真目标
我们的目标是通过无人机扔的烟幕弹来遮蔽真目标，不让导弹发现真目标
本人遮蔽思路如下
将导弹与真目标间的连线看作导弹的视线，如果这一时刻这条线经过烟幕，那么判断这一时刻遮蔽有效
这个思路是全文的核心遮蔽有效判断逻辑
第一问
第一问算是定下来你对整个问题的思考方向，比较重要
进行求解时首先将真目标进行采样，因为题目里给了真目标的参数，所以肯定不能视作质点来计算，但也不能全计算，计算量太大也算不出来答案，我们选择把圆柱抽象化成一堆点进行采样，我采样了圆柱上下底面圆周与数条高。对于有效遮蔽时长策略的制定，我们将导弹轨迹与烟幕干扰弹爆弹后的烟幕轨迹抽象为数学模型，导弹轨迹实际上是一条由导弹原点直指假目标的直线，烟幕干扰弹的抛掷实际上是一次平抛运动，而烟幕干扰弹爆弹后的烟云轨迹是一条竖直向下的球状轨迹。而对于有效遮蔽时长的策略，我们将导弹与采样点间的连线定义为导弹视线方程，若导弹视线在这一时刻穿过了烟幕云轨迹，我们认为此刻烟幕的遮蔽是有效的，基于此我们可建立出问题一的求解方程，固定计算
本人代码如下
import numpy as np
import math
def calculate_problem1_maximum_precision():
g_val = 9.8
v_missile = 300.0
v_sink = 3.0
r_cloud = 10.0
t_effective = 20.0
p_fy1 = np.array([17800.0, 0.0, 1800.0])
p_m1 = np.array([20000.0, 0.0, 2000.0])
p_fake = np.array([0.0, 0.0, 0.0])
p_target = np.array([0.0, 200.0, 0.0])
vec_missile = p_fake - p_m1
dist_missile = np.linalg.norm(vec_missile)
u_missile = vec_missile / dist_missile
print(f"导弹单位方向向量: [{u_missile[0]:.10f}, {u_missile[1]:.10f}, {u_missile[2]:.10f}]")
t_drop = 1.5
t_delay = 3.6
v_uav = 120.0
p_drop = p_fy1 + np.array([-v_uav * t_drop, 0.0, 0.0])
dx = -v_uav * t_delay
dz = -0.5 * g_val * t_delay ** 2
p_burst = p_drop + np.array([dx, 0.0, dz])
t_burst = t_drop + t_delay
print(f"起爆点: ({p_burst[0]:.10f}, {p_burst[1]:.10f}, {p_burst[2]:.10f})")
print(f"起爆时间: {t_burst:.10f} s")
r_cyl = 7.0
h_cyl = 10.0
points_cyl = []
points_cyl.append(p_target.copy())
angles = np.linspace(0, 2 * np.pi, 72, endpoint=False)
heights = np.linspace(-h_cyl / 2, h_cyl / 2, 21)
for h in heights:
for a in angles:
x = p_target[0] + r_cyl * np.cos(a)
y = p_target[1] + r_cyl * np.sin(a)
z = p_target[2] + h
points_cyl.append(np.array([x, y, z]))
print(f"圆柱采样点总数: {len(points_cyl)}")
def check_intersection(p1, p2, center, radius):
d_vec = p2 - p1
to_center = p1 - center
a_val = np.dot(d_vec, d_vec)
b_val = 2 * np.dot(d_vec, to_center)
c_val = np.dot(to_center, to_center) - radius ** 2
disc = b_val ** 2 - 4 * a_val * c_val
if disc < 0:
return False
sqrt_disc = math.sqrt(disc)
t1 = (-b_val + sqrt_disc) / (2 * a_val)
t2 = (-b_val - sqrt_disc) / (2 * a_val)
return (0 <= t1 <= 1) or (0 <= t2 <= 1)
dt = 0.001
total_shield = 0.0
t_start = t_burst
t_end = t_burst + t_effective
print(f"时间窗口: [{t_start:.3f}s, {t_end:.3f}s]")
print(f"时间步长: {dt}s")
print(f"预计计算点数: {int((t_end - t_start) / dt)}")
all_points = []
intervals = []
current_start = None
count = 0
t_current = t_start
while t_current <= t_end:
pos_missile = p_m1 + v_missile * t_current * u_missile
sink_dist = v_sink * (t_current - t_burst)
cloud_center = p_burst + np.array([0.0, 0.0, -sink_dist])
shielded = False
for pt in points_cyl:
if check_intersection(pos_missile, pt, cloud_center, r_cloud):
shielded = True
break
if shielded and current_start is None:
current_start = t_current
elif not shielded and current_start is not None:
intervals.append((current_start, t_current - dt))
current_start = None
if shielded:
total_shield += dt
all_points.append({
't': t_current,
'mx': pos_missile[0],
'my': pos_missile[1],
'mz': pos_missile[2],
'cx': cloud_center[0],
'cy': cloud_center[1],
'cz': cloud_center[2],
'shield': shielded
})
if count % 1000 == 0:
progress = (t_current - t_start) / (t_end - t_start) * 100
print(f"进度: {progress:.1f}% (t={t_current:.3f}s, 累计遮蔽: {total_shield:.6f}s)")
count += 1
t_current += dt
if current_start is not None:
intervals.append((current_start, t_end))
print(f"\n有效遮蔽区间详情:")
if intervals:
total_interval = 0
for idx, (s, e) in enumerate(intervals):
dur = e - s
total_interval += dur
print(f"区间{idx + 1}: [{s:.6f}s, {e:.6f}s], 时长: {dur:.6f}s")
print(f"区间总时长验证: {total_interval:.8f}s")
else:
print("未发现有效遮蔽区间")
print(f"{'时间(s)':<10} {'导弹位置':<30} {'云团位置':<30} {'状态'}")
print("-" * 80)
for i in range(0, len(all_points), 100):
d = all_points[i]
m_str = f"({d['mx']:.1f},{d['my']:.1f},{d['mz']:.1f})"
c_str = f"({d['cx']:.1f},{d['cy']:.1f},{d['cz']:.1f})"
status = '遮蔽' if d['shield'] else '无遮蔽'
print(f"{d['t']:<10.3f} {m_str:<30} {c_str:<30} {status}")
print(f"有效遮蔽总时长: {total_shield:.10f} 秒")
print(f"遮蔽区间数量: {len(intervals)}")
print(f"计算精度: {dt} 秒/步")
print(f"总计算步数: {len(all_points)}")
print(f"数据完整性: 100%")
return total_shield, all_points
if __name__ == "__main__":
result, data = calculate_problem1_maximum_precision()
第二问
针对问题二,问题二是在问题一的基础上由原来四项固定的数据变量变为了未知量，我们需要自己来优化单机单弹的投放策略，我们需要同时优化四个参数即无人机飞行速度、无人机飞行角度、无人机投放时间、烟幕弹起爆延迟、四个参数相互耦合来寻找最长遮蔽时间的最优组合，因为搜索的数量非常庞大，我们采用网格粗略搜索+细致搜索+局部优化搜索的方式，寻找最优参数，寻找最长遮蔽时间的计算过程同问题一
import numpy as np
import math
import time
from numba import jit, njit
import warnings
from tqdm import tqdm, trange
warnings.filterwarnings('ignore')
@njit
def check_intersect(p1, p2, center, r):
d = p2 - p1
v = p1 - center
a = np.dot(d, d)
if a < 1e-12:
return np.linalg.norm(p1 - center) <= r
b = 2 * np.dot(d, v)
c = np.dot(v, v) - r * r
disc = b * b - 4 * a * c
if disc < 0:
return False
sqrt_disc = math.sqrt(disc)
t1 = (-b + sqrt_disc) / (2 * a)
t2 = (-b - sqrt_disc) / (2 * a)
return (0 <= t1 <= 1) or (0 <= t2 <= 1) or (t1 < 0 and t2 > 1)
@njit
def calc_shield(v_uav, angle, t_drop, t_delay, p_fy, p_m, u_dir, p_tgt, points, t_flight):
if not (70 <= v_uav <= 140 and 0 <= t_drop <= 50 and 0.01 <= t_delay <= 25):
return 0.0
dir_vec = np.array([math.cos(angle), math.sin(angle), 0.0])
p_drop = p_fy + v_uav * t_drop * dir_vec
dx = v_uav * t_delay * dir_vec
dz = np.array([0.0, 0.0, -0.5 * 9.8 * t_delay * t_delay])
p_burst = p_drop + dx + dz
t_burst = t_drop + t_delay
if t_burst >= t_flight or t_burst < 0:
return 0.0
dt = 0.005
shield_total = 0.0
t_start = t_burst
t_end = min(t_burst + 20.0, t_flight)
if t_start >= t_end:
return 0.0
t_cur = t_start
while t_cur <= t_end:
p_missile = p_m + 300.0 * t_cur * u_dir
sink_dist = 3.0 * (t_cur - t_burst)
p_cloud = p_burst + np.array([0.0, 0.0, -sink_dist])
shielded = False
for i in range(0, len(points), 3):
pt = points[i]
if check_intersect(p_missile, pt, p_cloud, 10.0):
shielded = True
break
if shielded:
shield_total += dt
t_cur += dt
return shield_total
class FastSmokeOptimizer:
def __init__(self):
self.p_fy = np.array([17800.0, 0.0, 1800.0])
self.p_m = np.array([20000.0, 0.0, 2000.0])
self.p_fake = np.array([0.0, 0.0, 0.0])
self.p_tgt = np.array([0.0, 200.0, 0.0])
dir_vec = self.p_fake - self.p_m
self.u_dir = dir_vec / np.linalg.norm(dir_vec)
self.t_flight = np.linalg.norm(dir_vec) / 300.0
self.points = self._gen_points()
print(f"导弹飞行时间: {self.t_flight:.3f}s")
print(f"采样点数: {len(self.points)}")
def _gen_points(self):
r = 7.0
h = 10.0
pts = []
pts.append(self.p_tgt.copy())
angles = np.linspace(0, 2 * np.pi, 36, endpoint=False)
heights = np.linspace(-h / 2, h / 2, 11)
for z in heights:
for a in angles:
x = self.p_tgt[0] + r * np.cos(a)
y = self.p_tgt[1] + r * np.sin(a)
z_val = self.p_tgt[2] + z
pts.append(np.array([x, y, z_val]))
return np.array(pts)
def calc_shield_time(self, params):
return calc_shield(
params[0], params[1], params[2], params[3],
self.p_fy, self.p_m, self.u_dir,
self.p_tgt, self.points, self.t_flight
)
def smart_search(self):
best_params = None
best_score = 0
speeds = np.linspace(70, 140, 8)
angles = np.linspace(0, 2 * np.pi, 12)
drop_times = np.linspace(0.5, 35, 12)
delays = np.linspace(0.1, 20, 12)
total_pts = len(speeds) * len(angles) * len(drop_times) * len(delays)
good_regions = []
with tqdm(total=total_pts, desc="第一层搜索") as pbar:
for v in speeds:
for a in angles:
for t_d in drop_times:
for t_del in delays:
p = [v, a, t_d, t_del]
s = self.calc_shield_time(p)
if s > best_score:
best_score = s
best_params = p.copy()
pbar.set_postfix(最佳=f"{best_score:.3f}s")
if s > 2.0:
good_regions.append({
'p': p.copy(),
's': s,
'r': [v, a, t_d, t_del]
})
pbar.update(1)
print(f"第一层完成，最佳: {best_score:.6f}s，良好区域: {len(good_regions)}")
if not good_regions:
return best_params, best_score
good_regions.sort(key=lambda x: x['s'], reverse=True)
top_regions = good_regions[:max(1, len(good_regions) // 5)]
total_pts2 = len(top_regions) * 8 * 8 * 8 * 8
with tqdm(total=total_pts2, desc="第二层搜索") as pbar:
for reg in top_regions:
c = reg['r']
v_range = np.linspace(max(70, c[0] - 15), min(140, c[0] + 15), 8)
a_range = np.linspace(c[1] - 0.8, c[1] + 0.8, 8)
d_range = np.linspace(max(0.1, c[2] - 5), c[2] + 5, 8)
del_range = np.linspace(max(0.1, c[3] - 3), c[3] + 3, 8)
for v in v_range:
for a in a_range:
for t_d in d_range:
for t_del in del_range:
p = [v, a, t_d, t_del]
s = self.calc_shield_time(p)
if s > best_score:
best_score = s
best_params = p.copy()
pbar.set_postfix(最佳=f"{best_score:.3f}s")
pbar.update(1)
return best_params, best_score
def refine_search(self, start_p, steps):
print("开始局部优化...")
cur_p = np.array(start_p)
cur_s = self.calc_shield_time(cur_p)
dirs = np.array([
[1, 0, 0, 0], [-1, 0, 0, 0],
[0, 1, 0, 0], [0, -1, 0, 0],
[0, 0, 1, 0], [0, 0, -1, 0],
[0, 0, 0, 1], [0, 0, 0, -1]
])
improved = True
iter_count = 0
max_iter = 500
with tqdm(total=max_iter, desc="局部优化") as pbar:
while improved and iter_count < max_iter:
improved = False
iter_count += 1
for d in dirs:
new_p = cur_p + d * steps
new_p[0] = max(70, min(140, new_p[0]))
new_p[1] = new_p[1] % (2 * np.pi)
new_p[2] = max(0.1, min(50, new_p[2]))
new_p[3] = max(0.01, min(25, new_p[3]))
new_s = self.calc_shield_time(new_p)
if new_s > cur_s:
cur_p = new_p
cur_s = new_s
improved = True
pbar.set_postfix(最佳=f"{cur_s:.3f}s")
break
if not improved:
steps *= 0.8
if np.max(steps) < 1e-3:
break
pbar.update(1)
print(f"局部优化完成，迭代{iter_count}次，结果: {cur_s:.6f}s")
return cur_p, cur_s
def analyze_result(self, params, score):
v, a, t_drop, t_del = params
print(f"遮蔽时间: {score:.6f}s")
print(f"速度: {v:.2f} m/s")
print(f"角度: {math.degrees(a):.2f}°")
print(f"投放时间: {t_drop:.3f}s")
print(f"延迟: {t_del:.3f}s")
dir_vec = np.array([math.cos(a), math.sin(a), 0])
p_drop = self.p_fy + v * t_drop * dir_vec
p_burst = p_drop + v * t_del * dir_vec + np.array([0, 0, -0.5 * 9.8 * t_del ** 2])
print(f"起爆点: ({p_burst[0]:.0f}, {p_burst[1]:.0f}, {p_burst[2]:.0f})")
print(f"起爆时间: {t_drop + t_del:.3f}s")
def main():
start_t = time.time()
opt = FastSmokeOptimizer()
grid_p, grid_s = opt.smart_search()
t1 = time.time() - start_t
print(f"第一阶段用时: {t1:.1f}s，结果: {grid_s:.6f}s")
if grid_p is None:
print("无有效解")
return
steps = np.array([5.0, 0.2, 2.0, 1.0])
final_p, final_s = opt.refine_search(grid_p, steps)
total_t = time.time() - start_t
opt.analyze_result(final_p, final_s)
print(f"\n性能:")
print(f"总时间: {total_t:.1f}s")
print(f"最终结果: {final_s:.6f}s")
base_p = [120.0, math.pi, 1.5, 3.6]
base_s = opt.calc_shield_time(base_p)
ratio = final_s / base_s if base_s > 0 else 0
return opt, final_p, final_s
if __name__ == "__main__":
main()问题三
针对问题三，问题三是在问题二与问题一基础上由优化单机单弹策略扩展为优化单机多弹，问题三需要同时在时间与空间上优化八个参数，且有多个约束条件，对此我们抽象为一个约束优化来寻找最佳策略的最长遮蔽时间，对此我们在最优控制理论问题的研究上，选择了遗传算法作为解决本问题的主要算法核心原理，将时间与空间参数视作遗传算法的基因变量，结合最优控制问题的基本思想，将烟幕弹轨迹视为控制变量，用适应度函数考虑遮蔽时间和时间连续性，用遗传算法求解高维非线性优化。
import numpy as np
import math
import sys
# -----------------------
# 可选依赖：pandas（导出Excel），不可用则回退 openpyxl
# -----------------------
HAS_PANDAS = True
try:
import pandas as pd
except Exception as e:
print("Pandas 不可用，自动回退到 openpyxl 写 Excel：", e)
HAS_PANDAS = False
try:
from openpyxl import Workbook
HAS_OPENPYXL = True
except Exception as e2:
HAS_OPENPYXL = False
print("openpyxl 也不可用，无法写 Excel：", e2)
# -----------------------
# 可选依赖：matplotlib（可视化）
# -----------------------
HAS_MPL = True
try:
import matplotlib.pyplot as plt
except Exception as e:
HAS_MPL = False
print("Matplotlib 不可用，跳过绘图：", e)
# -----------------------
# 可选依赖：Numba（加速几何判定）
# -----------------------
try:
from numba import njit
HAS_NUMBA = True
except Exception as e:
HAS_NUMBA = False
print("Numba 不可用，将使用纯Python判定（会慢一些）：", e)
def njit(func=None, **kwargs):
if func is None:
def wrapper(f):
return f
return wrapper
return func
# -----------------------
# 常量与几何配置
# -----------------------
g = 9.8
MISSILE_SPEED = 300.0
SMOKE_RADIUS = 10.0
SMOKE_SINK = 3.0
SMOKE_LIFE = 20.0
V_MIN, V_MAX = 70.0, 140.0
FY1_INIT = np.array([17800.0, 0.0, 1800.0])
M1_INIT = np.array([20000.0, 0.0, 2000.0])
FAKE_TGT = np.array([0.0, 0.0, 0.0])
REAL_TGT_CENTER = np.array([0.0, 200.0, 0.0])
missile_dir = FAKE_TGT - M1_INIT
missile_unit_dir = missile_dir / np.linalg.norm(missile_dir)
T_FLIGHT = np.linalg.norm(missile_dir) / MISSILE_SPEED # 导弹到达假目标时间
# 目标圆柱采样（72×21 + 中心）与低密度版本
def generate_cylinder_points(high_density=True):
points = []
cylinder_radius = 7.0
cylinder_height = 10.0
points.append(REAL_TGT_CENTER.copy())
if high_density:
angles = np.linspace(0, 2*np.pi, 72, endpoint=False)
heights = np.linspace(-cylinder_height/2, cylinder_height/2, 21)
else:
angles = np.linspace(0, 2*np.pi, 36, endpoint=False)
heights = np.linspace(-cylinder_height/2, cylinder_height/2, 11)
for h in heights:
for a in angles:
x = REAL_TGT_CENTER[0] + cylinder_radius*np.cos(a)
y = REAL_TGT_CENTER[1] + cylinder_radius*np.sin(a)
z = REAL_TGT_CENTER[2] + h
points.append(np.array([x,y,z]))
return np.array(points, dtype=float)
CYL_POINTS_FAST = generate_cylinder_points(high_density=False) # GA评估
CYL_POINTS_FULL = generate_cylinder_points(high_density=True) # 复核
# 线段-球体相交（Numba可选）
@njit
def line_intersects_sphere_numba(line_start, line_end, sphere_center, sphere_radius):
d = line_end - line_start
f = line_start - sphere_center
a = d[0]*d[0] + d[1]*d[1] + d[2]*d[2]
if a < 1e-12:
dx = line_start[0] - sphere_center[0]
dy = line_start[1] - sphere_center[1]
dz = line_start[2] - sphere_center[2]
return (dx*dx + dy*dy + dz*dz) <= sphere_radius * sphere_radius
b = 2.0 * (d[0]*f[0] + d[1]*f[1] + d[2]*f[2])
c = (f[0]*f[0] + f[1]*f[1] + f[2]*f[2]) - sphere_radius * sphere_radius
disc = b*b - 4*a*c
if disc < 0.0:
return False
sqrt_disc = math.sqrt(disc)
t1 = (-b - sqrt_disc) / (2*a)
t2 = (-b + sqrt_disc) / (2*a)
return (0.0 <= t1 <= 1.0) or (0.0 <= t2 <= 1.0) or (t1 < 0.0 and t2 > 1.0)
# 策略解码（排序+间隔约束）
def decode_strategy(chrom):
# chrom: [v, theta, t1_raw, t2_raw, t3_raw, d1_raw, d2_raw, d3_raw]
v = float(np.clip(chrom[0], V_MIN, V_MAX))
theta = float(chrom[1] % (2*np.pi))
ts_raw = np.clip(chrom[2:5], 0.0, T_FLIGHT)
ts = np.sort(ts_raw)
t1 = float(ts[0])
t2 = float(max(ts[1], t1 + 1.0))
t3 = float(max(ts[2], t2 + 1.0))
t3 = float(min(t3, T_FLIGHT))
d1 = float(np.clip(chrom[5], 0.01, SMOKE_LIFE))
d2 = float(np.clip(chrom[6], 0.01, SMOKE_LIFE))
d3 = float(np.clip(chrom[7], 0.01, SMOKE_LIFE))
return v, theta, np.array([t1, t2, t3], dtype=float), np.array([d1, d2, d3], dtype=float)
def burst_from(v, theta, t, d):
dir2 = np.array([math.cos(theta), math.sin(theta), 0.0], dtype=float)
p_drop = FY1_INIT + v * t * dir2
tau = t + d
p_burst = p_drop + v * d * dir2 + np.array([0.0, 0.0, -0.5*g*d*d], dtype=float)
return tau, p_burst
# 遮蔽并集（快速评估）
def compute_occlusion_union(v, theta, ts, ds, dt=0.005, cyl_points=None, stride=3):
if cyl_points is None:
cyl_points = CYL_POINTS_FAST
bursts = []
for i in range(3):
tau, p_burst = burst_from(v, theta, ts[i], ds[i])
if tau > T_FLIGHT:
bursts.append((np.inf, p_burst))
else:
bursts.append((tau, p_burst))
t0 = min([b[0] for b in bursts] + [T_FLIGHT])
if not np.isfinite(t0):
return 0.0, [], {}
t_start = max(0.0, t0)
t_end = T_FLIGHT
if t_start >= t_end:
return 0.0, [], {}
times = np.arange(t_start, t_end + 1e-12, dt, dtype=float)
occl_mask = np.zeros(times.shape[0], dtype=bool)
per_bomb_mask = [np.zeros_like(occl_mask) for _ in range(3)]
for idx, t in enumerate(times):
pm = M1_INIT + MISSILE_SPEED * t * missile_unit_dir
any_hit = False
for j in range(3):
tau, pb = bursts[j]
if not np.isfinite(tau):
continue
if t < tau or t > min(tau + SMOKE_LIFE, T_FLIGHT):
continue
c = pb + np.array([0.0, 0.0, -SMOKE_SINK*(t - tau)], dtype=float)
hit = False
for k in range(0, cyl_points.shape[0], stride):
if line_intersects_sphere_numba(pm, cyl_points[k], c, SMOKE_RADIUS):
hit = True
break
if hit:
any_hit = True
per_bomb_mask[j][idx] = True
occl_mask[idx] = any_hit
total = float(occl_mask.sum() * dt)
intervals = []
in_seg = False
seg_start = 0.0
for i in range(len(times)):
if occl_mask[i] and not in_seg:
in_seg = True
seg_start = times[i]
if (not occl_mask[i] and in_seg) or (i == len(times)-1 and in_seg):
seg_end = times[i] if not occl_mask[i] else times[i]
intervals.append((round(seg_start,5), round(seg_end,5)))
in_seg = False
details = dict(times=times, mask=occl_mask, per_bomb=per_bomb_mask, bursts=bursts)
return total, intervals, details
# 遮蔽并集（高精度复核）
def compute_occlusion_union_full(v, theta, ts, ds, dt=0.001):
return compute_occlusion_union(v, theta, ts, ds, dt=dt, cyl_points=CYL_POINTS_FULL, stride=1)
# -----------------------
# 遗传算法
# -----------------------
def init_population(pop_size=120):
pop = []
for _ in range(pop_size):
v = np.random.uniform(V_MIN, V_MAX)
theta = np.random.uniform(0, 2*np.pi)
t1 = np.random.uniform(0.0, max(0.01, T_FLIGHT - 2.0))
t2 = t1 + np.random.uniform(1.0, 3.0)
t3 = t2 + np.random.uniform(1.0, 3.0)
d1 = np.random.uniform(0.2, 10.0)
d2 = np.random.uniform(0.2, 10.0)
d3 = np.random.uniform(0.2, 10.0)
pop.append(np.array([v, theta, t1, t2, t3, d1, d2, d3], dtype=float))
return np.array(pop, dtype=float)
def fitness(chrom, fast_dt=0.005):
v, th, ts, ds = decode_strategy(chrom)
taus = ts + ds
if np.min(taus) > T_FLIGHT:
return -1e6
total, _, _ = compute_occlusion_union(v, th, ts, ds, dt=fast_dt)
reg = -0.001 * np.sum(taus) # 轻微鼓励早起爆
return total + reg
def select_parents(pop, fits, k=2):
idx = np.random.choice(len(pop), size=min(6, len(pop)), replace=False)
best_two = idx[np.argsort(fits[idx])][-k:]
return pop[best_two[0]], pop[best_two[1]]
def crossover(p1, p2, pc=0.9):
if np.random.rand() > pc:
return p1.copy(), p2.copy()
alpha = np.random.uniform(-0.2, 1.2, size=p1.shape)
c1 = alpha * p1 + (1-alpha) * p2
c2 = alpha * p2 + (1-alpha) * p1
return c1, c2
def mutate(c, pm=0.2, scales=None):
if scales is None:
scales = np.array([5.0, 0.2, 1.0, 1.0, 1.0, 0.6, 0.6, 0.6], dtype=float)
for i in range(len(c)):
if np.random.rand() < pm:
c[i] += np.random.randn() * float(scales[i])
return c
def run_ga(gens=60, pop_size=120, elite=6, fast_dt=0.005, verbose=True):
pop = init_population(pop_size)
fits = np.array([fitness(ind, fast_dt=fast_dt) for ind in pop], dtype=float)
best_hist = []
for g in range(gens):
new_pop = []
elite_idx = np.argsort(fits)[-elite:]
for ei in elite_idx:
new_pop.append(pop[ei].copy())
while len(new_pop) < pop_size:
p1, p2 = select_parents(pop, fits, k=2)
c1, c2 = crossover(p1, p2, pc=0.9)
c1 = mutate(c1, pm=0.2)
c2 = mutate(c2, pm=0.2)
new_pop.append(c1)
if len(new_pop) < pop_size:
new_pop.append(c2)
pop = np.array(new_pop, dtype=float)
fits = np.array([fitness(ind, fast_dt=fast_dt) for ind in pop], dtype=float)
b_idx = int(np.argmax(fits))
best_hist.append(float(fits[b_idx]))
if verbose and (g % 5 == 0 or g == gens-1):
print(f"Gen {g+1}/{gens} best≈{fits[b_idx]:.4f} s")
best_idx = int(np.argmax(fits))
return pop[best_idx], float(fits[best_idx]), best_hist
# -----------------------
# 主流程：优化、复核、导出、可视化
# -----------------------
def main_problem3():
print("问题三：FY1 对 M1 投放3枚烟幕弹的最优控制 + 遗传算法求解")
print(f"导弹总飞行时间 T_flight={T_FLIGHT:.3f} s")
# 1) GA 快速搜索
best_chrom, best_fit_fast, hist = run_ga(gens=60, pop_size=120, elite=8, fast_dt=0.005, verbose=True)
v, theta, ts, ds = decode_strategy(best_chrom)
print(f"\nGA-FAST 最佳：v={v:.2f} m/s, θ={math.degrees(theta):.2f}°, t={ts}, δ={ds}, fit≈{best_fit_fast:.4f} s")
# 2) 高精度复核
total_full, intervals_full, details = compute_occlusion_union_full(v, theta, ts, ds, dt=0.001)
print(f"高精度复核并集遮蔽时长：{total_full:.5f} s")
print(f"遮蔽区间（前若干）：{intervals_full[:5]} ... 共 {len(intervals_full)} 段")
# 3) 计算每弹起爆信息
rows = []
for i in range(3):
tau, p_burst = burst_from(v, theta, ts[i], ds[i])
rows.append({
'UAV_ID':'FY1','Missile_ID':'M1',
'speed_mps': round(v,5),
'heading_deg': round(math.degrees(theta),5),
'drop_time_s': round(float(ts[i]),5),
'burst_delay_s': round(float(ds[i]),5),
'burst_time_s': round(float(tau),5),
'burst_x': round(float(p_burst[0]),5),
'burst_y': round(float(p_burst[1]),5),
'burst_z': round(float(p_burst[2]),5)
})
# 4) 导出到 result1.xlsx
try:
if HAS_PANDAS:
df = pd.DataFrame(rows)
df_attr = pd.DataFrame([{
'total_union_occlusion_s': round(float(total_full),5),
'T_flight_s': round(float(T_FLIGHT),5),
'dt_eval_s': 0.001
}])
with pd.ExcelWriter('result1.xlsx') as writer:
df.to_excel(writer, sheet_name='strategy', index=False)
df_attr.to_excel(writer, sheet_name='summary', index=False)
elif 'HAS_OPENPYXL' in globals() and HAS_OPENPYXL:
wb = Workbook()
ws1 = wb.active
ws1.title = 'strategy'
headers = list(rows[0].keys())
ws1.append(headers)
for r in rows:
ws1.append([r[h] for h in headers])
ws2 = wb.create_sheet('summary')
ws2.append(['total_union_occlusion_s','T_flight_s','dt_eval_s'])
ws2.append([round(float(total_full),5), round(float(T_FLIGHT),5), 0.001])
wb.save('result1.xlsx')
else:
print("未能导出 Excel：既无 pandas 也无 openpyxl。")
except Exception as e:
print("导出 result1.xlsx 失败：", e)
print("导出完成（若无错误信息则已生成 result1.xlsx）。")
# 5) 可视化：遮蔽时间轴
if HAS_MPL:
try:
times = details['times']
mask = details['mask']
plt.figure(figsize=(10,2.4))
plt.plot(times, mask.astype(int), drawstyle='steps-mid')
plt.yticks([0,1], ['No','Occluded'])
plt.xlabel('Time (s)')
plt.title('Occlusion Timeline (Union of 3 Clouds)')
plt.grid(True, axis='x', alpha=0.3)
plt.tight_layout()
plt.savefig('occlusion_timeline.png', dpi=150)
print("已保存 occlusion_timeline.png")
except Exception as e:
print("保存遮蔽时间轴失败：", e)
# 6) 3D几何示意
try:
from mpl_toolkits.mplot3d import Axes3D # noqa: F401
fig = plt.figure(figsize=(6,5))
ax = fig.add_subplot(111, projection='3d')
t_show = np.linspace(0, T_FLIGHT, 100)
pm = M1_INIT + np.outer(t_show*MISSILE_SPEED, missile_unit_dir)
ax.plot(pm[:,0], pm[:,1], pm[:,2], 'r-', label='Missile Path')
for i in range(3):
tau, p_burst = burst_from(v, theta, ts[i], ds[i])
ax.scatter([p_burst[0]],[p_burst[1]],[p_burst[2]], s=40, label=f'Burst{i+1}')
ax.scatter([REAL_TGT_CENTER[0]],[REAL_TGT_CENTER[1]],[REAL_TGT_CENTER[2]], c='g', s=40, label='Real Target Center')
ax.set_xlabel('X'); ax.set_ylabel('Y'); ax.set_zlabel('Z')
ax.set_title('Geometry Overview')
ax.legend()
plt.tight_layout()
plt.savefig('geometry_overview.png', dpi=150)
print("已保存 geometry_overview.png")
except Exception as e:
print("3D绘图失败（可忽略）：", e)
return dict(v=v, theta=theta, ts=ts, ds=ds, total=total_full, intervals=intervals_full)
if __name__ == '__main__':
np.set_printoptions(precision=5, suppress=True)
res = main_problem3()
# 可在此打印关键参数
if res:
print("\n最优参数（复核）：")
print(f"speed={res['v']:.3f} m/s, heading={math.degrees(res['theta']):.3f} deg")
print(f"drop_times(s)={res['ts']}, delays(s)={res['ds']}")
print(f"union_occlusion={res['total']:.5f} s")
问题四
针对问题四，问题四是在前三个问题的基础上由优化单机多弹策略转变为优化多机单弹，是承上启下的，因为FY1，FY2，FY3在x,y,z上的坐标是逐层递减的，可粗略视作在一条直线上，而导弹的飞行速度极快，所以这里我们可以分区间来讨论，依据导弹飞行速度与方向，我们将三个无人机的行进方向都定于朝向原点方向，并以朝向原点方向为基准适当的扩展几度的角度，这样可以大幅减少角度方向的计算量，遗传算法应用时对于无人机航向角就不用随机产生了而是只产生无人机朝向原点角度±2度的，而无人机的速度与掷蛋和爆弹时间则应与导弹此刻的位置产生一些关联，不能随意随机，我的策略是等前一个烟雾已经无法产生有效遮蔽时间时且导弹快追上下一个无人机时，无人机掷弹并优化爆弹时间来最大化有效遮蔽时间，我们对此设置一个合理的区间可以大大减少运算量。剩余的计算问题则继续采用问题二与问题三的思路与方式。
import numpy as np
import math
import time
from numba import njit
import random
from typing import List, Tuple, Dict
# 系统参数
FY1_INIT = np.array([17800.0, 0.0, 1800.0])
FY2_INIT = np.array([9000.0, 0.0, 900.0])
FY3_INIT = np.array([4500.0, 0.0, 450.0])
M1_INIT = np.array([20000.0, 0.0, 2000.0])
M2_INIT = np.array([15000.0, 5000.0, 1500.0])
ORIGIN = np.array([0.0, 0.0, 0.0])
REAL_TARGET = np.array([0.0, 200.0, 0.0])
MISSILE_SPEED = 300.0
SMOKE_RADIUS = 10.0
SMOKE_DURATION = 20.0
SINK_SPEED = 3.0
g = 9.8
# 问题二的最优FY1参数（固定）
FY1_OPTIMAL = {
'speed': 140.0, # m/s
'angle': math.radians(6.55), # 转换为弧度
'drop_time': 0.1, # s
'burst_delay': 0.529 # s
}
class ImprovedMultiUAVOptimizer:
def __init__(self):
self.uav_positions = {
'FY1': FY1_INIT,
'FY2': FY2_INIT,
'FY3': FY3_INIT
}
# 导弹参数
self.missiles = {
'M1': {'pos': M1_INIT, 'target': ORIGIN},
'M2': {'pos': M2_INIT, 'target': ORIGIN}
}
# 计算导弹飞行参数
self.missile_params = {}
for m_id, m_info in self.missiles.items():
direction = m_info['target'] - m_info['pos']
flight_time = np.linalg.norm(direction) / MISSILE_SPEED
unit_dir = direction / np.linalg.norm(direction)
self.missile_params[m_id] = {
'flight_time': flight_time,
'unit_direction': unit_dir,
'initial_pos': m_info['pos']
}
# 目标采样点
self.target_points = self.generate_target_points()
# 计算FY1的遮蔽时间窗口
self.fy1_coverage = self.calculate_fy1_coverage()
print("系统初始化完成:")
print(f"M1飞行时间: {self.missile_params['M1']['flight_time']:.2f}s")
print(f"M2飞行时间: {self.missile_params['M2']['flight_time']:.2f}s")
print(f"FY1遮蔽窗口: {self.fy1_coverage}")
def generate_target_points(self):
"""生成目标区域采样点"""
points = [REAL_TARGET.copy()]
radius = 7.0
height = 10.0
angles = np.linspace(0, 2 * np.pi, 20, endpoint=False)
heights = np.linspace(-height / 2, height / 2, 7)
for h in heights:
for a in angles:
x = REAL_TARGET[0] + radius * np.cos(a)
y = REAL_TARGET[1] + radius * np.sin(a)
z = REAL_TARGET[2] + h
points.append(np.array([x, y, z]))
return np.array(points)
def calculate_fy1_coverage(self):
"""计算FY1的遮蔽时间窗口"""
# 根据FY1参数计算起爆位置和时间
speed = FY1_OPTIMAL['speed']
angle = FY1_OPTIMAL['angle']
drop_time = FY1_OPTIMAL['drop_time']
burst_delay = FY1_OPTIMAL['burst_delay']
flight_dir = np.array([math.cos(angle), math.sin(angle), 0.0])
drop_pos = FY1_INIT + speed * drop_time * flight_dir
burst_pos = drop_pos + speed * burst_delay * flight_dir + \
np.array([0.0, 0.0, -0.5 * g * burst_delay ** 2])
burst_time = drop_time + burst_delay
end_time = burst_time + SMOKE_DURATION
return {
'start_time': burst_time,
'end_time': end_time,
'burst_pos': burst_pos
}
def calculate_strategic_windows(self):
"""计算FY2和FY3的战略投放窗口"""
fy1_end = self.fy1_coverage['end_time']
windows = {}
# FY2窗口：FY1结束前2-5秒开始，确保连续覆盖
windows['FY2'] = {
'drop_start': max(0.1, fy1_end - 8.0),
'drop_end': fy1_end - 1.0,
'optimal_burst_time': fy1_end - 2.0 # 理想起爆时间
}
# FY3窗口：FY2预期结束前开始
fy2_expected_end = windows['FY2']['optimal_burst_time'] + SMOKE_DURATION
windows['FY3'] = {
'drop_start': max(0.1, fy2_expected_end - 10.0),
'drop_end': fy2_expected_end - 1.0,
'optimal_burst_time': fy2_expected_end - 2.0
}
return windows
@njit
def line_intersects_sphere_fast(line_start, line_end, sphere_center, sphere_radius):
"""快速线-球相交检测"""
d = line_end - line_start
f = line_start - sphere_center
a = np.dot(d, d)
if a < 1e-12:
return np.linalg.norm(line_start - sphere_center) <= sphere_radius
b = 2.0 * np.dot(d, f)
c = np.dot(f, f) - sphere_radius * sphere_radius
discriminant = b * b - 4 * a * c
if discriminant < 0:
return False
sqrt_disc = math.sqrt(discriminant)
t1 = (-b - sqrt_disc) / (2 * a)
t2 = (-b + sqrt_disc) / (2 * a)
return (0.0 <= t1 <= 1.0) or (0.0 <= t2 <= 1.0) or (t1 < 0.0 and t2 > 1.0)
class ConstrainedGeneticAlgorithm:
"""约束遗传算法 - 只优化FY2和FY3"""
def __init__(self, optimizer):
self.optimizer = optimizer
self.pop_size = 60
self.generations = 100
self.mutation_rate = 0.2
# 计算战略窗口
self.windows = optimizer.calculate_strategic_windows()
# 基因组结构: [FY2_speed, FY2_angle, FY2_drop_time, FY2_burst_delay,
# FY3_speed, FY3_angle, FY3_drop_time, FY3_burst_delay]
self.genome_length = 8
print("优化参数空间:")
for uav_id, window in self.windows.items():
print(f"{uav_id}: 投放窗口 [{window['drop_start']:.1f}s, {window['drop_end']:.1f}s]")
def initialize_population(self):
"""智能初始化种群"""
population = []
for _ in range(self.pop_size):
genome = []
# FY2参数
fy2_window = self.windows['FY2']
fy2_speed = random.uniform(80, 140)
# FY2角度：朝向原点方向±3度
fy2_to_origin = ORIGIN - FY2_INIT
fy2_base_angle = math.atan2(fy2_to_origin[1], fy2_to_origin[0])
fy2_angle = random.uniform(fy2_base_angle - math.radians(3),
fy2_base_angle + math.radians(3))
# FY2时间参数
fy2_drop_time = random.uniform(fy2_window['drop_start'], fy2_window['drop_end'])
# 根据目标起爆时间反推延迟
target_burst = fy2_window['optimal_burst_time']
ideal_delay = target_burst - fy2_drop_time
fy2_burst_delay = max(0.1, min(15.0, random.gauss(ideal_delay, 2.0)))
genome.extend([fy2_speed, fy2_angle, fy2_drop_time, fy2_burst_delay])
# FY3参数
fy3_window = self.windows['FY3']
fy3_speed = random.uniform(80, 140)
# FY3角度：朝向原点方向±3度
fy3_to_origin = ORIGIN - FY3_INIT
fy3_base_angle = math.atan2(fy3_to_origin[1], fy3_to_origin[0])
fy3_angle = random.uniform(fy3_base_angle - math.radians(3),
fy3_base_angle + math.radians(3))
# FY3时间参数
fy3_drop_time = random.uniform(fy3_window['drop_start'], fy3_window['drop_end'])
target_burst = fy3_window['optimal_burst_time']
ideal_delay = target_burst - fy3_drop_time
fy3_burst_delay = max(0.1, min(15.0, random.gauss(ideal_delay, 2.0)))
genome.extend([fy3_speed, fy3_angle, fy3_drop_time, fy3_burst_delay])
population.append(np.array(genome))
return population
def calculate_total_shielding(self, genome):
"""计算总遮蔽时间（包含FY1固定贡献）"""
try:
# 解码FY2和FY3参数
fy2_params = {
'speed': float(genome[0]),
'angle': float(genome[1]),
'drop_time': float(genome[2]),
'burst_delay': float(genome[3])
}
fy3_params = {
'speed': float(genome[4]),
'angle': float(genome[5]),
'drop_time': float(genome[6]),
'burst_delay': float(genome[7])
}
# 计算所有烟云事件
smoke_events = []
# FY1固定事件
fy1_coverage = self.optimizer.fy1_coverage
smoke_events.append({
'start_time': fy1_coverage['start_time'],
'end_time': fy1_coverage['end_time'],
'burst_pos': fy1_coverage['burst_pos']
})
# FY2事件
fy2_dir = np.array([math.cos(fy2_params['angle']), math.sin(fy2_params['angle']), 0.0])
fy2_drop_pos = FY2_INIT + fy2_params['speed'] * fy2_params['drop_time'] * fy2_dir
fy2_burst_pos = fy2_drop_pos + fy2_params['speed'] * fy2_params['burst_delay'] * fy2_dir + \
np.array([0.0, 0.0, -0.5 * g * fy2_params['burst_delay'] ** 2])
fy2_burst_time = fy2_params['drop_time'] + fy2_params['burst_delay']
if fy2_burst_pos[2] >= 0: # 地面以上
smoke_events.append({
'start_time': fy2_burst_time,
'end_time': fy2_burst_time + SMOKE_DURATION,
'burst_pos': fy2_burst_pos
})
# FY3事件
fy3_dir = np.array([math.cos(fy3_params['angle']), math.sin(fy3_params['angle']), 0.0])
fy3_drop_pos = FY3_INIT + fy3_params['speed'] * fy3_params['drop_time'] * fy3_dir
fy3_burst_pos = fy3_drop_pos + fy3_params['speed'] * fy3_params['burst_delay'] * fy3_dir + \
np.array([0.0, 0.0, -0.5 * g * fy3_params['burst_delay'] ** 2])
fy3_burst_time = fy3_params['drop_time'] + fy3_params['burst_delay']
if fy3_burst_pos[2] >= 0:
smoke_events.append({
'start_time': fy3_burst_time,
'end_time': fy3_burst_time + SMOKE_DURATION,
'burst_pos': fy3_burst_pos
})
if len(smoke_events) == 0:
return 0.0
# 计算总遮蔽时间
total_shielding = 0.0
dt = 0.01
min_time = min([event['start_time'] for event in smoke_events])
max_time = max([event['end_time'] for event in smoke_events])
# 对每个导弹计算
for m_id, m_params in self.optimizer.missile_params.items():
flight_time = m_params['flight_time']
unit_dir = m_params['unit_direction']
init_pos = m_params['initial_pos']
t = max(0.0, min_time)
while t <= min(max_time, flight_time):
missile_pos = init_pos + MISSILE_SPEED * t * unit_dir
# 检查是否被任一烟云遮蔽
is_shielded = False
for event in smoke_events:
if event['start_time'] <= t <= event['end_time']:
sink_distance = SINK_SPEED * (t - event['start_time'])
current_center = event['burst_pos'] + np.array([0.0, 0.0, -sink_distance])
# 检查遮蔽（采样点稀疏化以提高速度）
for i in range(0, len(self.optimizer.target_points), 2):
point = self.optimizer.target_points[i]
if line_intersects_sphere_fast(missile_pos, point, current_center, SMOKE_RADIUS):
is_shielded = True
break
if is_shielded:
break
if is_shielded:
total_shielding += dt
t += dt
return total_shielding
except Exception as e:
return 0.0
def fitness_function(self, individual):
"""适应度函数"""
# 基本遮蔽时间
shielding_time = self.calculate_total_shielding(individual)
# 约束惩罚
penalty = 0.0
# 参数边界约束
if not (70 <= individual[0] <= 140): # FY2速度
penalty += 2.0
if not (70 <= individual[4] <= 140): # FY3速度
penalty += 2.0
# 时间约束
fy2_window = self.windows['FY2']
if not (fy2_window['drop_start'] <= individual[2] <= fy2_window['drop_end']):
penalty += 3.0
fy3_window = self.windows['FY3']
if not (fy3_window['drop_start'] <= individual[6] <= fy3_window['drop_end']):
penalty += 3.0
# 时序协调奖励：鼓励连续覆盖
fy1_end = self.optimizer.fy1_coverage['end_time']
fy2_start = individual[2] + individual[3] # FY2起爆时间
fy3_start = individual[6] + individual[7] # FY3起爆时间
# 连续性奖励
if abs(fy2_start - fy1_end) < 2.0: # FY2与FY1连接好
shielding_time += 0.5
if abs(fy3_start - (fy2_start + SMOKE_DURATION)) < 2.0: # FY3与FY2连接好
shielding_time += 0.5
return max(0.0, shielding_time - penalty)
def tournament_selection(self, population, fitness_scores, k=3):
"""锦标赛选择"""
indices = random.choices(range(len(population)), k=k)
winner_idx = max(indices, key=lambda i: fitness_scores[i])
return population[winner_idx].copy()
def crossover(self, parent1, parent2):
"""交叉操作"""
child1 = parent1.copy()
child2 = parent2.copy()
# 按无人机分段交叉
if random.random() < 0.5:
# 交换FY2参数 (前4个基因)
child1[:4] = parent2[:4]
child2[:4] = parent1[:4]
if random.random() < 0.5:
# 交换FY3参数 (后4个基因)
child1[4:] = parent2[4:]
child2[4:] = parent1[4:]
return child1, child2
def mutate(self, individual):
"""变异操作"""
mutated = individual.copy()
for i in range(len(mutated)):
if random.random() < self.mutation_rate:
if i in [0, 4]: # 速度
mutated[i] += random.gauss(0, 10)
mutated[i] = np.clip(mutated[i], 70, 140)
elif i in [1, 5]: # 角度
mutated[i] += random.gauss(0, math.radians(5))
elif i in [2, 6]: # 投放时间
mutated[i] += random.gauss(0, 2.0)
else: # 起爆延迟
mutated[i] += random.gauss(0, 1.0)
mutated[i] = max(0.1, mutated[i])
return mutated
def optimize(self):
"""主优化循环"""
print("开始约束遗传算法优化...")
population = self.initialize_population()
best_individual = None
best_fitness = 0.0
history = []
for generation in range(self.generations):
# 计算适应度
fitness_scores = [self.fitness_function(ind) for ind in population]
# 更新最优解
current_best_idx = np.argmax(fitness_scores)
current_best = fitness_scores[current_best_idx]
if current_best > best_fitness:
best_fitness = current_best
best_individual = population[current_best_idx].copy()
history.append(best_fitness)
# 进度输出
if (generation + 1) % 15 == 0:
avg_fitness = np.mean(fitness_scores)
valid_count = sum(1 for f in fitness_scores if f > 0)
print(f"第{generation + 1}代: 最佳={best_fitness:.4f}s, "
f"平均={avg_fitness:.4f}s, 有效解={valid_count}/{self.pop_size}")
# 生成新一代
new_population = []
# 精英保留
elite_size = max(1, self.pop_size // 10)
elite_indices = np.argsort(fitness_scores)[-elite_size:]
for idx in elite_indices:
new_population.append(population[idx].copy())
# 生成后代
while len(new_population) < self.pop_size:
parent1 = self.tournament_selection(population, fitness_scores)
parent2 = self.tournament_selection(population, fitness_scores)
child1, child2 = self.crossover(parent1, parent2)
child1 = self.mutate(child1)
child2 = self.mutate(child2)
new_population.extend([child1, child2])
population = new_population[:self.pop_size]
return best_individual, best_fitness, history
def analyze_solution(self, solution, fitness):
"""分析最优解"""
print(f"\n=== 多机协同策略分析 ===")
print(f"总遮蔽时间: {fitness:.6f} 秒")
# FY1 (固定)
print(f"\nFY1 (固定最优解):")
print(f" 速度: {FY1_OPTIMAL['speed']:.1f} m/s")
print(f" 角度: {math.degrees(FY1_OPTIMAL['angle']):.2f}°")
print(f" 投放时间: {FY1_OPTIMAL['drop_time']:.3f}s")
print(f" 起爆延迟: {FY1_OPTIMAL['burst_delay']:.3f}s")
print(f" 起爆时间: {FY1_OPTIMAL['drop_time'] + FY1_OPTIMAL['burst_delay']:.3f}s")
# FY2
print(f"\nFY2 (优化结果):")
print(f" 速度: {solution[0]:.1f} m/s")
print(f" 角度: {math.degrees(solution[1]):.2f}°")
print(f" 投放时间: {solution[2]:.3f}s")
print(f" 起爆延迟: {solution[3]:.3f}s")
print(f" 起爆时间: {solution[2] + solution[3]:.3f}s")
# FY3
print(f"\nFY3 (优化结果):")
print(f" 速度: {solution[4]:.1f} m/s")
print(f" 角度: {math.degrees(solution[5]):.2f}°")
print(f" 投放时间: {solution[6]:.3f}s")
print(f" 起爆延迟: {solution[7]:.3f}s")
print(f" 起爆时间: {solution[6] + solution[7]:.3f}s")
# 协同分析
fy1_burst = FY1_OPTIMAL['drop_time'] + FY1_OPTIMAL['burst_delay']
fy2_burst = solution[2] + solution[3]
fy3_burst = solution[6] + solution[7]
print(f"\n=== 协同时序分析 ===")
print(f"FY1起爆: {fy1_burst:.3f}s → {fy1_burst + SMOKE_DURATION:.3f}s")
print(f"FY2起爆: {fy2_burst:.3f}s → {fy2_burst + SMOKE_DURATION:.3f}s")
print(f"FY3起爆: {fy3_burst:.3f}s → {fy3_burst + SMOKE_DURATION:.3f}s")
gap1 = fy2_burst - (fy1_burst + SMOKE_DURATION)
gap2 = fy3_burst - (fy2_burst + SMOKE_DURATION)
print(f"FY1→FY2间隙: {gap1:.3f}s {'(连续)' if abs(gap1) < 2 else '(有间隙)' if gap1 > 0 else '(重叠)'}")
print(f"FY2→FY3间隙: {gap2:.3f}s {'(连续)' if abs(gap2) < 2 else '(有间隙)' if gap2 > 0 else '(重叠)'}")
return solution
def main():
"""问题四主程序"""
print("问题四改进版: 基于FY1最优解的多机协同策略")
print("=" * 60)
print("约束条件:")
print(f"FY1参数固定: 速度{FY1_OPTIMAL['speed']}m/s, 角度{math.degrees(FY1_OPTIMAL['angle']):.2f}°")
print(f"优化目标: 最大化FY1+FY2+FY3的总遮蔽时间")
print("=" * 60)
start_time = time.time()
# 初始化
optimizer = ImprovedMultiUAVOptimizer()
# 创建约束遗传算法
ga = ConstrainedGeneticAlgorithm(optimizer)
# 执行优化
best_solution, best_fitness, history = ga.optimize()
total_time = time.time() - start_time
if best_fitness > 0:
print(f"\n优化成功完成！用时: {total_time:.1f}s")
result = ga.analyze_solution(best_solution, best_fitness)
# 评估结果
baseline_single = 1.39 # 问题一单弹结果
improvement = best_fitness / baseline_single
print(f"\n=== 性能评估 ===")
print(f"相比单弹基准提升: {improvement:.2f}倍")
print(f"三机协同效率: {best_fitness / 3:.2f}s/机")
if best_fitness > 10:
print("协同效果卓越!")
elif best_fitness > 6:
print("协同效果良好!")
elif best_fitness > 3:
print("协同效果一般，仍有提升空间")
else:
print("协同效果有限，需要重新分析策略")
return best_solution, best_fitness
else:
print("优化失败，未找到有效解")
return None, 0
if __name__ == "__main__":
solution, fitness = main()
if solution is not None:
print(f"\n最终结果:")
print(f"多机协同总遮蔽时间: {fitness:.6f}s")
print(f"策略特点: FY1固定最优 + FY2/FY3协同优化")
第五问
第五问没做出来，涉及的参数和变量实在太多了，我只有一些思路，代码纯ai跑的也算不出来答案，也算是本人今年一个遗憾吧，
针对问题五，问题五是综合前四个问题的精华，我们要优化多机多弹对多导弹的优化问题。因为第五问涉及四十个决策变量，十八个轨迹方程，在这里我们不采用广义的遗传算法与变量约束，我们通过物理直觉来定性分析，将每个变量控制在一个合理且较少的范围，这样可以在算法优化运算量情况下取得一个十分接近最优解的答案
类似第四问的优化思路，我们将所有无人机的方向指定为飞向原点（0,0,0）但留有一定角度的空间供我们细化，无人机的速度应当是FY1>FY5>FY2>FY4>FY3，因为在水平面xy投影下的位置，无人机对于假目标原点的水平距离是FY1>FY5>FY2>FY4>FY3，距离导弹近的应该飞的快一点早一点扔烟幕弹，距离远的应该飞的慢一点晚一点扔烟幕弹，这样能保证导弹在飞来的时候，在飞出前一个烟幕弹的遮蔽距离后能及时被后一个烟幕弹遮蔽，同时我们注意到三枚导弹的路径是在一条直线上的，我认为这提示我们在确定好无人机的方向和速度后，在引爆第一枚烟幕弹后，后两枚烟幕弹的引爆时间是有一定规律的，同时我们可以采用前几问已经得出的结论来进行问题五的分析
代码这里就不放了，实在想看看的可以下载我的论文看看
比赛参考了以下一些算法思路
公众号数学建模老哥给予了本蒟蒻莫大的帮助，在里面有很多很全的资料，不用费心整理了光学就行了，在此向hero们致敬，本文算法源代码都是从公众号里的小程序里的资料再学习的，想下载大家可以自己去关注下载，免费的，
如果文章涉及侵权请联系删除，这里只做分享学习用
梯度下降法（找多元函数最小值）
import matplotlib as mpl
import numpy as np
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']
def f(x):
return x[0] ** 2 + x[1] ** 2 + x[0] + x[0] + 1#目标函数：x² + y² + x + x + 1
start = [11.0, 2.0] # 初始点
lr = 0.1 # 学习率
max_iter = 1000 # 最大迭代次数
tol = 1e-6 # 可允许的误差
def grad(x):
dx = 2 * x[0] - 2 * x[1]
dy = 6 * x[1] - 2 * x[0]
return np.array([dx, dy])
def sd(start, lr, max_iter, tol):
x = np.array(start, dtype=np.float64)
path = [x.copy()]
for i in range(max_iter):
g = grad(x)
if np.linalg.norm(g) < tol:
break
x = x - lr * g
path.append(x.copy())
return x, path, i + 1
opt_point, path, iters = sd(start, lr, max_iter, tol)
opt_value = f(opt_point)
print(f"初始点: {start}")
print(f"下降次数: {iters}")
print(f"最优解: x = {opt_point[0]:.6f}, y = {opt_point[1]:.6f}")
print(f"最优值: f(x,y) = {opt_value:.6f}")
print(f"结果梯度模长: {np.linalg.norm(grad(opt_point)):.6e}")
#画图
if True:
x_vals = np.linspace(-1, 6, 100)
y_vals = np.linspace(-1, 6, 100)
X, Y = np.meshgrid(x_vals, y_vals)
Z = f([X, Y])
plt.figure(figsize=(8, 6))
plt.contour(X, Y, Z, levels=30, cmap='viridis')
plt.colorbar(label='函数值')
path_arr = np.array(path)
plt.plot(path_arr[:, 0], path_arr[:, 1], 'ro-', markersize=5, label='迭代路径')
plt.plot(start[0], start[1], 'go', markersize=8, label='初始点')
plt.plot(opt_point[0], opt_point[1], 'bo', markersize=8, label='最优解')
plt.xlabel('x')
plt.ylabel('y')
plt.title('最速下降法迭代路径')
plt.legend()
plt.grid(True)
plt.show()
秩比综合评价法
将效益型指标从小到大排序进行排名、成本型指标从大到小排序进行排名，再计算秩和比，最后统计回归、分档排序
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
秩和比综合评价法案例说明：
本案例对8所高校的综合实力进行评价，考虑3项指标：
1. 科研经费(万元) 2. 师资力量(教授人数) 3. 就业率(%)
所有指标均为效益型指标（数值越大越好），通过秩和比法计算RSR值，
并进行分档（1-4档，1档最优）和排序，评价各高校的综合实力。
"""
import numpy as np
# 原始数据矩阵（8个样本×3个指标）
# 行：高校1-8，列：科研经费、师资力量、就业率
data = np.array([
[8500, 90, 78], # 高校1
[9200, 85, 90], # 高校2
[7800, 92, 86], # 高校3
[8800, 88, 94], # 高校4
[9000, 82, 89], # 高校5
[8200, 87, 91], # 高校6
[8600, 91, 83], # 高校7
[8900, 84, 87] # 高校8
])
# 各指标权重
weights = np.array([0.3, 0.4, 0.3])
# 指标类型（True表示效益型，False表示成本型）
is_benefit = np.array([True, True, True])
m, n = data.shape # m=样本数, n=指标数
# 1. 编秩（对每个指标进行排序并赋值秩次）
R = np.zeros((m, n)) # 秩# 矩阵
for j in range(n):
# 按降序排序并获取索引（从大到小）
idx = np.argsort(-data[:, j])
r = np.zeros(m)
for i in range(m):
r[idx[i]] = i + 1 # 秩从1到m（最大值秩为1）
# 对于成本型指标，反转秩次（最小值秩为1）
if not is_benefit[j]:
r = m + 1 - r
R[:, j] = r
# 2. 计算加权秩和比RSR
rsr = np.dot(R, weights) / m # 加权求和后除以样本数
# 3. 排序（从优到劣）
sorted_indices = np.argsort(-rsr)
rank = np.zeros(m, dtype=int)
for i in range(m):
rank[sorted_indices[i]] = i + 1 # 排名从1开始
# 4. 分档（分为4档，1档最优）
grade = np.zeros(m, dtype=int)
quarter = np.ceil(m / 4) # 每档的大致数量
grade[sorted_indices[:int(quarter)]] = 1
grade[sorted_indices[int(quarter):int(2*quarter)]] = 2
grade[sorted_indices[int(2*quarter):int(3*quarter)]] = 3
grade[sorted_indices[int(3*quarter):]] = 4
# 输出结果
print("秩和比综合评价法计算结果：")
print(f"样本数: {m}, 指标数: {n}")
print(f"各指标权重: {weights.round(4)}")
print("\n各高校秩和比(RSR)值、排名及分档：")
for i in range(m):
print(f"高校{i+1}: RSR值={rsr[i]:.4f}, 排名={rank[i]}, 档级={grade[i]}")
整数规划（决策变量要求是整数的问题）
整数规划主要有三种
1. 纯整数规划：指全部决策变量都必须取整数的整数规划。
2. 混合整数规划：指决策变量有一部分必须取整数值的整数规划。
3. 0-1 整数规划：指决策变量只能取值 0 或 1 的整数规划。
这里给一个0-1规划的例子
import pulp as pl
# 创建问题实例：最大化总预期收益
prob = pl.LpProblem("投资组合问题", pl.LpMaximize)
# 定义决策变量：5个项目是否投资（0-1整数）
projects = range(5)
x = pl.LpVariable.dicts("投资项目", projects, cat=pl.LpBinary)
# 项目数据：[初始投资(万元), 预期收益(万元)]
project_data = {
0: [100, 150],
1: [120, 180],
2: [80, 120],
3: [150, 220],
4: [90, 130]
}
# 定义目标函数：最大化总预期收益
prob += sum(project_data[i][1] * x[i] for i in projects), "总预期收益"
# 定义约束条件：总投资不超过预算
prob += sum(project_data[i][0] * x[i] for i in projects) <= 300, "预算约束"
# 求解问题
prob.solve()
# 输出结果
print(f"求解状态: {pl.LpStatus[prob.status]}")
print("\n最优投资组合:")
total_invest = 0
total_profit = 0
for i in projects:
if pl.value(x[i]) == 1:
print(f"投资项目{i}: 投资{project_data[i][0]}万元, 预期收益{project_data[i][1]}万元")
total_invest += project_data[i][0]
total_profit += project_data[i][1]
print(f"\n总投资: {total_invest}万元, 总预期收益: {total_profit}万元")
print(f"预算剩余: {300 - total_invest}万元")
遗传算法（存在大量局部最优但是要全局最优的问题）
简单来说就是复杂曲面函数求最值
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
遗传算法案例说明：
本案例使用遗传算法求解Griewank函数的最小值。
Griewank函数是一个具有大量局部最优解的多峰函数，全局最小值为0，位于(0,0,...,0)。
遗传算法模拟生物进化过程，通过选择、交叉和变异操作寻找最优解，
适合求解复杂的全局优化问题。
"""
import numpy as np
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']
# Griewank函数：一个复杂的多峰函数
def griewank(x):
sum_part = sum([xi**2 for xi in x]) / 4000
prod_part = np.prod([np.cos(xi / np.sqrt(i+1)) for i, xi in enumerate(x)])
return sum_part - prod_part + 1
# 初始化种群
def initialize_population(pop_size, dim, bounds):
population = []
for _ in range(pop_size):
individual = [np.random.uniform(low, high) for low, high in bounds]
population.append(individual)
return np.array(population)
# 选择操作（轮盘赌选择）
def select(population, fitness, num_parents):
# 适应度越小越好，转换为选择概率（取倒数）
fitness = np.max(fitness) - fitness + 1e-10 # 确保非负
probabilities = fitness / np.sum(fitness)
# 选择父代
parents = []
for _ in range(num_parents):
idx = np.random.choice(len(population), p=probabilities)
parents.append(population[idx])
return np.array(parents)
# 交叉操作（单点交叉）
def crossover(parents, offspring_size):
offspring = []
crossover_rate = 0.8 # 交叉概率
for i in range(offspring_size[0]):
# 随机选择两个父代
parent1_idx = i % len(parents)
parent2_idx = (i + 1) % len(parents)
parent1 = parents[parent1_idx]
parent2 = parents[parent2_idx]
# 以一定概率进行交叉
if np.random.rand() < crossover_rate:
# 随机选择交叉点
crossover_point = np.random.randint(1, len(parent1))
# 生成子代
child = np.concatenate([parent1[:crossover_point], parent2[crossover_point:]])
offspring.append(child)
else:
# 不交叉，直接复制父代
offspring.append(parent1.copy())
return np.array(offspring)
# 变异操作
def mutate(offspring, bounds, mutation_rate):
for i in range(len(offspring)):
for j in range(len(offspring[i])):
# 以一定概率进行变异
if np.random.rand() < mutation_rate:
# 高斯变异
offspring[i][j] += np.random.normal(0, 0.5)
# 边界处理
offspring[i][j] = max(bounds[j][0], min(offspring[i][j], bounds[j][1]))
return offspring
# 遗传算法
def genetic_algorithm(objective_func, bounds, pop_size, num_generations):
dim = len(bounds)
# 初始化种群
population = initialize_population(pop_size, dim, bounds)
# 评估初始种群
fitness = np.array([objective_func(ind) for ind in population])
# 记录最优解
best_idx = np.argmin(fitness)
best_solution = population[best_idx].copy()
best_value = fitness[best_idx]
best_history = [best_value]
# 主循环
for gen in range(num_generations):
# 选择父代
num_parents = pop_size // 2
parents = select(population, fitness, num_parents)
# 交叉产生子代
offspring_size = (pop_size - num_parents, dim)
offspring = crossover(parents, offspring_size)
# 变异
mutation_rate = 0.1 # 变异概率
offspring = mutate(offspring, bounds, mutation_rate)
# 形成新种群
population = np.concatenate([parents, offspring])
# 评估新种群
fitness = np.array([objective_func(ind) for ind in population])
# 更新最优解
current_best_idx = np.argmin(fitness)
current_best_value = fitness[current_best_idx]
if current_best_value < best_value:
best_solution = population[current_best_idx].copy()
best_value = current_best_value
# 记录历史
best_history.append(best_value)
return best_solution, best_value, best_history
# 参数设置
bounds = [(-600, 600), (-600, 600)] # 变量范围
pop_size = 50 # 种群大小
num_generations = 100 # 进化代数
# 运行遗传算法
best_solution, best_value, best_history = genetic_algorithm(
griewank, bounds, pop_size, num_generations
)
# 输出结果
print("遗传算法求解Griewank函数结果：")
print(f"最优解: x = {best_solution[0]:.6f}, y = {best_solution[1]:.6f}")
print(f"最优值: f(x,y) = {best_value:.6f}")
print(f"理论最优解: (0, 0)，最优值: 0")
# 可视化结果（可选）
if True:
# 绘制函数曲面
fig = plt.figure(figsize=(15, 6))
# 3D曲面图
ax1 = fig.add_subplot(121, projection='3d')
x = np.linspace(bounds[0][0], bounds[0][1], 50)
y = np.linspace(bounds[1][0], bounds[1][1], 50)
X, Y = np.meshgrid(x, y)
Z = np.array([griewank([x, y]) for x, y in zip(np.ravel(X), np.ravel(Y))]).reshape(X.shape)
ax1.plot_surface(X, Y, Z, cmap='viridis', alpha=0.7)
ax1.scatter(best_solution[0], best_solution[1], best_value, color='red', s=100, label='最优解')
ax1.set_xlabel('x')
ax1.set_ylabel('y')
ax1.set_zlabel('f(x,y)')
ax1.set_title('Griewank函数曲面与最优解')
ax1.legend()
# 收敛曲线图
ax2 = fig.add_subplot(122)
ax2.plot(best_history)
ax2.set_xlabel('进化代数')
ax2.set_ylabel('最优值')
ax2.set_title('收敛曲线')
ax2.grid(True)
plt.tight_layout()
plt.show()
箱型图异常值检测
把数据四等分，每一部分包含25%的数据点。
把数据集中间50%的部分作为正常的数据，设置一个阈值，超出阈值的部分检测为异常值
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
箱型图检测异常值案例说明：
本案例使用箱型图方法检测数据中的异常值。箱型图基于数据的四分位数，
通过计算上下限（Q1-1.5*IQR和Q3+1.5*IQR）来识别异常值，其中IQR是四分位距。
示例中生成包含异常值的随机数据，使用箱型图进行可视化并标记异常值。
"""
import numpy as np
import matplotlib.pyplot as plt
import pandas as pd
plt.rcParams['font.sans-serif'] = ['Microsoft YaHei']
# 生成示例数据：主要分布在[10, 50]之间，添加一些异常值
np.random.seed(42) # 设置随机种子，确保结果可复现
normal_data = np.random.normal(loc=30, scale=8, size=200) # 正常数据
outliers = np.array([5, 60, 65, 70, -2, 75]) # 异常值
data = np.concatenate([normal_data, outliers]) # 合并数据
# 计算四分位数和异常值边界
q1 = np.percentile(data, 25) # 第一四分位数
q3 = np.percentile(data, 75) # 第三四分位数
iqr = q3 - q1 # 四分位距
lower_bound = q1 - 1.5 * iqr # 下限
upper_bound = q3 + 1.5 * iqr # 上限
# 检测异常值
outlier_indices = np.where((data < lower_bound) | (data > upper_bound))[0]
outlier_values = data[outlier_indices]
# 输出结果
print("箱型图异常值检测结果：")
print(f"数据总量: {len(data)} 个")
print(f"第一四分位数(Q1): {q1:.2f}")
print(f"第三四分位数(Q3): {q3:.2f}")
print(f"四分位距(IQR): {iqr:.2f}")
print(f"异常值下限: {lower_bound:.2f}")
print(f"异常值上限: {upper_bound:.2f}")
print(f"检测到异常值数量: {len(outlier_values)} 个")
print(f"异常值: {', '.join([f'{x:.2f}' for x in outlier_values])}")
# 可视化箱型图
plt.figure(figsize=(10, 6))
box = plt.boxplot(data, patch_artist=True,
boxprops=dict(facecolor='lightblue', color='blue'),
capprops=dict(color='blue'),
whiskerprops=dict(color='blue'),
flierprops=dict(marker='o', color='red', markersize=8),
medianprops=dict(color='green', linewidth=2))
# 添加文本说明
plt.text(1.1, q1, f'Q1: {q1:.2f}', verticalalignment='center')
plt.text(1.1, q3, f'Q3: {q3:.2f}', verticalalignment='center')
plt.text(1.1, lower_bound, f'下限: {lower_bound:.2f}', verticalalignment='center')
plt.text(1.1, upper_bound, f'上限: {upper_bound:.2f}', verticalalignment='center')
plt.title('箱型图异常值检测')
plt.ylabel('数据值')
plt.grid(axis='y', linestyle='--', alpha=0.7)
plt.show()
线性规划（线性目标函数的极值问题）
给的源代码访问约束条件上界的方式不正确，这个修改后的可正常运行
#!/usr/bin/env python
# -*- coding: utf-8 -*-
"""
线性规划模型案例说明：
本案例使用线性规划解决生产资源分配问题。某工厂生产两种产品A和B，
每种产品需要消耗原材料和工时，有最大供应量和市场需求限制。
目标是最大化总利润。线性规划适用于目标函数和约束条件均为线性的问题。
"""
import pulp as pl
# 创建问题实例：最大化总利润
prob = pl.LpProblem("生产资源分配问题", pl.LpMaximize)
# 定义决策变量：产品A和B的生产量（非负）
x1 = pl.LpVariable("产品A产量", lowBound=0)
x2 = pl.LpVariable("产品B产量", lowBound=0)
# 定义目标函数：最大化总利润（A利润50元/件，B利润60元/件）
prob += 50 * x1 + 60 * x2, "总利润"
# 定义约束条件
prob += 2 * x1 + 3 * x2 <= 100, "原材料约束" # 原材料最大供应量100单位
prob += 4 * x1 + 2 * x2 <= 120, "工时约束" # 总工时不超过120小时
prob += x1 <= 25, "产品A需求约束" # 产品A最大需求25件
prob += x2 <= 30, "产品B需求约束" # 产品B最大需求30件
# 求解问题
prob.solve()
# 输出结果
print(f"求解状态: {pl.LpStatus[prob.status]}")
print(f"最优解: 产品A生产 {pl.value(x1):.2f} 件, 产品B生产 {pl.value(x2):.2f} 件")
print(f"最大总利润: {pl.value(prob.objective):.2f} 元")
# 输出各约束条件的使用情况
print("\n约束条件使用情况:")
for name, constraint in prob.constraints.items():
actual_value = constraint.value()
# 获取约束的右端项（限制值）
# 对于不等式约束，右端常数项就是限制值
rhs = -constraint.constant # 注意：PuLP中常数项存储为负数形式
print(f"{name}: 实际使用 {actual_value:.2f}, 限制 {rhs:.2f}")
DAMO开发者矩阵，由阿里巴巴达摩院和中国互联网协会联合发起，致力于探讨最前沿的技术趋势与应用成果，搭建高质量的交流与分享平台，推动技术创新与产业应用链接，围绕“人工智能与新型计算”构建开放共享的开发者生态。
更多推荐
5
0- 0
已为社区贡献4条内容
所有评论(0)