如果模型能够顺利求解,通过一组数据测试,通常是需要烧高香的事情。
更常发生的事情是模型解不了,然后我们需要挠破头皮去找哪些约束出了问题,还是输入数据出现异常,还是变量上下界不对,还是遇到了玄学。
CPLEX提供了冲突检查的功能,能够在某些情况下稍微帮助我们定位可能的问题所在。
这里我们通过一个简单小例子看一下冲突检测怎么整。
例子和建模
临近下班了,急着回家,这里就简单一点,直接拿CPLEX manual里的一个模型举例。
# 导入库
from docplex.mp.model import Model
from docplex.mp.model_reader import ModelReader
from collections import defaultdict
import os
# 直接从lp文件读取模型
lp_string = """
Minimize
obj: cost
Subject To
c1: - cost + 80 x1 + 60 x2 + 55 x3 + 30 x4 + 25 x5 + 80 x6 + 60 x7 + 35 x8
+ 80 x9 + 55 x10 = 0
c2: x1 + x2 + 0.8 x3 + 0.6 x4 + 0.4 x5 >= 2.1
c3: x6 + 0.9 x7 + 0.5 x8 >= 1.2
c4: x9 + 0.9 x10 >= 0.8
c5: 0.2 x2 + x3 + 0.5 x4 + 0.5 x5 + 0.2 x7 + 0.5 x8 + x10 - service = 0
c6: x1 + x6 + x9 >= 1
c7: x1 + x2 + x3 + x6 + x7 + x9 >= 2
c8: x2 + x3 + x4 + x5 <= 0
c9: x4 + x5 + x8 <= 1
c10: x1 + x10 <= 1
Bounds
service >= 3.2
Binaries
x1 x2 x3 x4 x5 x6 x7 x8 x9 x10
End
"""
# cplex没有提供直接读取string的接口,不得不进行个文件暂存操作
temp_input_file_name = "temp.lp"
with open(temp_input_file_name, "w") as f:
f.write(lp_string)
mdl = ModelReader.read(temp_input_file_name, model_name = "InfeasibelLP")
# 清除临时文件
if os.path.exists(temp_input_file_name):
os.remove(temp_input_file_name)
建完模型之后我们就可以尝试求解了,为了对比一下求解中输出的信息和冲突检测输出的信息,我们打印求解过程:
# 尝试求解
sol = mdl.solve(log_output=True)
得到的求解信息如下:
Version identifier: 12.10.0.0 | 2019-11-26 | 843d4de
CPXPARAM_Read_DataCheck 1
Row 'c8' infeasible, all entries at implied bounds.
Presolve time = 0.00 sec. (0.00 ticks)
Root node processing (before b&c):
Real time = 0.00 sec. (0.01 ticks)
Parallel b&c, 16 threads:
Real time = 0.00 sec. (0.00 ticks)
Sync time (average) = 0.00 sec.
Wait time (average) = 0.00 sec.
------------
Total (root+branch&cut) = 0.00 sec. (0.01 ticks)
求解器提示我们 c8
不可行,我们可以打印所有约束出来看一眼:
print("----------------")
print("模型中共有{}条约束".format(mdl.number_of_constraints))
for ct in mdl.iter_constraints():
print(ct)
----------------
模型中共有10条约束
c1: -cost+80x1+60x2+55x3+30x4+25x5+80x6+60x7+35x8+80x9+55x10 == 0
c2: x1+x2+0.800x3+0.600x4+0.400x5 >= 2.1
c3: x6+0.900x7+0.500x8 >= 1.2
c4: x9+0.900x10 >= 0.8
c5: 0.200x2+x3+0.500x4+0.500x5+0.200x7+0.500x8+x10-service == 0
c6: x1+x6+x9 >= 1.0
c7: x1+x2+x3+x6+x7+x9 >= 2.0
c8: x2+x3+x4+x5 <= 0
c9: x4+x5+x8 <= 1.0
c10: x1+x10 <= 1.0
于是我们看到了c8:x2+x3+x4+x5 <= 0
,然后呢?为啥不可行??它和哪个其他约束冲突了???
这个模型比较小我们还可以手动试算一下,工业模型的规模能让人原地升天。
冲突检测
笔者没有在docplex模块中找到冲突检测的API,因此我们还是不得不借助cplex类提供的接口来进行冲突检测,示例代码如下:
# 获取cplex.Cplex()类对象
c = mdl.cplex
# 进行冲突检测
c.conflict.refine(c.conflict.all_constraints())
# 输出检测信息,再重新读入并在控制台中输出
# 需要吐槽的是Cplex并不支持以IOSteam为对象输入输出,因此不得不反复建立临时文件
output_fname = r"conflict.txt"
c.conflict.write(output_fname)
with open(output_fname, "r") as f:
print(f.read())
得到如下的检测信息:
\ENCODING=ISO-8859-1
\Problem name: _conflict
Minimize
obj1:
Subject To
c2: x1 + x2 + 0.8 x3 + 0.6 x4 + 0.4 x5 >= 2.1
c8: x2 + x3 + x4 + x5 <= 0
Bounds
0 <= x1 <= 1
0 <= x2 <= 1
0 <= x3 <= 1
0 <= x4 <= 1
0 <= x5 <= 1
Binaries
x1 x2 x3 x4 x5
End
这里我们看到,是约束c2
和c8
之间产生了冲突,我们的变量x1
- x5
都是binary variable,因此要满足约束c8
,只能让x2, x3, x4, x5
都取0,而x1
的上界是1,因此不能满足c2: x1>=2.1
。
急着回家就先这样吧,实际使用的时候通常还需要根据约束类型、变量类型、变量上下界等等对conflict groups进行筛选,再在各个组内进行检测,得到更加精细的冲突信息,读者可以自行研究。