python+gurobi复现代码MIP求解infeasible

Question

复现文章：Layout Design of Parking Lots with Mathematical Programming
Konrad Stephan, Felix Weidinger, and Nils Boysen
Transportation Science 2021 55:4, 930-945

背景是停车场布局设计，我正在实现第一个最简单的MIP模型，符号、约束如下图

用computeIIS()查找矛盾约束，感觉问题出现在约束7-10上，又有大M约束，又有绝对值，但是笔算一通也没看出来：

我的约束7-10的代码：

for i in range(m1 + 1):
    for j in range(n1 + 1):
        model_r1.addConstr(fr[i,j] <= M * y[i,j])
        model_r1.addConstr(fr[i,j] >= -M * y[i,j])
        model_r1.addConstr(fr[i,j] <= M * y[i,j+1])
        model_r1.addConstr(fr[i,j] >= -M * y[i,j+1])
        model_r1.addConstr(fd[i,j] <= M * y[i,j])
        model_r1.addConstr(fd[i,j] >= -M * y[i,j])
        model_r1.addConstr(fd[i,j] <= M * y[i+1,j])
        model_r1.addConstr(fd[i,j] >= -M * y[i+1,j])

我的尝试：

1. 如上所述，参考csdn里的文章，用computeIIS()查找矛盾约束，没找出来
2. 又参考另一篇文章，判断是否为大M取值的问题，从1e1试到1e10发现都不行
3. 把约束注释掉｜fr[i,j]｜<= M * y[i,j+1]和｜fd[i,j]｜<= M * y[i,j+1]中的任意两条后发现模型可行

求指点！！

Answer 1

引自免费微信小程序：皆我百晓生

根据您提供的符号定义和约束条件，以及所描述的问题，模型中出现的不可行性可能源于流量变量（fr[i,j] 和 fd[i,j]）与街道区域标记变量（y[i,j]）之间的大M约束。当这些约束导致了不等式系统的不可行时，您可以尝试以下方法来调试：

1. 检查大M值选取：虽然您已经尝试了不同的大M值范围，但有时问题可能在于特定约束组合下的M值太小或太大。对于此类问题，一个常见的调试步骤是仔细分析每个约束，并尝试找出理论上合理的最小和最大流量值，然后选择适当的大M值以确保松弛足够。

2. 简化约束：

   • 由于您提到注释掉部分约束后模型变为可行，可以逐个分析这些约束是否在某些特殊情况下产生矛盾。
   • 对于约束 fr[i,j] <= M * y[i,j+1] 和 fr[i,j] >= -M * y[i,j+1]，考虑是否可以在不损失模型表达力的前提下将其合并为单个绝对值约束，例如：|fr[i,j]| <= M * y[i,j+1]。Gurobi 支持使用特殊的绝对值约束构造器直接添加这种形式的约束。

3. 修正约束方向：

   • 检查流量变量的方向性是否与实际停车布局中的车辆行驶方向一致。如果模型中存在逻辑错误，可能导致流向约束的方向设置有误。

4. 分块调试：

   • 尝试仅针对某一行或某一列应用网络流约束，观察是否能够找到某个局部子集引起的不可行性。

5. IIS（不可行解解释集）的精确解读：

   • 虽然您用 computeIIS() 方法没有直接找到矛盾约束，但有时候可能需要更深入地理解输出结果。确保正确识别并解析IIS报告中给出的约束集合，它们可能隐含着模型结构上的问题。

6. 代码实现修正： 根据您给出的约束代码片段，调整如下，将每对流量变量与相邻街道变量之间的约束整合在一个循环内，并采用绝对值约束：

import gurobipy as gp

# 假设已有的模型、变量及参数初始化
model_r1 = gp.Model()
M = model_r1.addVar(name="BigM", lb=0, ub=gp.GRB.INFINITY) # 大M变量
for i in range(m1 + 1):
    for j in range(n1 + 1):
        # 其他变量定义略...
        # 添加绝对值约束
        model_r1.addConstr(gp.abs_(fr[i,j]) <= M * y[i,j])
        model_r1.addConstr(gp.abs_(fr[i,j]) <= M * y[i,j+1]) if j < n1 else None
        model_r1.addConstr(gp.abs_(fd[i,j]) <= M * y[i,j])
        model_r1.addConstr(gp.abs_(fd[i,j]) <= M * y[i+1,j]) if i < m1 else None

请注意，在实际实现过程中，Gurobi 直接支持绝对值约束，因此应使用 gp.abs_() 函数而不是分开处理大于等于和小于等于约束。同时，请确认边界条件（如 j < n1 和 i < m1）以避免数组越界。

通过上述建议进行细致排查和修正，应该有助于定位并解决模型不可行的问题。在实施这些方案时，请密切监控模型状态和求解日志，以获得更多的线索和调试信息。