基于Gekko的混合整数非线性规划求解提速策略咨询

问题概述 #

当前求解的制造设施设计优化问题属于混合整数非线性规划（MINLP）类问题，核心要素如下：

• 离散设备选型参数α：需从集合{0.5, 1, 2, 5}中选取唯一值，为特殊有序集（SOS1）类型变量
• 零件生产数量β：共n个离散整数变量，取值范围为1~5的整数
• 固定输入：成本矩阵M，待建设单元总规模n
• 约束逻辑：对β做对数变换得到δ=log(β)，计算矩阵乘积Mδ，需满足所有维度下α*Mδ_i > 50
• 优化目标：最小化α取值

问题最简可复现代码如下：

import numpy as np
from gekko import GEKKO
m = GEKKO()
n = 100
# select one of the Special Ordered Set
α = m.sos1([0.5, 1, 2, 5])
# integer variables
β = m.Array(m.Var,n,integer=True,lb=1,ub=5)
# log transform
δ = [m.log(b) for b in β] 
# matrix
M = np.random.rand(n,n)
# constraints
Mδ = M.dot(δ)
m.Equations([α*x>50 for x in Mδ])
# objective
m.Minimize(α)
# solve
m.solve()
print('α: ', α.value[0])
print('β: ', β)

当n=100时，问题可行解总规模达到5^100 * 4 ≈ 3e70，完全不具备穷举可行性。基于GEKKO+APOPT求解器的基准测试中，n=100规模算例求解耗时约33.4秒，得到最优解α=1.0，求解日志如下：

Number of state variables:            205
 Number of total equations: -          102
 Number of slack variables: -          100

 ---------------------------------------------------
 Solver         :  APOPT (v1.0)
 Solution time  :    33.4135999999999      sec
 Objective      :    1.00000000000000     
 Successful solution
 ---------------------------------------------------
 
α:  1.0
β:  [[3.0] [4.0] [4.0] [5.0] [4.0] [3.0] [3.0] [3.0] [3.0] [5.0] [4.0] [4.0]
 [4.0] [4.0] [3.0] [3.0] [4.0] [4.0] [3.0] [3.0] [4.0] [3.0] [3.0] [4.0]
 [4.0] [4.0] [5.0] [3.0] [1.0] [4.0] [4.0] [4.0] [3.0] [4.0] [4.0] [3.0]
 [3.0] [3.0] [4.0] [4.0] [3.0] [4.0] [1.0] [3.0] [4.0] [3.0] [3.0] [4.0]
 [3.0] [3.0] [4.0] [4.0] [3.0] [1.0] [3.0] [3.0] [3.0] [3.0] [3.0] [3.0]
 [4.0] [3.0] [3.0] [4.0] [4.0] [4.0] [2.0] [3.0] [1.0] [4.0] [4.0] [3.0]
 [3.0] [3.0] [4.0] [4.0] [3.0] [3.0] [4.0] [3.0] [3.0] [4.0] [3.0] [3.0]
 [3.0] [2.0] [4.0] [5.0] [4.0] [2.0] [3.0] [4.0] [3.0] [4.0] [1.0] [3.0]
 [5.0] [3.0] [4.0] [5.0]]

可行提速策略 #

模型重构类（提速效果最显著，可实现量级提升） #

• 按目标单调性枚举α取值，拆分问题为独立可行性校验问题
  由于优化目标是最小化α，且α只有4个按从小到大排序的固定可选值，完全不需要将α作为SOS1变量放入模型联合求解：按0.5、1、2、5的顺序依次固定α取值，单独求解对应β的整数可行性问题即可。只要找到某个α下存在可行β，该解就是全局最优解，直接终止计算即可。
  *该方案不存在次优问题，比“找到首个可行解就终止”的通用策略可靠性更高，同时由于拆分后的问题不需要处理SOS1变量的分支逻辑，求解速度会明显提升。
• 消除模型非线性项，将MINLP转化为MILP
  原模型的非线性完全来自β的对数变换，但β的取值只有1、2、3、4、5这5个固定整数，对应的log值是可以提前计算的常数，不需要求解器在迭代中做非线性对数运算。
  重构方法：对每个β_i引入5个0-1变量y_i1,y_i2,y_i3,y_i4,y_i5，添加约束sum(y_ik for k in 1~5) = 1，则β_i = sum(k*y_ik)，δ_i = sum(log(k)*y_ik)，此时δ是0-1变量的线性组合，Mδ和所有约束都变为线性形式，原问题完全转化为混合整数线性规划（MILP）问题，求解效率远高于MINLP。
• 提前剪枝不可行的α取值
  结合β的取值范围（1~5），可以提前计算每个Mδ_i的理论上下界，直接排除不可能满足约束的α取值。例如α=0.5时，约束要求所有Mδ_i>100，如果Mδ_i的理论最大值都小于100，直接跳过α=0.5的求解流程，不需要进入求解器。

求解器参数调优类 #

• 调整分支优先级
  将对约束和目标影响大的变量设为高分支优先级：首先是α相关变量（如果保留联合求解逻辑），其次是M矩阵中权重高的位置对应的β变量，让分支定界过程优先搜索关键变量，大幅减少需要探索的分支节点数量。
• 设置合理的终止规则
  如果工程场景允许一定的最优性间隙，可以直接设置APOPT的最优性容差m.options.OTOL为可接受的值（例如0.05代表接受5%以内的间隙），也可以设置m.options.MAX_TIME为可接受的最长求解时间，到达阈值后直接输出当前找到的最优解。
• 换用适配MILP的求解器
  完成线性化重构后，可换用CBC、Gurobi（商业授权）等专门面向MILP的求解器，这类求解器内置的割平面、分支启发式逻辑比APOPT更成熟，n=100规模的问题通常可实现秒级求解。

工程近似类 #

• 利用矩阵结构分块求解
  如果成本矩阵M存在稀疏性、块对角结构，可以将100维的原问题拆分为多个独立的小规模子问题分别求解，求解复杂度会随拆分规模指数下降。
• 注入高质量初值
  求解前给β赋一个接近可行域的经验初值（例如所有β先取4或5），让求解器从可行域附近开始搜索，减少前期对不可行节点的探索时间。

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