diff --git a/README.md b/README.md index da4e0d2..e6b6697 100644 --- a/README.md +++ b/README.md @@ -1,85 +1,73 @@ -# 海上风电场集电线路设计优化工具 +# 海上风电场集电系统设计工具 -## 项目简介 +一个用于设计和优化海上风电场集电系统的Python工具，支持多种布局算法和电缆优化方案。 -这是一个用于海上风电场集电线路拓扑设计和优化的Python工具。它专注于解决大规模海上风电场的集电系统规划问题，通过算法比较不同设计方案的经济性和技术指标。 +## 功能特性 -本项目特别针对**海上风电**场景进行了优化，考虑了海缆的高昂成本、大功率风机（6-10MW）以及严格的电缆载流量约束。 +- 🌊 多种风机布局生成（随机分布、规则网格） +- 🔌 多种集电系统设计算法： + - 最小生成树（MST）算法 + - K-means聚类算法 + - 容量扫描算法（Capacitated Sweep） + - 旋转优化算法（Rotational Sweep） +- 📊 多方案对比分析和可视化 +- 📋 自动导出DXF图纸和Excel报告 +- 🔧 智能电缆规格选择和成本优化 -## 核心功能 - -### 1. 多种布局生成与导入 -- **自动生成**：支持生成规则的矩阵式（Grid）风机布局，模拟海上风电场常见排布。 -- **Excel导入**：支持从 `coordinates.xlsx` 导入自定义的风机和升压站坐标。 - - 格式要求：包含 `Type` (Turbine/Substation), `ID`, `X`, `Y`, `Power` 列。 - -### 2. 智能拓扑优化算法 -- **最小生成树 (MST)**： - - 计算全局最短路径长度。 - - *注意*：在大规模风电场中，纯MST往往会导致根部电缆严重过载，仅作为理论最短路径参考。 -- **扇区聚类 (Angular K-means)**： - - **无交叉设计**：基于角度（扇区）进行聚类，从几何上杜绝不同回路间的电缆交叉。 - - **容量约束**：自动计算所需的最小回路数（Clusters），确保每条集电线路的总功率不超过海缆极限。 - -### 3. 精细化电气计算与选型 -- **动态电缆选型**： - - 基于实际潮流计算（Power Flow），为每一段线路选择最经济且满足载流量的电缆。 - - 规格库：覆盖 35mm² 至 400mm² 海缆。 - - 参数：电压等级 **66kV**，功率因数 0.95。 -- **成本与损耗评估**： - - 考虑海缆材料及敷设成本（约为陆缆的5倍）。 - - 计算全场集电线路的 $I^2R$ 损耗。 - -### 4. 工程级可视化与输出 -- **可视化图表**： - - 生成直观的拓扑连接图。 - - **颜色编码**：使用不同颜色和粗细区分不同截面的电缆（如绿色细线为35mm²，红色粗线为400mm²）。 - - 自动保存为高清 PNG 图片。 -- **CAD (DXF) 导出**： - - 使用 `ezdxf` 生成 `.dxf` 文件。 - - 分层管理：风机、升压站、各规格电缆分层显示，可直接导入 AutoCAD 进行后续工程设计。 - -## 安装说明 - -### 环境要求 -- Python >= 3.10 -- 推荐使用 [uv](https://github.com/astral-sh/uv) 进行依赖管理。 - -### 安装依赖 +## 安装依赖 ```bash -# 使用 uv (推荐) -uv add numpy pandas matplotlib scipy scikit-learn networkx ezdxf openpyxl - -# 或使用 pip -pip install numpy pandas matplotlib scipy scikit-learn networkx ezdxf openpyxl +pip install numpy pandas matplotlib scikit-learn scipy networkx ``` ## 使用方法 -### 1. 运行主程序 +### 基本用法 ```bash -# 使用 uv -uv run main.py - -# 或直接运行 python main.py ``` -### 2. 数据输入模式 +### 指定数据文件 -程序会自动检测当前目录下是否存在 `coordinates.xlsx`： +```bash +python main.py --excel wind_farm_coordinates.xlsx +``` -- **存在**：优先读取 Excel 文件中的坐标数据进行计算。 -- **不存在**：自动生成 30 台风机的规则布局（Grid Layout）进行演示。 +### 覆盖默认簇数 -### 3. 结果输出 +```bash +python main.py --clusters 20 +``` -程序运行结束后会： -1. 在终端打印详细的成本、损耗及电缆统计数据。 -2. 弹窗显示拓扑对比图，并保存为 `wind_farm_design_imported.png` (或 `offshore_...png`)。 -3. 生成 CAD 图纸文件 `wind_farm_design.dxf`。 +## 算法说明 + +### 1. MST Method（最小生成树） +- 使用最小生成树连接所有风机到海上变电站 +- 简单高效，适合初步设计 + +### 2. K-means Clustering +- 将风机分组到多个回路中 +- 平衡每回路的功率分配 + +### 3. Capacitated Sweep（容量扫描） +- 考虑电缆容量约束的智能分组 +- 支持多种电缆规格自动选择 + +### 4. Rotational Sweep（旋转优化） +- 在容量扫描基础上进行旋转优化 +- 进一步降低总成本和损耗 + +## 输出文件 + +1. **可视化图片**：`wind_farm_design_comparison.png` + - 不同算法的设计方案对比图 + +2. **CAD图纸**：`wind_farm_design.dxf` + - 可导入CAD软件的详细设计图纸 + +3. **数据报告**：`wind_farm_design.xlsx` + - 包含所有方案的详细技术参数和成本分析 ## 关键参数说明 @@ -91,18 +79,36 @@ POWER_FACTOR = 0.95 # 功率因数 cost_multiplier = 5.0 # 海缆相对于陆缆的成本倍数 ``` +## 电缆规格配置 + +项目支持多种电缆规格，可在 `generate_template.py` 中配置： + +| 截面积(mm²) | 容量(MW) | 电阻(Ω/km) | 成本(元/m) | +|-------------|----------|------------|------------| +| 35 | 150 | 0.524 | 80 | +| 70 | 215 | 0.268 | 120 | +| 95 | 260 | 0.193 | 150 | +| 120 | 295 | 0.153 | 180 | +| 150 | 330 | 0.124 | 220 | +| 185 | 370 | 0.0991 | 270 | +| 240 | 425 | 0.0754 | 350 | +| 300 | 500 | 0.0601 | 450 | +| 400 | 580 | 0.0470 | 600 | + ## 输出示例 ```text -系统设计参数: 总功率 2000.0 MW, 单回路最大容量 50.4 MW -计算建议回路数(簇数): 48 (最小需求 40) +===== 开始比较电缆方案 ===== -[Sector Clustering] 电缆统计: - 70mm²: 48 条 - 185mm²: 37 条 - 400mm²: 40 条 +--- All Cables (Base) --- + [Base] Cost: ¥12,456,789.12 | Loss: 234.56 kW + [Rotational] Cost: ¥12,234,567.89 | Loss: 223.45 kW -成功导出DXF文件: wind_farm_design.dxf +--- High Current (Base) --- + [Base] Cost: ¥11,987,654.32 | Loss: 245.67 kW + [Rotational] Cost: ¥11,876,543.21 | Loss: 234.56 kW + +推荐方案: High Current (Rotational) (默认) ``` ## 许可证 diff --git a/esau_williams.py b/esau_williams.py new file mode 100644 index 0000000..b38aa4c --- /dev/null +++ b/esau_williams.py @@ -0,0 +1,241 @@ +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 diff --git a/main.py b/main.py index 484ded5..0dc06b9 100644 --- a/main.py +++ b/main.py @@ -9,6 +9,7 @@ import networkx as nx import math import argparse import os +from esau_williams import design_with_esau_williams # 设置matplotlib支持中文显示 plt.rcParams['font.sans-serif'] = ['Microsoft YaHei', 'SimHei', 'Arial'] @@ -1150,9 +1151,32 @@ def compare_design_methods(excel_path=None, n_clusters_override=None): }) print(f" [Rotational] Cost: ¥{eval_rot['total_cost']:,.2f} | Loss: {eval_rot['total_loss']:.2f} kW") + # --- Run 3: Esau-Williams Algorithm --- + ew_name = f"{name} (Esau-Williams)" + conns_ew, turbines_ew = design_with_esau_williams( + turbines.copy(), substation, max_cable_mw + ) + eval_ew = evaluate_design( + turbines, conns_ew, substation, cable_specs=current_specs, + is_offshore=is_offshore, method_name=ew_name + ) + + comparison_results.append({ + 'name': ew_name, + 'cost': eval_ew['total_cost'], + 'loss': eval_ew['total_loss'], + 'eval': eval_ew, + 'turbines': turbines_ew, + 'specs': current_specs + }) + print(f" [Esau-Williams] Cost: ¥{eval_ew['total_cost']:,.2f} | Loss: {eval_ew['total_loss']:.2f} kW") + # 记录最佳 if eval_rot['total_cost'] < best_cost: best_cost = eval_rot['total_cost'] + + if eval_ew['total_cost'] < best_cost: + best_cost = eval_ew['total_cost'] # best_result 不再需要单独维护，最后遍历 comparison_results 即可 if eval_base['total_cost'] < best_cost: