使用Pulp处理整数线性规划中的大参数问题

Pulp官方文档指出：

检查数字精度。若使用极大数值（高精度），通常会给求解器带来问题。例如，请勿在问题中使用100000000000这样的参数。

我的整数线性规划（ILP）中存在这类超大参数，尝试过统一缩放所有参数，但仅约束右侧数值过大，导致部分参数变为极小小数（比如约束3*x + 4*y == 10000000000，x、y为整数变量），求可行的问题转换方法。

我的代码 #

f = open('data\input13_test.txt', 'r')
lines = f.readlines()
cost = 0
for i in range(0, len(lines), 4):
    a_coeff_eq_1 = int(lines[i].split(':')[1].split(',')[0].split('+')[1])
    a_coeff_eq_2 = int(lines[i].split(':')[1].split(',')[1].split('+')[1])
    b_coeff_eq_1 = int(lines[i+1].split(':')[1].split(',')[0].split('+')[1])
    b_coeff_eq_2 = int(lines[i+1].split(':')[1].split(',')[1].split('+')[1])
    c_eq_1 = int(lines[i+2].split(':')[1].split(',')[0].split('=')[1]) + 10000000000000
    c_eq_2 = int(lines[i+2].split(':')[1].split(',')[1].split('=')[1]) + 10000000000000
    
    a_coeff_eq_1
    a_coeff_eq_2
    b_coeff_eq_1 
    b_coeff_eq_2
    c_eq_1
    c_eq_2
    
    prob = LpProblem("MyProb", LpMinimize)
    a = LpVariable("a", 0, None, LpInteger)
    b = LpVariable("b", 0, None, LpInteger)
    prob += a_coeff_eq_1 * a + b_coeff_eq_1 * b == c_eq_1
    prob += a_coeff_eq_2 * a + b_coeff_eq_2 * b == c_eq_2
    prob += 3*a + b
    prob.solve(PULP_CBC_CMD(msg=1))
    if LpStatus[prob.status] == 'Optimal':
        for v in prob.variables():
            if v.name == 'a':
                cost += 3*v.varValue
            elif v.name == 'b':
                cost += v.varValue
            print(v.name, "=", v.varValue)
            
print(cost)

输入文件示例 #

Button A: X+51, Y+11
Button B: X+38, Y+78
Prize: X=16146, Y=3706

可行转换方法 #

1. 变量拆分替换法（最适配你的场景）

你的问题核心是约束右侧加了固定超大常数10^13，变量系数都是小整数。可以通过拆分变量消除超大常数：

• 设原约束为：a_coeff * a + b_coeff * b = C + K（K=10^13，C是原输入的小数值）
• 拆分变量：a = a0 + a_k，b = b0 + b_k，其中a_k和b_k是满足a_coeff * a_k + b_coeff * b_k = K的整数解
• 替换后约束变为：a_coeff * a0 + b_coeff * b0 = C，所有参数回到小数值范围
• 目标函数变为：3*(a0+a_k) + (b0+b_k) = (3a0 + b0) + (3a_k + b_k)，其中3a_k + b_k是常数，不影响优化，只需优化3a0 + b0

修改后的代码示例

import math
from pulp import LpProblem, LpVariable, LpInteger, LpMinimize, PULP_CBC_CMD, LpStatus

f = open('data/input13_test.txt', 'r')
lines = f.readlines()
cost = 0
K = 10**13  # 提取固定超大常数

def extended_gcd(a, b):
    if b == 0:
        return (a, 1, 0)
    else:
        g, x, y = extended_gcd(b, a % b)
        return (g, y, x - (a // b) * y)

for i in range(0, len(lines), 4):
    # 读取系数与原C值
    a_coeff_eq_1 = int(lines[i].split(':')[1].split(',')[0].split('+')[1])
    a_coeff_eq_2 = int(lines[i].split(':')[1].split(',')[1].split('+')[1])
    b_coeff_eq_1 = int(lines[i+1].split(':')[1].split(',')[0].split('+')[1])
    b_coeff_eq_2 = int(lines[i+1].split(':')[1].split(',')[1].split('+')[1])
    c_eq_1 = int(lines[i+2].split(':')[1].split(',')[0].split('=')[1])
    c_eq_2 = int(lines[i+2].split(':')[1].split(',')[1].split('=')[1])
    
    # 求解a_k和b_k：满足两个方程的整数解
    a_k, b_k = None, None
    det = a_coeff_eq_1 * b_coeff_eq_2 - a_coeff_eq_2 * b_coeff_eq_1
    
    if det != 0:
        # 克莱姆法则求解，需确保K能被行列式整除
        if K % det != 0:
            print("当前组无解，跳过")
            continue
        a_k = (b_coeff_eq_2 * K - b_coeff_eq_1 * K) // det
        b_k = (a_coeff_eq_1 * K - a_coeff_eq_2 * K) // det
        # 验证解的正确性
        if not (a_coeff_eq_1 * a_k + b_coeff_eq_1 * b_k == K and a_coeff_eq_2 * a_k + b_coeff_eq_2 * b_k == K):
            print("解验证失败，跳过")
            continue
    else:
        # 系数矩阵线性相关，用扩展欧几里得找第一个方程的解
        g = math.gcd(a_coeff_eq_1, b_coeff_eq_1)
        if K % g != 0:
            print("当前组无解，跳过")
            continue
        g, x0, y0 = extended_gcd(a_coeff_eq_1, b_coeff_eq_1)
        a_k = x0 * (K // g)
        b_k = y0 * (K // g)
        # 验证第二个方程是否满足
        if a_coeff_eq_2 * a_k + b_coeff_eq_2 * b_k != K:
            print("当前组无解，跳过")
            continue
    
    # 构建新ILP问题，变量为a0、b0
    prob = LpProblem("MyProb", LpMinimize)
    a0 = LpVariable("a0", 0, None, LpInteger)
    b0 = LpVariable("b0", 0, None, LpInteger)
    
    # 新约束：系数*a0 + 系数*b0 = C
    prob += a_coeff_eq_1 * a0 + b_coeff_eq_1 * b0 == c_eq_1
    prob += a_coeff_eq_2 * a0 + b_coeff_eq_2 * b0 == c_eq_2
    
    # 目标函数：仅优化变量部分，常数项最后累加
    prob += 3*a0 + b0
    prob.solve(PULP_CBC_CMD(msg=1))
    
    if LpStatus[prob.status] == 'Optimal':
        a0_val = a0.varValue
        b0_val = b0.varValue
        # 还原原变量值
        a_val = a0_val + a_k
        b_val = b0_val + b_k
        cost += 3*a_val + b_val
        print(f"a = {a_val}, b = {b_val}")

print(f"总成本：{cost}")

注意事项

• 拆分后需保证a0 = a - a_k ≥ 0、b0 = b - b_k ≥0，如果初始找到的a_k或b_k为负，可以通过添加齐次解的倍数调整解的大小，确保替换后的变量非负。
• 若原问题无解，拆分后的问题也会无解，需提前判断。

2. 公约数提取法（适配部分场景）

如果约束左右两边存在共同公约数d，可将整个约束除以d，转换为整数约束（因x、y是整数，两边必须能被d整除）。例如约束2x +4y=1e10，可除以2得到x+2y=5e9，缩小数值规模缓解精度问题。但你的场景中系数与超大常数公约数为1，此方法不适用。

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