441 lines
14 KiB
Python
441 lines
14 KiB
Python
import numpy as np
|
||
import pandas as pd
|
||
from scipy.spatial import distance_matrix
|
||
from scipy.sparse.csgraph import minimum_spanning_tree
|
||
from collections import defaultdict
|
||
import random
|
||
|
||
try:
|
||
import pulp
|
||
|
||
pulp_available = True
|
||
except ImportError:
|
||
pulp = None
|
||
pulp_available = False
|
||
|
||
try:
|
||
import pyomo.environ as pyo_env
|
||
|
||
pyomo_available = True
|
||
except (ImportError, AttributeError):
|
||
pyomo_available = False
|
||
print("Pyomo not available, falling back to PuLP")
|
||
|
||
|
||
def design_with_pyomo(
|
||
turbines,
|
||
substation,
|
||
cable_specs=None,
|
||
voltage=66000,
|
||
power_factor=0.95,
|
||
system_params=None,
|
||
max_clusters=None,
|
||
time_limit=300,
|
||
evaluate_func=None,
|
||
total_invest_func=None,
|
||
get_max_capacity_func=None,
|
||
):
|
||
"""
|
||
使用Pyomo求解器优化集电线路布局
|
||
:param turbines: 风机DataFrame
|
||
:param substation: 升压站坐标
|
||
:param cable_specs: 电缆规格
|
||
:param system_params: 系统参数（用于NPV计算）
|
||
:param max_clusters: 最大簇数，默认基于功率计算
|
||
:param time_limit: 求解时间限制（秒）
|
||
:param evaluate_func: 评估函数
|
||
:param total_invest_func: 总投资计算函数
|
||
:param get_max_capacity_func: 获取最大容量函数
|
||
:return: 连接列表和带有簇信息的turbines
|
||
"""
|
||
if get_max_capacity_func:
|
||
max_mw = get_max_capacity_func(cable_specs, voltage, power_factor)
|
||
else:
|
||
max_mw = 100.0
|
||
total_power = turbines["power"].sum()
|
||
if max_clusters is None:
|
||
max_clusters = int(np.ceil(total_power / max_mw))
|
||
n_turbines = len(turbines)
|
||
|
||
all_coords = np.vstack([substation, turbines[["x", "y"]].values])
|
||
dist_matrix_full = distance_matrix(all_coords, all_coords)
|
||
|
||
# Simple fallback for now - use PuLP instead
|
||
print("Pyomo not fully implemented, falling back to PuLP")
|
||
return design_with_mip(
|
||
turbines,
|
||
substation,
|
||
cable_specs,
|
||
voltage,
|
||
power_factor,
|
||
system_params,
|
||
max_clusters,
|
||
time_limit,
|
||
evaluate_func,
|
||
total_invest_func,
|
||
get_max_capacity_func,
|
||
)
|
||
|
||
|
||
def design_with_mip(
|
||
turbines,
|
||
substation,
|
||
cable_specs=None,
|
||
voltage=66000,
|
||
power_factor=0.95,
|
||
system_params=None,
|
||
max_clusters=None,
|
||
time_limit=300,
|
||
evaluate_func=None,
|
||
total_invest_func=None,
|
||
get_max_capacity_func=None,
|
||
):
|
||
"""
|
||
使用混合整数规划(MIP)优化集电线路布局
|
||
:param turbines: 风机DataFrame
|
||
:param substation: 升压站坐标
|
||
:param cable_specs: 电缆规格
|
||
:param system_params: 系统参数（用于NPV计算）
|
||
:param max_clusters: 最大簇数，默认基于功率计算
|
||
:param time_limit: 求解时间限制（秒）
|
||
:param evaluate_func: 评估函数
|
||
:param total_invest_func: 总投资计算函数
|
||
:param get_max_capacity_func: 获取最大容量函数
|
||
:return: 连接列表和带有簇信息的turbines
|
||
"""
|
||
if not pulp_available:
|
||
print(
|
||
"WARNING: PuLP library not available. MIP optimization skipped, falling back to MST."
|
||
)
|
||
from main import design_with_mst
|
||
|
||
connections = design_with_mst(turbines, substation)
|
||
return connections, turbines
|
||
|
||
if get_max_capacity_func:
|
||
max_mw = get_max_capacity_func(cable_specs, voltage, power_factor)
|
||
else:
|
||
max_mw = 100.0
|
||
if max_clusters is None:
|
||
max_clusters = int(np.ceil(turbines["power"].sum() / max_mw))
|
||
n_turbines = len(turbines)
|
||
|
||
print(
|
||
f"MIP Model Setup: n_turbines={n_turbines}, max_clusters={max_clusters}, max_mw={max_mw:.2f} MW"
|
||
)
|
||
|
||
all_coords = np.vstack([substation, turbines[["x", "y"]].values])
|
||
dist_matrix_full = distance_matrix(all_coords, all_coords)
|
||
|
||
prob = pulp.LpProblem("WindFarmCollectorMIP", pulp.LpMinimize)
|
||
|
||
# Create all decision variables upfront to avoid duplicates
|
||
assign_vars = {}
|
||
for i in range(n_turbines):
|
||
for k in range(max_clusters):
|
||
assign_vars[(i, k)] = pulp.LpVariable(f"assign_{i}_{k}", cat="Binary")
|
||
|
||
cluster_vars = {}
|
||
for k in range(max_clusters):
|
||
cluster_vars[k] = pulp.LpVariable(f"cluster_{k}", cat="Binary")
|
||
|
||
# Helper functions to access variables
|
||
def assign_var(i, k):
|
||
return assign_vars[(i, k)]
|
||
|
||
def cluster_var(k):
|
||
return cluster_vars[k]
|
||
|
||
# Simplified objective function: minimize total distance
|
||
prob += pulp.lpSum(
|
||
[
|
||
dist_matrix_full[0, i + 1] * assign_var(i, k)
|
||
for i in range(n_turbines)
|
||
for k in range(max_clusters)
|
||
]
|
||
)
|
||
|
||
for i in range(n_turbines):
|
||
prob += pulp.lpSum([assign_var(i, k) for k in range(max_clusters)]) == 1
|
||
|
||
for k in range(max_clusters):
|
||
cluster_power = pulp.lpSum(
|
||
[turbines.iloc[i]["power"] * assign_var(i, k) for i in range(n_turbines)]
|
||
)
|
||
prob += cluster_power <= max_mw * 1.2 * cluster_var(k)
|
||
|
||
for k in range(max_clusters):
|
||
for i in range(n_turbines):
|
||
prob += assign_var(i, k) <= cluster_var(k)
|
||
|
||
print(
|
||
f"MIP Model: {len(prob.variables())} variables, {len(prob.constraints)} constraints"
|
||
)
|
||
|
||
# Debug: Print model structure
|
||
print("MIP model structure check:")
|
||
print(f" Variables: {len(prob.variables())}")
|
||
print(f" Constraints: {len(prob.constraints)}")
|
||
print(f" Time limit: {time_limit}s")
|
||
print(f" Turbines: {n_turbines}, Clusters: {max_clusters}")
|
||
|
||
# Test solver availability
|
||
try:
|
||
import subprocess
|
||
|
||
test_solver = subprocess.run(
|
||
[
|
||
r"D:\code\windfarm\.venv\Lib\site-packages\pulp\apis\..\solverdir\cbc\win\i64\cbc.exe",
|
||
"-version",
|
||
],
|
||
capture_output=True,
|
||
text=True,
|
||
timeout=5,
|
||
)
|
||
print(
|
||
f"CBC solver test: {test_solver.stdout[:100] if test_solver.stdout else 'No output'}"
|
||
)
|
||
except Exception as solver_test_error:
|
||
print(f"CBC solver test failed: {solver_test_error}")
|
||
|
||
print("MIP: Starting to solve...")
|
||
try:
|
||
# Try to use CBC solver with different configurations
|
||
solver = pulp.PULP_CBC_CMD(
|
||
timeLimit=time_limit,
|
||
msg=False,
|
||
warmStart=False,
|
||
)
|
||
print(f"Using CBC solver with time limit: {time_limit}s")
|
||
status = prob.solve(solver)
|
||
print(
|
||
f"MIP: Solver status={pulp.LpStatus[prob.status]}, Objective value={pulp.value(prob.objective):.4f}"
|
||
)
|
||
except Exception as e:
|
||
print(f"MIP: CBC solver execution failed: {e}")
|
||
# Try alternative solver configurations
|
||
try:
|
||
print("MIP: Trying alternative solver configuration...")
|
||
solver = pulp.PULP_CBC_CMD(
|
||
msg=True, # Enable messages for debugging
|
||
threads=1, # Single thread
|
||
timeLimit=time_limit,
|
||
)
|
||
status = prob.solve(solver)
|
||
print(
|
||
f"MIP: Alternative solver status={pulp.LpStatus[prob.status]}, Objective value={pulp.value(prob.objective):.4f}"
|
||
)
|
||
except Exception as e2:
|
||
print(f"MIP: All solver attempts failed: {e2}, falling back to MST")
|
||
from main import design_with_mst
|
||
|
||
connections = design_with_mst(turbines, substation)
|
||
return connections, turbines
|
||
|
||
if pulp.LpStatus[prob.status] != "Optimal":
|
||
print(
|
||
f"MIP solver status: {pulp.LpStatus[prob.status]}, solution not found, falling back to MST"
|
||
)
|
||
print("Model feasibility check:")
|
||
print(f"Total power: {turbines['power'].sum():.2f} MW")
|
||
print(f"Max cluster capacity: {max_mw:.2f} MW")
|
||
print(f"Number of clusters: {max_clusters}, Number of turbines: {n_turbines}")
|
||
|
||
for k in range(max_clusters):
|
||
cluster_power = pulp.value(
|
||
pulp.lpSum(
|
||
[
|
||
turbines.iloc[i]["power"] * assign_var(i, k)
|
||
for i in range(n_turbines)
|
||
]
|
||
)
|
||
)
|
||
cluster_used = pulp.value(cluster_var(k))
|
||
print(
|
||
f"Cluster {k}: Power={cluster_power:.2f} MW (max {max_mw * 1.2:.2f}), Used={cluster_used}"
|
||
)
|
||
|
||
from main import design_with_mst
|
||
|
||
connections = design_with_mst(turbines, substation)
|
||
return connections, turbines
|
||
|
||
cluster_assign = [-1] * n_turbines
|
||
active_clusters = []
|
||
for k in range(max_clusters):
|
||
if pulp.value(cluster_var(k)) > 0.5:
|
||
active_clusters.append(k)
|
||
|
||
for i in range(n_turbines):
|
||
assigned = False
|
||
for k in active_clusters:
|
||
if pulp.value(assign_var(i, k)) > 0.5:
|
||
cluster_assign[i] = k
|
||
assigned = True
|
||
break
|
||
if not assigned:
|
||
dists = [dist_matrix_full[0, i + 1] for k in active_clusters]
|
||
cluster_assign[i] = active_clusters[np.argmin(dists)]
|
||
|
||
clusters = defaultdict(list)
|
||
for i, c in enumerate(cluster_assign):
|
||
clusters[c].append(i)
|
||
|
||
connections = []
|
||
for c, members in clusters.items():
|
||
if len(members) == 0:
|
||
continue
|
||
coords = turbines.iloc[members][["x", "y"]].values
|
||
if len(members) > 1:
|
||
dm = distance_matrix(coords, coords)
|
||
mst = minimum_spanning_tree(dm).toarray()
|
||
for i in range(len(members)):
|
||
for j in range(len(members)):
|
||
if mst[i, j] > 0:
|
||
connections.append(
|
||
(
|
||
f"turbine_{members[i]}",
|
||
f"turbine_{members[j]}",
|
||
mst[i, j],
|
||
)
|
||
)
|
||
dists = [dist_matrix_full[0, m + 1] for m in members]
|
||
closest = members[np.argmin(dists)]
|
||
connections.append((f"turbine_{closest}", "substation", min(dists)))
|
||
|
||
turbines["cluster"] = cluster_assign
|
||
|
||
# Check cluster distances
|
||
min_cluster_distance = check_cluster_distances(clusters, turbines)
|
||
if min_cluster_distance is not None:
|
||
print(
|
||
f"Cluster validation: Minimum distance between clusters = {min_cluster_distance:.2f} m"
|
||
)
|
||
if min_cluster_distance < 1000:
|
||
print(
|
||
f"WARNING: Clusters are very close to each other ({min_cluster_distance:.2f} m < 1000 m)"
|
||
)
|
||
elif min_cluster_distance < 2000:
|
||
print(
|
||
f"NOTICE: Clusters are relatively close ({min_cluster_distance:.2f} m)"
|
||
)
|
||
|
||
# Check for cable crossings
|
||
cable_crossings = check_cable_crossings(connections, turbines, substation)
|
||
if cable_crossings:
|
||
print(
|
||
f"WARNING: Found {len(cable_crossings)} cable crossing(s) in the solution"
|
||
)
|
||
for i, (idx1, idx2, p1, p2, p3, p4) in enumerate(cable_crossings):
|
||
conn1 = connections[idx1]
|
||
conn2 = connections[idx2]
|
||
print(
|
||
f" Crossing {i + 1}: Connection {conn1[0]}-{conn1[1]} crosses {conn2[0]}-{conn2[1]}"
|
||
)
|
||
else:
|
||
print("No cable crossings detected in the solution")
|
||
|
||
print(
|
||
f"MIP optimization completed successfully, {len(connections)} connections generated"
|
||
)
|
||
return connections, turbines
|
||
|
||
|
||
def calculate_cluster_centroids(clusters, turbines):
|
||
"""Calculate the centroid coordinates for each cluster."""
|
||
centroids = {}
|
||
for c, members in clusters.items():
|
||
if len(members) == 0:
|
||
centroids[c] = (0, 0)
|
||
else:
|
||
coords = turbines.iloc[members][["x", "y"]].values
|
||
centroid_x = np.mean(coords[:, 0])
|
||
centroid_y = np.mean(coords[:, 1])
|
||
centroids[c] = (centroid_x, centroid_y)
|
||
return centroids
|
||
|
||
|
||
def check_cluster_distances(clusters, turbines, min_distance_threshold=1000):
|
||
"""Check if any clusters are too close to each other."""
|
||
if len(clusters) < 2:
|
||
return None
|
||
|
||
centroids = calculate_cluster_centroids(clusters, turbines)
|
||
active_clusters = [c for c, members in clusters.items() if len(members) > 0]
|
||
|
||
min_distance = float("inf")
|
||
min_pair = None
|
||
|
||
for i in range(len(active_clusters)):
|
||
for j in range(i + 1, len(active_clusters)):
|
||
c1, c2 = active_clusters[i], active_clusters[j]
|
||
centroid1 = np.array(centroids[c1])
|
||
centroid2 = np.array(centroids[c2])
|
||
distance = np.linalg.norm(centroid1 - centroid2)
|
||
|
||
if distance < min_distance:
|
||
min_distance = distance
|
||
min_pair = (c1, c2)
|
||
|
||
return min_distance
|
||
|
||
|
||
def check_cable_crossings(connections, turbines, substation):
|
||
"""Check if there are cable crossings in the solution."""
|
||
crossings = []
|
||
|
||
def line_intersection(p1, p2, p3, p4):
|
||
"""Check if line segments (p1,p2) and (p3,p4) intersect."""
|
||
x1, y1 = p1
|
||
x2, y2 = p2
|
||
x3, y3 = p3
|
||
x4, y4 = p4
|
||
|
||
denom = (y4 - y3) * (x2 - x1) - (x4 - x3) * (y2 - y1)
|
||
|
||
if abs(denom) < 1e-10:
|
||
return False
|
||
|
||
ua = ((x4 - x3) * (y1 - y3) - (y4 - y3) * (x1 - x3)) / denom
|
||
ub = ((x2 - x1) * (y1 - y3) - (y2 - y1) * (x1 - x3)) / denom
|
||
|
||
return 0 <= ua <= 1 and 0 <= ub <= 1
|
||
|
||
def get_turbine_coord(connection_part):
|
||
"""Get coordinates from connection part (turbine_# or substation)."""
|
||
if connection_part == "substation":
|
||
# Ensure substation is returned as a proper tuple for unpacking
|
||
if isinstance(substation, (list, np.ndarray)):
|
||
return (substation[0], substation[1])
|
||
else:
|
||
return (substation[0], substation[1])
|
||
else:
|
||
turbine_idx = int(connection_part.split("_")[1])
|
||
return (
|
||
turbines.iloc[turbine_idx]["x"],
|
||
turbines.iloc[turbine_idx]["y"],
|
||
)
|
||
|
||
for i in range(len(connections)):
|
||
for j in range(i + 1, len(connections)):
|
||
conn1 = connections[i]
|
||
conn2 = connections[j]
|
||
|
||
p1 = get_turbine_coord(conn1[0])
|
||
p2 = get_turbine_coord(conn1[1])
|
||
p3 = get_turbine_coord(conn2[0])
|
||
p4 = get_turbine_coord(conn2[1])
|
||
|
||
if (
|
||
np.array_equal(p1, p3)
|
||
or np.array_equal(p1, p4)
|
||
or np.array_equal(p2, p3)
|
||
or np.array_equal(p2, p4)
|
||
):
|
||
continue
|
||
|
||
if line_intersection(p1, p2, p3, p4):
|
||
crossings.append((i, j, p1, p2, p3, p4))
|
||
|
||
return crossings
|