# 题目重述
在边长为160公里的正方形空域（顶点为$(0,0)$、$(160,0)$、$(160,160)$、$(0,160)$）中，现有5架飞机正在飞行，另有一架新飞机从$(0,0)$以方向角$52^\circ$进入。所有飞机速度均为800 km/h。
要求：
1. 任意两架飞机之间的最小距离必须大于8 km；
2. 每架飞机的方向角调整幅度不超过$30^\circ$；
3. 新进入飞机在边界时与区域内飞机距离大于60 km（已满足）；
4. 最多考虑6架飞机；
5. 不考虑飞机离开区域后的情况；
6. 输出调整后的方向角，误差不超过$0.01^\circ$，且总调整量尽可能小。
请建立数学模型并使用群智能优化算法求解。
---
# 详解
本问题是一个典型的**动态避障约束优化问题**，目标是最小化方向角调整总量，在满足飞行安全的前提下实现无碰撞航路规划。
## 1. 数学建模
### （1）变量定义
设第$i$架飞机的调整量为$\Delta\theta_i$，则其实际飞行方向为：
$$
\theta_i' = \theta_i + \Delta\theta_i, \quad \text{其中 } |\Delta\theta_i| \leq 30^\circ
$$
### （2）位置函数（时间$t$小时后）
$$
x_i(t) = x_i + 800t \cdot \cos(\theta_i'\cdot\pi/180)
$$
$$
y_i(t) = y_i + 800t \cdot \sin(\theta_i'\cdot\pi/180)
$$
### （3）两机间距
对任意$i < j$，有：
$$
d_{ij}(t) = \sqrt{(x_i(t)-x_j(t))^2 + (y_i(t)-y_j(t))^2}
$$
我们需保证在预测时间段内（如$[0, T]$，取$T=1$小时），所有$d_{ij}(t) > 8$
### （4）目标函数
最小化总调整量（使用绝对值之和）：
$$
\min \sum_{i=1}^{6} |\Delta\theta_i|
$$
### （5）约束条件
- $|\Delta\theta_i| \leq 30^\circ$
- $\forall i<j, \forall t \in [0,1],\ d_{ij}(t) > 8$
- 初始时刻新飞机与其他飞机距离 $>60$ km（原始数据满足）
由于连续时间难以遍历，采用**离散化采样**：将$[0,1]$分为60个点（每分钟一次）
---
## 2. 解法选择：粒子群优化算法（PSO）
PSO适合解决此类高维、非线性、带约束的优化问题。每个粒子表示一组方向角调整方案。
---
# 代码实现（Python + NumPy + Matplotlib 可视化）
```python
import numpy as np
import matplotlib.pyplot as plt
from matplotlib.patches import Circle
# ================= 参数设置 =================
N_PARTICLES = 50 # 粒子数量
MAX_ITER = 200 # 最大迭代次数
DIM = 6 # 维度：6架飞机的角度调整
V_MIN, V_MAX = -5, 5 # 粒子速度范围
UB = 30 # 上界
LB = -30 # 下界
W = 0.7 # 惯性权重
C1 = 1.5 # 个体学习因子
C2 = 1.5 # 群体学习因子
T_STEPS = 60 # 时间步数（模拟未来1小时，每分钟一步）
PENALTY = 1e5 # 约束违反惩罚项
# 初始数据 [x, y, theta_deg]
planes = np.array([
[150, 40, 243],
[85, 85, 236],
[150, 155, 220.5],
[145, 50, 159],
[130, 150, 230],
[0, 0, 52] # 新进入飞机
])
speed = 800 # km/h
def simulate_trajectory(delta_theta, t_hours):
"""
计算所有飞机在t时刻的位置
delta_theta: (6,) 调整量数组
t_hours: float 或 array of floats
返回: (6, 2) 位置矩阵
"""
theta_prime = np.radians(planes[:, 2] + delta_theta)
cos_t = np.cos(theta_prime)
sin_t = np.sin(theta_prime)
x = planes[:, 0][:, None] + speed * t_hours * cos_t[:, None]
y = planes[:, 1][:, None] + speed * t_hours * sin_t[:, None]
return np.stack((x, y), axis=-1)
def fitness(individual):
"""
适应度函数
individual: (6,) 当前调整量
返回: scalar fitness value (越小越好)
"""
# 检查是否超出调整限制
if np.any(individual < LB) or np.any(individual > UB):
return np.inf
# 时间序列: 0 到 1 小时，共 T_STEPS 步
t_vec = np.linspace(0, 1, T_STEPS)
# 获取所有时刻所有飞机的位置
traj = simulate_trajectory(individual, t_vec) # shape: (6, T_STEPS, 2)
min_dist = np.inf
for i in range(6):
for j in range(i+1, 6):
dist = np.sqrt(np.sum((traj[i] - traj[j])**2, axis=1)) # 各时间点距离
min_dist = min(min_dist, np.min(dist))
# 目标：最小化调整量；约束：最小距离 > 8
obj = np.sum(np.abs(individual))
if min_dist <= 8:
penalty = PENALTY * (8 - min_dist)**2
else:
penalty = 0
return obj + penalty
# ================= PSO 主循环 =================
# 初始化粒子群
np.random.seed(42)
X = np.random.uniform(LB, UB, (N_PARTICLES, DIM)) # 位置
V = np.random.uniform(V_MIN, V_MAX, (N_PARTICLES, DIM)) # 速度
pbest_pos = X.copy()
pbest_fit = np.array([fitness(X[i]) for i in range(N_PARTICLES)])
gbest_idx = np.argmin(pbest_fit)
gbest_pos = pbest_pos[gbest_idx].copy()
gbest_fit = pbest_fit[gbest_idx]
# 迭代优化
for iter in range(MAX_ITER):
for i in range(N_PARTICLES):
# 更新速度
r1, r2 = np.random.rand(DIM), np.random.rand(DIM)
V[i] = W*V[i] + C1*r1*(pbest_pos[i] - X[i]) + C2*r2*(gbest_pos - X[i])
V[i] = np.clip(V[i], V_MIN, V_MAX)
# 更新位置
X[i] += V[i]
X[i] = np.clip(X[i], LB, UB)
# 更新个体最优
fit = fitness(X[i])
if fit < pbest_fit[i]:
pbest_pos[i] = X[i].copy()
pbest_fit[i] = fit
# 更新全局最优
if fit < gbest_fit:
gbest_pos = X[i].copy()
gbest_fit = fit
if iter % 50 == 0:
print(f"Iter {iter}, Best Fitness: {gbest_fit:.4f}")
# 输出最终结果
final_adjustment = np.round(gbest_pos, 2)
final_angles = planes[:, 2] + final_adjustment
print("\n=== 最终调整方案 ===")
print("编号\t原方向角\t调整量\t\t新方向角")
for idx in range(6):
print(f"{idx+1}\t{planes[idx,2]:.1f}\t\t{final_adjustment[idx]:+.2f}\t\t{final_angles[idx]:.2f}")
# 验证最小距离
def check_min_distance(delta_theta):
t_vec = np.linspace(0, 1, T_STEPS)
traj = simulate_trajectory(delta_theta, t_vec)
min_d = np.inf
for i in range(6):
for j in range(i+1, 6):
d = np.sqrt(np.sum((traj[i] - traj[j])**2, axis=1))
min_d = min(min_d, np.min(d))
return min_d
min_separation = check_min_distance(final_adjustment)
print(f"\n所有飞机间的最小距离: {min_separation:.2f} km (>8 km 才安全)")
assert min_separation > 8, "❌ 存在碰撞风险！"
print("✅ 安全飞行方案已找到")
# ================= 可视化轨迹 =================
def plot_trajectories(delta_theta):
plt.figure(figsize=(10, 10))
t_vec = np.linspace(0, 1, 60)
traj = simulate_trajectory(delta_theta, t_vec)
colors = ['r', 'g', 'b', 'c', 'm', 'orange']
labels = [f'Plane {i+1}' for i in range(6)]
for i in range(6):
x_traj = traj[i, :, 0]
y_traj = traj[i, :, 1]
plt.plot(x_traj, y_traj, c=colors[i], label=labels[i], linewidth=2)
plt.scatter(x_traj[0], y_traj[0], c=colors[i], s=50, zorder=5)
# 添加方向箭头
dx = x_traj[5] - x_traj[0]
dy = y_traj[5] - y_traj[0]
plt.arrow(x_traj[0], y_traj[0], dx*2, dy*2, head_width=3, fc=colors[i], ec=colors[i])
# 绘制安全圆（可选）
for i in range(6):
circle = Circle((traj[i, 0, 0], traj[i, 0, 1]), radius=4, color=colors[i], alpha=0.2)
plt.gca().add_patch(circle)
plt.xlim(0, 160)
plt.ylim(0, 160)
plt.xlabel('X (km)')
plt.ylabel('Y (km)')
plt.title('Flight Paths After Adjustment\n(Arrows show initial direction)')
plt.legend()
plt.grid(True, linestyle='--', alpha=0.5)
plt.axis('equal')
plt.tight_layout()
plt.show()
# 是否显示可视化？
# plot_trajectories(final_adjustment)
```
---
# 知识点
- **粒子群优化算法（PSO）**：一种基于群体寻优的启发式算法，通过跟踪个体与全局最优更新粒子状态。
- **运动学建模**：利用三角函数构建匀速直线运动轨迹，实现飞机位置随时间演化。
- **约束处理与罚函数法**：将硬约束转化为目标函数中的惩罚项，便于智能算法求解。
扫一扫
448
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
