Change MIP objective function to minimize total investment

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(