From 41ac6f3963a639da8023d3eb1ee6658fe1ff3352 Mon Sep 17 00:00:00 2001 From: dmy Date: Thu, 8 Jan 2026 13:01:36 +0800 Subject: [PATCH] Change MIP objective function to minimize total investment --- mip.py | 185 +++++++++++++++++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 179 insertions(+), 6 deletions(-) diff --git a/mip.py b/mip.py index a359f24..4854e3c 100644 --- a/mip.py +++ b/mip.py @@ -7,8 +7,167 @@ 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, + ) + y = pyo_env.Var( + max_clusters_list, + within=pyo_env.BinarySet, + initialize=False, + name="use_cluster", + ) + + prob.constraints += pyo_env.SummizationRule( + pyo_env.Summ(x[i, k] for i in n_turbines_list for k in max_clusters_list) == 1, + name="assign_one_turbine", + ) + + cluster_powers = [ + pyo_env.Summ([turbines.iloc[i]["power"] * x[i, k] for i in n_turbines_list]) + for k in max_clusters_list + ] + + prob += pyo_env.SummizationRule( + pyo_env.Summ(y[k] * max_mw for k in max_clusters_list), name="cluster_capacity" + ) + + for k in max_clusters_list: + prob += pyo_env.SummizationRule( + pyo_env.Summ(x[i, k] * turbines.iloc[i]["power"] for i in n_turbines_list), + name=f"turbine_assign_{k}", + ) + + for k in max_clusters_list: + prob += y[k] <= pyo_env.summation(x[i, k] for i in n_turbines_list) + + prob += pyo_env.Minimize(pyo_env.Summ([y[k] for k in max_clusters_list])) + + print( + f"Pyomo Model: {len(prob.variables)} variables, {len(prob.constraints)} constraints" + ) + + solver = pyo_env.SolverFactory("cbc") + print("Pyomo: Starting to solve...") + result = solver.solve(prob, time_limit=time_limit) + print( + f"Pyomo: Solver status={result.solver.status}, Termination condition={result.solver.termination_condition}" + ) + + if result.solver.status == pyo_env.SolverStatus.ok: + cluster_assign = [-1] * n_turbines + active_clusters = [] + + for k in range(max_clusters): + if pyo_env.value(y[k]) > 0.5: + active_clusters.append(k) + + for i in n_turbines_list: + assigned = False + for k in active_clusters: + if pyo_env.value(x[i, k]) > 0.5: + cluster_assign[i] = k + assigned = True + break + if not assigned and active_clusters: + 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 + print( + f"Pyomo optimization completed successfully, {len(connections)} connections generated" + ) + return connections, turbines def design_with_mip( @@ -37,7 +196,7 @@ def design_with_mip( :param get_max_capacity_func: 获取最大容量函数 :return: 连接列表和带有簇信息的turbines """ - if pulp is None: + if not pulp_available: print( "WARNING: PuLP library not available. MIP optimization skipped, falling back to MST." ) @@ -69,7 +228,14 @@ def design_with_mip( def cluster_var(k): return pulp.LpVariable(f"cluster_{k}", cat="Binary") - prob += pulp.lpSum([cluster_var(k) for k in range(max_clusters)]) + # 目标函数：最小化总投资（近似为风机到升压站的总距离） + 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 @@ -90,10 +256,17 @@ def design_with_mip( print("MIP: Starting to solve...") solver = pulp.PULP_CBC_CMD(timeLimit=time_limit, msg=0, warmStart=False, path=None) - status = prob.solve(solver) - print( - f"MIP: Solver status={pulp.LpStatus[prob.status]}, Objective value={pulp.value(prob.objective):.4f}" - ) + try: + 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: Solver execution failed: {e}, 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(