PuLP中约束取反后解不一致的技术咨询

我使用PuLP编写了以最小化无约束gamma变量为目标的线性规划模型，原模型包含≤和≥混合约束时求解显示不可行；但将目标1、2、4、5的≥约束取反转换为≤约束后，却得到了可行解。从数学角度看这两种约束形式等价，想了解PuLP运行机制中导致该差异的原因。

原约束代码 #

import pulp as lp

# Define the problem
prob = lp.LpProblem("Annual_usage", lp.LpMinimize)

# Decision variables
# x_i is the project type
x = {
    'X1': lp.LpVariable('x_1', lowBound=0, cat='Binary'),
    'X2': lp.LpVariable('x_2', lowBound=0, cat='Binary'),
    'X3': lp.LpVariable('x_3', lowBound=0, cat='Binary'),
    'X4': lp.LpVariable('x_4', lowBound=0, cat='Binary'),
    'X5': lp.LpVariable('x_5', lowBound=0, cat='Binary'),
    'X6': lp.LpVariable('x_6', lowBound=0, cat='Binary'),
    'X7': lp.LpVariable('x_7', lowBound=0, cat='Binary'),
    'X8': lp.LpVariable('x_8', lowBound=0, cat='Binary'),
    'gamma': lp.LpVariable('gamma', lowBound=None, cat='Continuous')
}

# Usage
usage = {
    'X1': 4.7, 'X2': 12.5, 'X3': 3.2, 'X4': 7.5, 'X5': 41,  'X6': 47,   'X7': 23,  'X8': 16,
    'gamma': 0
}

# Costs
cost = {
    'X1': 75,  'X2': 180,  'X3': 350, 'X4': 45,  'X5': 120, 'X6': 80,   'X7': 115, 'X8': 210,
    'gamma': 0
}

# Acreage
acreage = {
    'X1': 7,   'X2': 12,   'X3': 20,  'X4': 6, 'X5': 3,   'X6': 25,   'X7': 5,   'X8': 8,
    'gamma': 0
}

goal1 = { # at least 1 of X3 and X6
    'X1': 0,   'X2': 0,   'X3': 1,  'X4': 0, 'X5': 0,   'X6':  1,   'X7': 0,   'X8': 0,
    'gamma': 16/6
}

goal2 = { # at least 130
    'X1': 4.7, 'X2': 12.5, 'X3': 3.2, 'X4': 7.5, 'X5': 41,  'X6': 47,   'X7': 23,  'X8': 16,
    'gamma': 16/3
}

goal3 = { #sum is at most 7
    'X1': 3,   'X2': 2,    'X3': 1,   'X4': 3, 'X5': 2,   'X6':  1,    'X7': 2,   'X8': 3,
    'gamma': 16/2
}

goal4 = { #at least 495)
    'X1': 75,  'X2': 180,  'X3': 350, 'X4': 45,  'X5': 120, 'X6': 80,   'X7': 115, 'X8': 210,
    'gamma': 16
}

goal5 = { # sum is at least 5
    'X1': 1,   'X2': 1,   'X3': 1,  'X4': 1, 'X5': 1,   'X6':  1,   'X7': 1,   'X8': 1,
    'gamma': 4
}

gamma_objective = { 
    'X1': 0,   'X2': 0,   'X3': 0,  'X4': 0, 'X5': 0,   'X6':  0,   'X7': 0,   'X8': 0,
    'gamma': 1
}


# Objective function
prob += lp.lpSum(gamma_objective[i] * x[i] for i in x)

# Constraints
prob += lp.lpSum(x[i] * goal1[i] for i in x) >= 1, "goal 1"
prob += lp.lpSum(x[i] * goal2[i] for i in x) >= 130, "goal 2"
prob += lp.lpSum(x[i] * goal3[i] for i in x) <= 7,  "goal 3"
prob += lp.lpSum(x[i] * goal4[i] for i in x) >= 495, "goal 4"
prob += lp.lpSum(x[i] * goal5[i] for i in x) >= 5, "goal 5"

prob += lp.lpSum(x[i] * cost[i] for i in x) <=550, "# cost constraint"
prob += lp.lpSum(x[i] * acreage[i] for i in x) <=50, "# acreage constraint"


# Solve the problem
prob.solve()

# Print the solution
for v in prob.variables():
    print(v.name, "=", v.varValue)

print("Total Cost =", lp.value(prob.objective))

status = prob.status
print("Status:", lp.LpStatus[status])

约束取反后的代码 #

import pulp as lp

# Define the problem
prob = lp.LpProblem("Annual_usage", lp.LpMinimize)

# Decision variables
# x_i is the project type
x = {
    'X1': lp.LpVariable('x_1', lowBound=0, cat='Binary'),
    'X2': lp.LpVariable('x_2', lowBound=0, cat='Binary'),
    'X3': lp.LpVariable('x_3', lowBound=0, cat='Binary'),
    'X4': lp.LpVariable('x_4', lowBound=0, cat='Binary'),
    'X5': lp.LpVariable('x_5', lowBound=0, cat='Binary'),
    'X6': lp.LpVariable('x_6', lowBound=0, cat='Binary'),
    'X7': lp.LpVariable('x_7', lowBound=0, cat='Binary'),
    'X8': lp.LpVariable('x_8', lowBound=0, cat='Binary'),
    'gamma': lp.LpVariable('gamma', lowBound=None, cat='Continuous')
}

