242 lines
8.8 KiB
Python
242 lines
8.8 KiB
Python
import numpy as np
|
||
import pandas as pd
|
||
from scipy.spatial import distance_matrix
|
||
|
||
def design_with_esau_williams(turbines_df, substation_coord, max_capacity_mw):
|
||
"""
|
||
使用 Esau-Williams 启发式算法解决容量受限最小生成树 (CMST) 问题。
|
||
|
||
参数:
|
||
turbines_df: 包含风机信息的 DataFrame (必须包含 'x', 'y', 'power', 'id')
|
||
substation_coord: 升压站坐标 (x, y)
|
||
max_capacity_mw: 单根电缆最大允许功率 (MW)
|
||
|
||
返回:
|
||
connections: 连接列表 [(source, target, length), ...]
|
||
turbines_with_cluster: 带有 'cluster' 列的 turbines DataFrame (用于兼容性)
|
||
"""
|
||
|
||
# 数据准备
|
||
n_turbines = len(turbines_df)
|
||
coords = turbines_df[['x', 'y']].values
|
||
powers = turbines_df['power'].values
|
||
ids = turbines_df['id'].values
|
||
|
||
# 升压站坐标
|
||
if substation_coord.ndim > 1:
|
||
sx, sy = substation_coord[0][0], substation_coord[0][1]
|
||
else:
|
||
sx, sy = substation_coord[0], substation_coord[1]
|
||
|
||
# 1. 计算距离矩阵
|
||
# 风机到风机
|
||
dist_matrix = distance_matrix(coords, coords)
|
||
# 风机到升压站
|
||
dists_to_sub = np.sqrt((coords[:, 0] - sx)**2 + (coords[:, 1] - sy)**2)
|
||
|
||
# 2. 初始化组件 (Components)
|
||
# 初始状态下，每个风机是一个独立的组件，直接连接到升压站
|
||
# 为了方便查找，我们维护一个 components 字典
|
||
# key: component_root_id (代表该组件的唯一标识)
|
||
# value: {
|
||
# 'members': {node_idx, ...},
|
||
# 'total_power': float,
|
||
# 'gate_node': int (连接到升压站的节点索引),
|
||
# 'gate_cost': float (gate_node 到升压站的距离)
|
||
# }
|
||
|
||
components = {}
|
||
for i in range(n_turbines):
|
||
components[i] = {
|
||
'members': {i},
|
||
'total_power': powers[i],
|
||
'gate_node': i,
|
||
'gate_cost': dists_to_sub[i]
|
||
}
|
||
|
||
# 记录已经建立的连接 (不包括通往升压站的默认连接)
|
||
# 格式: (u, v, length)
|
||
established_edges = []
|
||
# 记录已建立连接的坐标，用于交叉检查: [((x1, y1), (x2, y2)), ...]
|
||
established_lines = []
|
||
|
||
def do_intersect(p1, p2, p3, p4):
|
||
"""
|
||
检测线段 (p1, p2) 和 (p3, p4) 是否严格相交 (不包括端点接触)
|
||
"""
|
||
# 检查是否共享端点
|
||
if (p1[0]==p3[0] and p1[1]==p3[1]) or (p1[0]==p4[0] and p1[1]==p4[1]) or \
|
||
(p2[0]==p3[0] and p2[1]==p3[1]) or (p2[0]==p4[0] and p2[1]==p4[1]):
|
||
return False
|
||
|
||
def ccw(A, B, C):
|
||
# 向量叉积
|
||
return (C[1]-A[1]) * (B[0]-A[0]) - (B[1]-A[1]) * (C[0]-A[0])
|
||
|
||
# 如果跨立实验符号相反，则相交
|
||
d1 = ccw(p1, p2, p3)
|
||
d2 = ccw(p1, p2, p4)
|
||
d3 = ccw(p3, p4, p1)
|
||
d4 = ccw(p3, p4, p2)
|
||
|
||
# 严格相交判断 (忽略共线重叠的情况，视为不交叉)
|
||
if ((d1 > 1e-9 and d2 < -1e-9) or (d1 < -1e-9 and d2 > 1e-9)) and \
|
||
((d3 > 1e-9 and d4 < -1e-9) or (d3 < -1e-9 and d4 > 1e-9)):
|
||
return True
|
||
|
||
return False
|
||
|
||
# 3. 迭代优化
|
||
while True:
|
||
# 预先收集当前所有组件的 Gate Edges (连接升压站的线段)
|
||
# 格式: {cid: (gate_node_coord, substation_coord)}
|
||
current_gate_lines = {}
|
||
sub_coord_tuple = (sx, sy)
|
||
for cid, data in components.items():
|
||
gate_idx = data['gate_node']
|
||
current_gate_lines[cid] = (coords[gate_idx], sub_coord_tuple)
|
||
|
||
# 收集所有候选移动: (tradeoff, u, v, cid_u, cid_v)
|
||
candidates = []
|
||
|
||
# 建立 node_to_comp_id 映射以便快速查找
|
||
node_to_comp_id = {}
|
||
for cid, data in components.items():
|
||
for member in data['members']:
|
||
node_to_comp_id[member] = cid
|
||
|
||
# 遍历所有边 (i, j)
|
||
for i in range(n_turbines):
|
||
cid_i = node_to_comp_id[i]
|
||
gate_cost_i = components[cid_i]['gate_cost']
|
||
|
||
for j in range(n_turbines):
|
||
if i == j: continue
|
||
|
||
cid_j = node_to_comp_id[j]
|
||
|
||
# 必须是不同组件
|
||
if cid_i == cid_j:
|
||
continue
|
||
|
||
# 检查容量约束
|
||
if components[cid_i]['total_power'] + components[cid_j]['total_power'] > max_capacity_mw:
|
||
continue
|
||
|
||
# 计算 Tradeoff
|
||
dist_ij = dist_matrix[i, j]
|
||
tradeoff = dist_ij - gate_cost_i
|
||
|
||
# 只有当 tradeoff < 0 时，合并才是有益的
|
||
if tradeoff < -1e-9:
|
||
candidates.append((tradeoff, i, j, cid_i, cid_j))
|
||
|
||
# 按 tradeoff 排序 (从小到大，越小越好)
|
||
candidates.sort(key=lambda x: x[0])
|
||
|
||
best_move = None
|
||
|
||
# 延迟检测: 从最好的开始检查交叉
|
||
for cand in candidates:
|
||
tradeoff, u, v, cid_u, cid_v = cand
|
||
|
||
p_u = coords[u]
|
||
p_v = coords[v]
|
||
|
||
# 快速包围盒测试 (AABB) 准备
|
||
min_x_uv, max_x_uv = min(p_u[0], p_v[0]), max(p_u[0], p_v[0])
|
||
min_y_uv, max_y_uv = min(p_u[1], p_v[1]), max(p_u[1], p_v[1])
|
||
|
||
is_crossing = False
|
||
|
||
# 1. 检查与已固定的内部边的交叉
|
||
for line in established_lines:
|
||
p_a, p_b = line[0], line[1]
|
||
|
||
if max(p_a[0], p_b[0]) < min_x_uv or min(p_a[0], p_b[0]) > max_x_uv or \
|
||
max(p_a[1], p_b[1]) < min_y_uv or min(p_a[1], p_b[1]) > max_y_uv:
|
||
continue
|
||
|
||
if do_intersect(p_u, p_v, p_a, p_b):
|
||
is_crossing = True
|
||
break
|
||
|
||
if is_crossing:
|
||
continue
|
||
|
||
# 2. 检查与所有活跃 Gate Edges 的交叉 (排除被移除的那个)
|
||
# 正在合并 cid_u -> cid_v，意味着 cid_u 的 Gate 将被移除。
|
||
# 但 cid_v 的 Gate 以及其他所有组件的 Gate 仍然存在。
|
||
for cid, gate_line in current_gate_lines.items():
|
||
if cid == cid_u:
|
||
continue # 这个 Gate 即将移除，不构成障碍
|
||
|
||
p_a, p_b = gate_line[0], gate_line[1]
|
||
|
||
if max(p_a[0], p_b[0]) < min_x_uv or min(p_a[0], p_b[0]) > max_x_uv or \
|
||
max(p_a[1], p_b[1]) < min_y_uv or min(p_a[1], p_b[1]) > max_y_uv:
|
||
continue
|
||
|
||
if do_intersect(p_u, p_v, p_a, p_b):
|
||
is_crossing = True
|
||
break
|
||
|
||
if not is_crossing:
|
||
best_move = (u, v, cid_u, cid_v)
|
||
break
|
||
|
||
# 如果没有找到有益的合并，或者所有可行合并都会增加成本，则停止
|
||
if best_move is None:
|
||
break
|
||
|
||
# 执行合并
|
||
u, v, cid_u, cid_v = best_move
|
||
|
||
# 将 cid_u 并入 cid_v
|
||
# 1. 记录新边
|
||
established_edges.append((u, v, dist_matrix[u, v]))
|
||
established_lines.append((coords[u], coords[v]))
|
||
|
||
# 2. 更新组件信息
|
||
comp_u = components[cid_u]
|
||
comp_v = components[cid_v]
|
||
|
||
# 合并成员
|
||
comp_v['members'].update(comp_u['members'])
|
||
comp_v['total_power'] += comp_u['total_power']
|
||
|
||
# Gate 节点和 Cost 保持 cid_v 的不变
|
||
# (因为我们将 U 接到了 V 上，U 的原 gate 被移除，V 的 gate 仍是通往升压站的路径)
|
||
|
||
# 3. 删除旧组件
|
||
del components[cid_u]
|
||
|
||
# 4. 构建最终结果
|
||
connections = []
|
||
|
||
# 添加内部边 (风机间)
|
||
for u, v, length in established_edges:
|
||
source = f'turbine_{u}'
|
||
target = f'turbine_{v}'
|
||
connections.append((source, target, length))
|
||
|
||
# 添加 Gate 边 (风机到升压站)
|
||
# 此时 components 中剩下的每个组件都有一个 gate_node 连接到 Substation
|
||
cluster_mapping = {} # node_id -> cluster_id (0..N-1)
|
||
|
||
for idx, (cid, data) in enumerate(components.items()):
|
||
gate_node = data['gate_node']
|
||
gate_cost = data['gate_cost']
|
||
|
||
connections.append((f'turbine_{gate_node}', 'substation', gate_cost))
|
||
|
||
# 记录 cluster id
|
||
for member in data['members']:
|
||
cluster_mapping[member] = idx
|
||
|
||
# 更新 turbines DataFrame
|
||
turbines_with_cluster = turbines_df.copy()
|
||
turbines_with_cluster['cluster'] = turbines_with_cluster['id'].map(cluster_mapping)
|
||
|
||
return connections, turbines_with_cluster
|