Fix MIP algorithm: simplify model formulation and add detailed debugging

mip.py

@@ -19,7 +19,7 @@ def design_with_mip(
     power_factor=0.95,
     system_params=None,
     max_clusters=None,
-    time_limit=300,  # seconds
+    time_limit=300,
     evaluate_func=None,
     total_invest_func=None,
     get_max_capacity_func=None,
@@ -49,45 +49,37 @@ def design_with_mip(
     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()
+        max_mw = 100.0
     if max_clusters is None:
-        max_clusters = int(np.ceil(total_power / max_mw))
+        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)
 
-    # MIP 模型
     prob = pulp.LpProblem("WindFarmCollectorMIP", pulp.LpMinimize)
 
-    # 创建所有变量和约束的辅助函数
     def assign_var(i, k):
         return pulp.LpVariable(f"assign_{i}_{k}", cat="Binary")
 
     def cluster_var(k):
         return pulp.LpVariable(f"cluster_{k}", cat="Binary")
 
-    # 目标函数：最小化使用的簇数（越少回路，成本越低）
     prob += pulp.lpSum([cluster_var(k) for k in range(max_clusters)])
 
-    # 约束1：每个风机分配到一个簇
     for i in range(n_turbines):
         prob += pulp.lpSum([assign_var(i, k) for k in range(max_clusters)]) == 1
 
-    # 约束2：簇功率约束（允许20%过载，增加求解器灵活性）
     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)
 
-    # 约束3：如果簇未使用，则无分配
     for k in range(max_clusters):
         for i in range(n_turbines):
             prob += assign_var(i, k) <= cluster_var(k)
@@ -96,9 +88,8 @@ def design_with_mip(
         f"MIP Model: {len(prob.variables())} variables, {len(prob.constraints)} constraints"
     )
 
-    # 求解
-    solver = pulp.PULP_CBC_CMD(timeLimit=time_limit, msg=0, warmStart=False)
     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}"
@@ -108,25 +99,36 @@ def design_with_mip(
         print(
             f"MIP solver status: {pulp.LpStatus[prob.status]}, solution not found, falling back to MST"
         )
-        print(f"Model infeasibility check:")
-        print(
-            f"Total power: {turbines['power'].sum():.2f} MW, Max cluster capacity: {max_mw:.2f} MW"
-        )
+        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:
@@ -134,12 +136,10 @@ def design_with_mip(
                 cluster_assign[i] = k
                 assigned = True
                 break
-        if not assigned and active_clusters:
-            # 如果没有分配，分配到距离最近的活跃簇
+        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)
@@ -162,7 +162,6 @@ def design_with_mip(
                                 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)))