# Usage
usage = {
    'X1': 4.7, 'X2': 12.5, 'X3': 3.2, 'X4': 7.5, 'X5': 41,  'X6': 47,   'X7': 23,  'X8': 16,
    'gamma': 0
}

# Costs
cost = {
    'X1': 75,  'X2': 180,  'X3': 350, 'X4': 45,  'X5': 120, 'X6': 80,   'X7': 115, 'X8': 210,
    'gamma': 0
}

# Acreage
acreage = {
    'X1': 7,   'X2': 12,   'X3': 20,  'X4': 6, 'X5': 3,   'X6': 25,   'X7': 5,   'X8': 8,
    'gamma': 0
}

goal1 = { # at least X3 and X6
    'X1': 0,   'X2': 0,   'X3': -1,  'X4': 0, 'X5': 0,   'X6':  -1,   'X7': 0,   'X8': 0,
    'gamma': -16/6
}

goal2 = { # at least 130
    'X1': -4.7, 'X2': -12.5, 'X3': -3.2, 'X4': -7.5, 'X5': -41,  'X6': -47,   'X7': -23,  'X8': -16,
    'gamma': -16/3
}

goal3 = { #sum is at most 7
    'X1': -3,   'X2': -2,    'X3': -1,   'X4': -3, 'X5': -2,   'X6':  -1,    'X7': -2,   'X8':-3,
    'gamma': -16/2
}

goal4 = { # At  least 495k
    'X1': -75,  'X2': -180,  'X3': -350, 'X4': -45,  'X5': -120, 'X6': -80,   'X7': -115, 'X8': -210,
    'gamma': -16
}

goal5 = { # sum is at least 5
    'X1': -1,   'X2': -1,   'X3': -1,  'X4': -1, 'X5': -1,   'X6':  -1,   'X7': -1,   'X8': -1,
    'gamma': -4
}

gamma_objective = { 
    'X1': 0,   'X2': 0,   'X3': 0,  'X4': 0, 'X5': 0,   'X6':  0,   'X7': 0,   'X8': 0,
    'gamma': 1
}


# Objective function
prob += lp.lpSum(gamma_objective[i] * x[i] for i in x)

# Constraints
prob += lp.lpSum(x[i] * goal1[i] for i in x) <= -1, "goal 1"
prob += lp.lpSum(x[i] * goal2[i] for i in x) <= -130, "goal 2"
prob += lp.lpSum(x[i] * goal3[i] for i in x) <= 7,  "goal 3"
prob += lp.lpSum(x[i] * goal4[i] for i in x) <= -495, "goal 4"
prob += lp.lpSum(x[i] * goal5[i] for i in x) <= -5, "goal 5"

prob += lp.lpSum(x[i] * cost[i] for i in x) <=550, "# cost constraint"
prob += lp.lpSum(x[i] * acreage[i] for i in x) <=50, "# acreage constraint"


# Solve the problem
prob.solve()

# Print the solution
for v in prob.variables():
    print(v.name, "=", v.varValue)

print("Total Cost =", lp.value(prob.objective))

status = prob.status
print("Status:", lp.LpStatus[status])

1. 求解器数值精度与约束表述的差异 #

PuLP仅作为建模接口，实际求解依赖底层求解器（如CBC、GLPK）。虽然两种约束形式数学等价，但求解器的可行性检查逻辑会因表述不同而产生数值误差差异：

• 对于A >= b形式的约束，求解器会验证A - b >= -ε（ε为求解器内置容差）；而转换为-A <= -b后，验证逻辑变为-A + b <= ε。模型中大量浮点数系数（如16/6、4.7）的舍入误差，在两种验证逻辑下的累积方向不同，可能导致一种情况被判定为不可行，另一种通过容差检查。

2. 无界gamma变量的迭代影响 #

你的gamma变量未设置上下界（lowBound=None），在最小化目标驱动下，求解器会尝试将gamma推向负无穷。当约束为>=形式时，gamma的极端取值可能导致约束计算出现数值溢出或精度丢失，进而让求解器误判模型不可行；而转换为<=约束后，约束的数值计算方向更贴合求解器的迭代逻辑，避免了这类极端值引发的问题。

3. 约束标准化的细微差异 #

PuLP会将模型转换为求解器要求的标准形式（通常统一为<=约束，无界变量拆分为两个非负变量的差）。手动转换约束与让PuLP自动转换，在标准化过程中可能产生细微的系数舍入差异，这也会影响求解器的可行性判断。

• 显式指定求解器并调整容差，例如使用prob.solve(lp.PULP_CBC_CMD(gapRel=1e-6, maxSeconds=10))，观察原模型求解结果是否变化。
• 查看可行解中gamma的取值，确认是否是极端值导致原模型约束计算出现数值问题。
• 将所有浮点数系数替换为分数形式（如用Fraction(16,6)替代16/6），减少浮点数精度误差。

内容的提问来源于stack exchange，提问作者gunsnfloyd