Rewrite MIP model formulation and add comprehensive debugging

mip.py

@@ -51,9 +51,14 @@ def design_with_mip(
     else:
         max_mw = 100.0  # 默认值
     total_power = turbines["power"].sum()
-    max_clusters = int(np.ceil(total_power / max_mw))
+    if max_clusters is None:
+        max_clusters = int(np.ceil(total_power / 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)
@@ -61,59 +66,78 @@ def design_with_mip(
     # MIP 模型
     prob = pulp.LpProblem("WindFarmCollectorMIP", pulp.LpMinimize)
 
-    # 决策变量：风机分配到簇 (binary)
-    x = pulp.LpVariable.dicts(
-        "assign", (range(n_turbines), range(max_clusters)), cat="Binary"
-    )
+    # 创建所有变量和约束的辅助函数
+    def assign_var(i, k):
+        return pulp.LpVariable(f"assign_{i}_{k}", cat="Binary")
 
-    # 簇使用变量 (binary)
-    y = pulp.LpVariable.dicts("use_cluster", range(max_clusters), cat="Binary")
+    def cluster_var(k):
+        return pulp.LpVariable(f"cluster_{k}", cat="Binary")
 
-    # 目标函数：最小化总成本 (线性简化版：到升压站距离总和)
-    # 由于MIP线性约束，简化目标为最小化风机到升压站的距离总和（通过簇）
-    prob += pulp.lpSum(
-        [
-            dist_matrix_full[0, i + 1] * x[i][k]
-            for i in range(n_turbines)
-            for k in range(max_clusters)
-        ]
-    )
+    # 目标函数：最小化使用的簇数（越少回路，成本越低）
+    prob += pulp.lpSum([cluster_var(k) for k in range(max_clusters)])
 
-    # 约束：每个风机分配到一个簇
+    # 约束1：每个风机分配到一个簇
     for i in range(n_turbines):
-        prob += pulp.lpSum([x[i][k] for k in range(max_clusters)]) == 1
+        prob += pulp.lpSum([assign_var(i, k) for k in range(max_clusters)]) == 1
 
-    # 簇功率约束
+    # 约束2：簇功率约束（允许20%过载，增加求解器灵活性）
     for k in range(max_clusters):
-        prob += (
-            pulp.lpSum([turbines.iloc[i]["power"] * x[i][k] for i in range(n_turbines)])
-            <= max_mw * y[k]
+        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 += x[i][k] <= y[k]
+            prob += assign_var(i, k) <= cluster_var(k)
+
+    print(
+        f"MIP Model: {len(prob.variables())} variables, {len(prob.constraints)} constraints"
+    )
 
     # 求解
-    solver = pulp.PULP_CBC_CMD(timeLimit=time_limit)
+    solver = pulp.PULP_CBC_CMD(timeLimit=time_limit, msg=0, warmStart=False)
+    print("MIP: Starting to solve...")
     status = prob.solve(solver)
+    print(
+        f"MIP: Solver status={pulp.LpStatus[prob.status]}, Objective value={pulp.value(prob.objective):.4f}"
+    )
 
     if pulp.LpStatus[prob.status] != "Optimal":
-        print(f"MIP not optimal: {pulp.LpStatus[prob.status]}")
-        # 返回默认方案，如 MST
+        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(f"Number of clusters: {max_clusters}, Number of turbines: {n_turbines}")
+
         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):
-        for k in range(max_clusters):
-            if pulp.value(x[i][k]) > 0.5:
+        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 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)
@@ -138,10 +162,13 @@ 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)))
 
     turbines["cluster"] = cluster_assign
+    print(
+        f"MIP optimization completed successfully, {len(connections)} connections generated"
+    )
     return connections, turbines