Gekko求解MINLP模型时约束不满足的技术问题咨询

问题背景与场景 #

现有一硬币支付场景：列表A包含3枚硬币，每枚硬币对应3种货币（v_0、v_1、v_2）的价值，初始每种货币总价值均为100000。每次随机生成3种货币的支付额payment，需通过Gekko求解0-1整数规划（MINLP）问题，选择硬币组合，满足选中硬币的各货币总价值≥对应支付额，并最大化目标函数：
$$\frac{1}{1 + total_input_v0 - payment[0]} + \frac{1}{1 + total_input_v1 - payment[1]} + \frac{1}{1 + total_input_v2 - payment[2]}$$
支付完成后，移除选中硬币，将找零作为新硬币加入A，同时添加随机存款硬币。

代码实现 #

主程序 #

if __name__ == "__main__":
    np.random.seed(2)

    coin1 = {
        'v_0': 80000,
        'v_1': 0,
        'v_2': 0
    }
    coin2 = {
        'v_0': 10000,
        'v_1': 100000,
        'v_2': 0
    }
    coin3 = {
        'v_0': 10000,
        'v_1': 0,
        'v_2': 100000
    }

    A = [coin1, coin2, coin3] # set of coins

    for it in range(2): # two iterations
        payment = abs(np.random.poisson(3000, 3)).tolist() # random payment values

        x, change = optimize(np.array(A), payment) # optimization routine

        A = [A[i] for i in range(len(x)) if x[i][0] != 1.0] # only keep not-selected coins in A

        A.append(change) # Add the change

        deposit = abs(np.random.poisson(3000, 3)).tolist()
        d = {
            'v_0': deposit[0],
            'v_1': deposit[1],
            'v_2': deposit[2],
        }
        A.append(d) # Add a deposit

优化函数 #

def optimize(A, payment):
    m = GEKKO()
    x = [m.Var(lb=0, ub=1, integer=True) for i in range(len(A))]

    input_tot = []
    input_tot_values = []

    # Calculate total input values for each currency
    for j in range(3):
        s = 0
        for i in range(len(A)):
            s = s + A[i][f'v_{j}'] * x[i]
        input_tot.append(s)
        input_tot_values.append(m.Intermediate(s)) # Intermediate s.t. we have evaluated expressions
    
    # Constraints that the total input should be larger than the payment
    m.Equation(equation=input_tot[0] > payment[0] + 1)
    m.Equation(equation=input_tot[1] > payment[1] + 1)
    m.Equation(equation=input_tot[2] > payment[2] + 1)

    # objective function
    obj = 1 / (1 + input_tot[0] - payment[0]) + 1 / (1 + input_tot[1] - payment[1]) + 1 / (1 + input_tot[2] - payment[2])

    m.Maximize(obj)
    m.options.SOLVER = 1
    m.solve()
    print(x)

    change = {
        'v_0': input_tot_values[0][0] - payment[0],
        'v_1': input_tot_values[1][0] - payment[1],
        'v_2': input_tot_values[2][0] - payment[2]
    }

    print('+++++++++++++++++++++++++++++++++++++++++')
    print(f'Payment: {payment}')
    print(f'Total Input: {input_tot_values}')
    print(f'Change: {change}')
    print('+++++++++++++++++++++++++++++++++++++++++')

    return x, change

问题现象 #

第二次迭代求解输出显示：

[[0.0], [0.0], [1.0]]
+++++++++++++++++++++++++++++++++++++++++
Payment: [2930, 2944, 3066]
Total Input: [[2949.0], [2967.0], [3019.0]]
Change: {'v_0': 19.0, 'v_1': 23.0, 'v_2': -47.0}
+++++++++++++++++++++++++++++++++++++++++

v_2货币的总输入（3019.0）小于支付额（3066），找零为负，说明Gekko求解器未满足约束条件。

问题排查与修复 #

核心原因 #

1. 约束定义错误：Gekko中>会被解析为非严格不等式≥，但代码中额外添加+1缩小了可行域，增加了无解概率。当求解器找不到可行解时，会返回违反约束的近似解。
2. 未设置严格约束容忍度：Gekko默认的容忍度允许微小约束违反，导致求解器返回不满足要求的解。
3. 硬币列表更新逻辑隐患：原代码中判断硬币是否被选中时，直接使用x[i][0]获取变量值，应该用x[i].value[0]，避免类型错误。
4. 迭代后硬币组合不足：第一次迭代后，原高额度硬币被移除，新增的找零和存款硬币额度较低，导致第二次迭代时没有足够的v_2货币满足支付需求，求解器只能返回违反约束的解。

修复措施 #

1. 修正约束条件：将约束改为符合业务需求的≥，去掉不必要的+1：

   # 替换原约束代码
   m.Equation(input_tot[0] >= payment[0])
   m.Equation(input_tot[1] >= payment[1])
   m.Equation(input_tot[2] >= payment[2])
   

2. 设置严格容忍度：在求解前添加以下代码，确保约束被严格满足：

   m.options.RTOL = 1e-6
   m.options.ATOL = 1e-6
   

3. 添加无解判断：求解后检查状态，若无解则抛出异常或处理：

   m.solve(disp=False)
   if m.options.APPSTATUS != 1:
       raise ValueError("No feasible solution available for this payment")
   

4. 修正硬币列表更新逻辑：正确获取Gekko变量的数值：

   A = [A[i] for i in range(len(x)) if x[i].value[0] != 1.0]
   

验证效果 #

修改后，求解器会严格满足约束条件，不会出现找零为负的情况。若确实无可行解，代码会抛出异常，便于后续处理。

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