PySCIPOpt 学习手册:Model、addCons 与求解诊断
简介:这份学习手册聚焦开源优化求解器 SCIP 的 Python 接口 PySCIPOpt,面向有一定 Python 基础、希望用代码求解混合整数规划与约束整数规划问题的运筹学学习者和工程实践者。内容以 Model 类的接口参考为主线,梳理初始化参数、对象比较与内存释放等方法,并分类整理了 addCons、addConsCardinality、addConsIndicator、addConsSOS1 等十余种约束添加接口的调用形式与 initial、separate、propagate 等可选参数含义,同时涉及 activateBenders、addBendersSubproblem 等 Benders 分解相关方法,便于按需查阅而非从头通读。压缩包内为 1 份 PDF 文档,体积约 28.69MB,内容以方法签名与参数说明的编排为主,适合当作手边的 API 速查材料。目前已有 1820 人学习,可帮助读者把建模思路顺利落地为可运行代码,减少在接口细节上反复试错的成本。
1. 从一份 Model 方法表说起:PySCIPOpt 解决的是哪一层问题
上个月有个做排产系统的朋友拿着两千个变量的模型来问:Python 业务层写得很顺,可自己搭的分支定界跑两小时还在原地打转。这类活儿不该用手写算法硬扛,PySCIPOpt 正是给这种场景准备的——它是 SCIP 求解器的 Python 接口,底层是一个支持约束整数规划的开源混合整数规划内核,上层让你用 addVar、addCons、optimize 这一套 API 把模型描述出来。这份学习手册的价值不在于教你线性规划入门,而在于把 Model 类的方法按功能拆开:构造与释放、变量与约束、求解状态查询,以及 include* 一整套插件注册。适合已经能用 Python 表达业务规则、需要把优化模型从“能跑”推到“能交付”的人;如果你还没写过一次 LP,先把建模教材过一遍再回来,收益大得多。
2. 环境到第一个可解模型:三步跑通 Model 生命周期
2.1 安装与版本对齐:pyscipopt 和 SCIP 内核的绑定关系
PySCIPOpt 不是纯 Python 包,它是对 SCIP 的 C 库做的 Cython 封装,所以“装上了”和“装对了”是两件事。现在多数平台直接一条命令就能装完,wheel 里已经带了编译好的内核;但如果你是从源码安装,或者本机之前装过独立版本的 SCIP,那 Python 包和 C 库的版本就可能对不上,症状通常不是导入时报错,而是运行时某个方法直接崩溃,这一类问题排查起来最费时间。
python -m pip install pyscipopt # 装完立刻验证:能构造 Model 且能报到内核版本,才算装对 python -c "from pyscipopt import Model; m=Model('probe'); m.printVersion()"我一般会再跑一段更细的自检,确认默认插件都注册上了:
from pyscipopt import Model m = Model("probe", defaultPlugins=True) print("stage:", m.getStage()) # 构造完成、尚未建模,预期为 'problem' print("readers:", m.getNReaders()) # 默认读取器数量,正常应大于 0 print("infinity:", m.infinity()) # 内核使用的无穷大表示getStage()返回的是模型对象所处的阶段:刚 new 出来是problem,读入或建好模型后进transformed,求解中在solving,结束后变solved。这个值后面会反复用到,因为很多接口只在特定阶段才允许调用。getNReaders()返回注册的读取器数量,如果它是 0,说明defaultPlugins没生效、插件注册被剥离过,那 addCons 之后的行为会很反直觉。
如果 import 阶段就报找不到 libscip,常见做法是清掉环境里重复的安装路径,用虚拟环境重新装一次,别在系统 Python 和 conda 环境之间来回混。
2.2 一个最小 MIP:变量、约束、目标构成的闭环
先看一个能从头跑到尾的完整例子,用一个五件物品的背包把 API 串起来。这里刻意选小规模,是为了让每一步的返回值都能手动核对。
from pyscipopt import Model, quicksum model = Model("knap_demo") # vtype: 'C' 连续 / 'I' 整数 / 'B' 二元;obj 是目标系数 x = {} for i in range(5): x[i] = model.addVar(name=f"x{i}", vtype="B", obj=-(i + 1) * 2) # 容量约束:sum((i+1) * x_i) <= 9 model.addCons(quicksum((i + 1) * x[i] for i in range(5)) <= 9, name="cap") model.setMaximize() # 把目标方向翻转为最大化 model.optimize() print("status:", model.getStatus()) print("obj:", model.getObjVal()) for i in range(5): print(i, model.getVal(x[i]))addVar(name='', vtype='C', lb=0.0, ub=None, obj=0.0, pricedVar=False)这几个参数的顺序值得记一下:ub=None表示上界取无穷,而不是取 0,这是新手最容易看反的地方;pricedVar是给列生成用的,普通建模保持 False,误设为 True 会让变量一开始就在列池里等着,求解行为会变。setMaximize()和setObjective(..., sense='maximize')效果上等价,区别在于前者只改方向、不动系数,后者会重写目标表达式,两者混用容易把系数搞丢。
getVal(x)只在有可行解时才给出数值,getObjVal()同理。如果getStatus()返回infeasible,这两个调用拿到的就是垃圾值或者直接抛异常,所以生产代码里必须先判状态再取值,别图省事一路 get 下去。
2.3 建模完成先做一次体检:getNVars / getNConss / checkSol
模型写完、还没 optimize 之前,值得花三十秒做个体检。数量统计是最便宜的哨兵:变量数、约束数、二元变量数对不上,说明循环或者条件分支写漏了。
print("vars:", model.getNVars(), "bin:", model.getNBinVars(), "int:", model.getNIntVars()) print("conss:", model.getNConss(), "rows:", model.getNLPRows()) if model.getNSols() > 0: sol = model.getBestSol() ok = model.checkSol(sol, printreason=True, original=True) print("feasible:", ok, "obj:", model.getSolObjVal(sol))checkSol的参数值得单独说:checkbounds校验变量界,checkintegrality校验整数性,checklprows校验 LP 行是否满足,completely=False表示遇到第一个越界就返回而不是全查一遍;original=True让你按原始问题(未经 presolve 转换)检查,跟业务口算的结果能对上。我用它的场景通常是:外部启发式算了一个解,想丢回模型里继续用,丢之前必须先 check,不然一个越界解进去会把后续所有节点带偏。
提示:
getNConss在 presolve 前后统计范围会变,别拿优化前的数字和优化后的报告直接对比,容易误判成“约束丢了”。
3. addCons 族的 11 个标志位:每段求解流程各管哪一部分
3.1 通用标志位:initial、separate、enforce、check、propagate 的分工
addCons(cons, name='', initial=True, separate=True, enforce=True, check=True, propagate=True, local=False, modifiable=False, dynamic=False, removable=False, stickingatnode=False)这一长串参数,第一次看确实像天书。它们其实是 SCIP 求解流程的分段开关:建模期、LP 构造期、分离期、强制期、检查期、传播期,每个标志位对应问一句“这条约束参不参加这一阶段”。
| 标志位 | 默认 | 作用阶段 | 什么时候需要改动 |
|---|---|---|---|
| initial | True | 初始 LP 构造 | 约束量极大、想延迟到需要时再加 |
| separate | True | 分离/割平面生成 | 只有自定义约束处理器才真正用得上 |
| enforce | True | 强制阶段(找可行解、证不可行) | 只希望它当割用、不参与完整性判定时 |
| check | True | 解的可行性检查 | 调试期千万别关,关掉等于蒙眼 |
| propagate | True | 域收缩传播 | 存在自写传播器时统一管理 |
| local | False | 作用域 | 分支决策后只在子树内成立的约束 |
| modifiable | False | 求解中修改 | 要调 chgLhs / chgRhs / addConsCoeff |
| dynamic | False | 可删除后重建 | 大量临时割、内存吃紧 |
| removable | False | 可删除 | 同上,割池管理的常用开关 |
| stickingatnode | False | 节点归属 | 与特定节点强绑定的约束 |
真正需要日常关心的是三个:initial、local、modifiable。initial=False在大模型上有效,能把 LP 规模压下来,让求解器从一个小松弛起步、逐步补约束;但副作用是首次 LP 松弛会很松,对偶界起步差,如果你的时间预算很短,反而更糟。local=True只对分支后的子节点有意义,写在根节点上等于白写——这是后面避坑章节里要展开的一条。
modifiable是最容易踩雷的一个。它标记的是“我打算在求解过程中改这条约束”。没标记就调用修改接口,行为是未定义的,轻则改动被忽略,重则内核状态错乱。
c = model.addCons(x[0] + x[1] <= 5, name="c_mod", modifiable=True) model.enableReoptimization() # 想改完再解,这个开关必须在首次 optimize 之前打开 model.optimize() model.chgRhs(c, 6) # 改右端项 model.chgReoptObjective({x[0]: 3.0}, sense="minimize") # 换目标 model.optimize() # 走重优化路径,通常比从头解快enableReoptimization()是这套用法成立的前提,而且必须在第一次optimize()之前调用,事后再补是无效的。重优化路径也不是万能的,它更适合“同一结构、系数小幅变动”的场景;如果每轮都新增删变量,收益会被抵消,还不如直接重建模型。
3.2 逻辑与结构约束:And / Or / Xor / Indicator / Cardinality / SOS
这批接口是手册里最省事的一节——很多建模书里要写三四行大 M 的约束,这里一行就够。但省事的代价是你要清楚它背后生成了什么。
a = model.addVar(vtype="B", name="a") b = model.addVar(vtype="B", name="b") r = model.addVar(vtype="B", name="r") model.addConsAnd([a, b], r, name="and1") # r = a AND b,双向等价 model.addConsOr([a, b], r, name="or1") # r = a OR b model.addConsXor([a, b], r, name="xor1") # r = a XOR b,奇偶异或 # 指示约束:binvar 取 activeone 时 cons 成立;取另一值时 cons 被放松,而不是取反 model.addConsIndicator(a + b >= 1, binvar=a, activeone=True, name="ind1") # 基数约束:{a,b} 中恰好 cardval 个为 1 model.addConsCardinality([a, b], 1, name="card1") # SOS1:集合中至多一个非零,由分支器直接处理 model.addConsSOS1([a, b], weights=[1, 2], name="sos1")三处容易被误解的地方。第一,And / Or / Xor 是等价关系,不是单向蕴含,内核会额外引入辅助变量和多个线性化约束,变量数会涨,别指望它零成本。第二,addConsIndicator里activeone=True表示binvar=1时约束生效;binvar=0时约束被放松(不再限制),很多人误以为取反后另一支也受约束,结果模型比脑补的松,解出来才发现在钻空子。第三,SOS1 / SOS2 在 LP 松弛里并不是线性行,靠的是分支器特殊处理,所以它的松弛界通常比等价的大 M 写法更紧,但前提是大 M 够小——这也是为什么很多人宁可用 SOS 不用大 M。
addConsLocal和addConsNode属于另一类:前者绑定当前节点的子树作用域,后者需要一个节点对象,只能在求解回调里构造。放在建模主流程里调用它们,拿不到预期效果。
3.3 增量修改的三条路:批量添加、改系数、改左右端
模型改动有三种粒度,选错了接口是常见的时间黑洞。
# 批量添加,返回约束对象列表,后面可以单独操作 cs = model.addConss([x[0] + x[1] <= 4, x[1] + x[2] <= 4], name="pair") # 只对线性约束有效:给第 0 条约束补上 x2 的系数 model.addConsCoeff(cs[0], x[2], 1.0) # 改左右端 model.chgLhs(cs[1], 1.0) model.chgRhs(cs[0], 5.0)addConss适合在循环里批量落约束,比逐条 addCons 少一层 Python 调用开销,更重要的是返回的对象列表能让你后续按索引定位,不用靠名字反查。addConsCoeff的前提是这条约束是线性的、且创建时标了modifiable=True,对二次约束或非线性表达式无效,硬用会抛类型错误。chgLhs/chgRhs则要留意方向:把一条<=约束的 lhs 从-inf改成具体值,相当于把它变成区间约束,逻辑上是新增了限制,不是放松。
改完约束之后必须重新optimize(),而且要靠enableReoptimization()才能复用之前的信息,否则内核会退回全量重解,这一点在批量调参的脚本里影响很大。
4. 求解诊断:用 get* 一族把 timelimit 和 infeasible 分开
4.1 三级状态判读:getStage / getStatus / getLPSolstat
结果读错,比模型建错更隐蔽。PySCIPOpt 里至少有三个“状态”,职责完全不同,混着看就会把“没算完”错认成“无解”。
| 接口 | 典型取值 | 说明 |
|---|---|---|
| getStage | problem / transformed / presolving / solving / solved | 模型对象的生命周期阶段 |
| getStatus | optimal / infeasible / unbounded / timelimit / nodelimit / unknown | 本次求解的结论 |
| getLPSolstat | optimal / infeasible / notsolved | 最近一次 LP 松弛的状态 |
我习惯封一个小函数把三个值一次打出来,并顺带决定后续动作:
def solve_report(model): stage, status = model.getStage(), model.getStatus() nsols = model.getNSols() print(f"stage={stage} status={status} nsols={nsols}") if nsols > 0: print("primal:", model.getPrimalbound(), "dual:", model.getDualbound(), "gap:", model.getGap()) else: print("no incumbent; lpstat:", model.getLPSolstat()) return status这里的关键判断是nsols:status为timelimit或unknown时,模型可能已经拿到了一个不错的可行解,只是没证明最优;此时不读解就白白浪费了。反过来,nsols == 0且status == 'infeasible'才是真不可行。getStatus()返回的是小写字符串,别拿它和数字常量做比较,这类笔误在脚本里特别常见。
4.2 上下界与 gap:getObjVal / getPrimalbound / getDualbound 的正确读法
getObjVal()是最优解的值,只有在status == 'optimal'时才等于最优;其余情况返回的是当前最好可行解的值。想区分“当前解”和“界”,要用getPrimalbound()和getDualbound()。二者的差就是 gap,getGap()给的是相对 gap。
pb, db = model.getPrimalbound(), model.getDualbound() if model.isInfinity(abs(pb)) or model.isInfinity(abs(db)): print("bound is infinite, keep solving or fix the model") else: print("abs gap:", abs(pb - db), "rel gap:", model.getGap())isInfinity(value)是官方提供的判定,比直接和1e20比较靠谱,因为内核的无穷表示是个大常数、不是float('inf')。isFeasEQ / isFeasLT / isFeasZero这一族同理,都是带容差的比较,业务代码里判“两个目标值是否相等”时应该用它们,而不是==。
想要提前收工,setObjlimit(objlimit)是常用手段:给一个“比当前可行解更好的目标阈值”,求解器一旦证明达不到就停。这在滚动排产里很实用——只关心有没有比现行方案更好的解,不关心最优解长什么样。
4.3 对偶、松弛与约化成本:先确认这些值什么时候才有效
这一节是整份手册里最容易误用的部分。getDualsolLinear(cons)取的是约束的对偶值,但它对应的是内核最后一次求解的 LP——在 MIP 场景下,那可能是某个节点上被割平面加强过的松弛,不是你在建模时写的那条原始约束。拿这个数做灵敏度分析,结论是不成立的。
# 只有在连续问题(全部变量 vtype='C')下,对偶值才对应原模型 if model.getNIntVars() == 0 and model.getStatus() == "optimal": for name, c in [("cap", cons_cap)]: print(name, "dual:", model.getDualsolLinear(c), "slack:", model.getSlack(c))需要 MIP 场景下的对偶信息时,我一般会复制一份模型、把变量类型全改成连续、单独解一次 LP,拿这份干净结果做分析;也可以调用getDualSolVal(cons, boundconstraint=False)取约束层面的对偶,boundconstraint=True时取的是变量界对应的那部分,做影子价格时这个区分很重要。
getVarRedcost(var)是约化成本,同样只在 LP 解有效时有意义。getLPBInvRow / getLPBInvARow这一族给的是基逆的行,做割的有效性判断和敏感性分析会用到,但它们按 LP 的行列索引取值,一旦重新 optimize 索引就变了,缓存下来是没用的。
5. 五类翻车现场排查:从 import 失败到目标值对不上
5.1 先分清是建模错还是求解错
排查的第一步不是看报错,而是看模型本身对不对。我通常按这个顺序:先把模型写成文件落盘,再打开统计信息,最后用一个小规模实例验证手算结果。落盘这一步很多人省掉,但它是最快的分界线——文件里能看出来的问题,都是建模问题;文件里看不出来的,才是求解问题。
model.writeProblem("case_a.lp") # 人类可读,适合肉眼核对 model.writeProblem("case_a.mps") # 通用格式,方便换求解器交叉验证 model.setLogfile("run_a.log") # 把日志固定下来,便于事后比对 # ... optimize ... model.printStatistics() # 节点数、割数、presolve 效果一览.lp文件里能直接看到每条约束的系数和左右端,很多时候“结果不对”就是某个循环里变量索引写错了一位,肉眼一扫就出来了。printStatistics()则告诉你 presolve 砍掉了多少行、根节点 LP 值是多少、割加了多少——如果 presolve 后约束数归零,说明模型退化,先回去查约束是不是全被吸收了。
5.2 五条现象、原因与解决办法
现象一:import pyscipopt 报找不到 libscip 或者段错误。原因几乎都是环境里存在多份来源不同的安装,Python 包找的是这份 C 库、系统里装的是另一份,版本不匹配。解决:在干净的虚拟环境里重装,pip show pyscipopt确认只有一处路径;conda 和 pip 不要混用同一个环境。装完立刻跑 2.1 里那段自检代码,别等到写完模型才发现。
现象二:约束加了,但解出来根本不满足它。最典型的原因是把约束表达式写进了条件语句,或者用 Python 的and/or连接了两个约束。x[0] + x[1] <= 5在 PySCIPOpt 里重载成了表达式对象,可一旦放进if的判断位置,Python 就会尝试求它的布尔值,结果要么直接抛异常,要么逻辑走偏——你写的是约束,解释器当成了真假判断。解决:所有约束表达式只作为addCons*的入参出现,绝不进条件语句;需要逻辑组合时用addConsAnd / addConsOr / addConsXor,不要用 Python 的布尔运算符。
现象三:optimize 之后 getVal 拿不到数,或者解释器直接崩。排查顺序是“状态、解、生命周期”。先确认getNSols() > 0,再取解;如果调过freeProb()或freeTransform(),之后再访问任何变量、约束对象都是悬空引用,崩溃是必然的。addSol(sol, free=True)也有同样的陷阱:free=True表示所有权交给内核,之后你再调freeSol(sol)就是重复释放。解决:把释放放进try/finally,释放后立刻把引用置None,不要在释放和退出之间夹带任何模型操作。
现象四:status 是 unknown,团队据此判断模型不可行。unknown通常意味着求解被限制截断了,可能是时间、节点、也可能是内存,而不是模型无解。判据是nsols:大于 0 说明有可行解,只是没证明最优;等于 0 且 LP 状态一路不可行,才谈得上无解。解决:把limits/time、limits/nodes、limits/gap三个值打进日志,和 status 一起记录;出报告时永远写“在 X 秒预算内未找到可行解”,而不是“模型无解”。
现象五:目标值和手算差一个常数。多半是目标偏移或者重复设目标造成的。addObjoffset(offset, solutions=False)加的是常量项,solutions=False时它不参与解的比较,容易被忘掉;而setObjective(coeffs, sense, clear='true')的clear参数控制是否先清空已有目标系数,连续调用两次而不清空,系数会叠加。解决:目标只在建模末尾设一次,需要改时显式传clear='true';用getObjoffset(original=True)把偏移单独打出来核对。
5.3 把一次求解钉死成可回放的文件
线上出的问题,复现不了就等于没修。我的习惯是每次关键求解都留下三样东西:模型文件、日志文件、参数快照。
import json model.writeProblem("repro/case.lp") model.setLogfile("repro/case.log") params = model.getParams() # 返回当前所有参数的字典 with open("repro/params.json", "w", encoding="utf-8") as f: json.dump({k: str(v) for k, v in params.items()}, f, ensure_ascii=False, indent=1)getParams()返回的是全量参数字典,条目很多,全量保存看起来笨,但恢复时最省心:用readParams(file)一次性灌回去,或者用setParams(dict)批量设,比逐条setParam可控。参数快照的另一个用途是版本对比——两次结果差异很大时,先 diff 参数文件,往往能直接看出是哪次调参改出了事。
6. 进阶:把 SCIP 当插件框架,用 include* 换掉一个求解环节
手册里includeHeur / includeSepa / includeBranchrule / includeEventhdlr / includeConshdlr这一大串,初看像是给内核开发者准备的,其实普通项目也用得上——最典型的是事件钩子:你想在每次找到更好的可行解时记录时间戳、上报监控、或者提前剪掉业务上不可接受的解,用事件比轮询强得多。
from pyscipopt import Model, Eventhdlr, SCIP_EVENTTYPE class IncumbentTracker(Eventhdlr): def eventinit(self): # 注册感兴趣的事件类型,必须和 eventexit 成对 self.model.catchEvent(SCIP_EVENTTYPE.BESTSOLFOUND, self) def eventexit(self): self.model.dropEvent(SCIP_EVENTTYPE.BESTSOLFOUND, self) def eventexec(self, event): print("new incumbent:", self.model.getPrimalbound(), "at", self.model.getSolvingTime()) m = Model("with_event") m.includeEventhdlr(IncumbentTracker(), "inctrk", "track incumbent improvements")三个要点:eventinit和eventexit必须成对,漏掉 drop 会在多次 optimize 之后泄漏处理器;eventexec里只做轻量动作,打印可以,重建模型不行,回调期间改结构会破坏内核状态;getSolvingTime()在回调里是安全的,但getObjVal()未必,取界用getPrimalbound()更稳。
参数这一侧,与其记一堆名字,不如记住五个真正影响时间预算的:
| 参数名 | 类型 | 常用取值 | 作用 |
|---|---|---|---|
| limits/time | real | 60 / 600 | 总墙钟时间上限 |
| limits/gap | real | 1e-4 | 相对 gap 达到即停 |
| limits/nodes | longint | 100000 | 节点数上限,压测时用来定规模 |
| presolving/maxrounds | int | 0 / -1 | 0 关预求解,-1 用默认轮数 |
| numerics/feastol | real | 1e-6 | 可行性容差,数值病态时适当放宽 |
注意:
limits/gap设成 0 不是“不限制”,而是要求相对 gap 严格为 0,等于逼求解器证明最优;压测时随手写 0,跑一晚上不出结果很常见。想让它在“够好就停”,写 1e-4 这类正数。
setEmphasis()是另一条捷径,feas偏向快速找可行解、opt偏向证明最优、numerics偏向数值稳定,适合在还没摸清参数之前先定个大方向,再逐条微调。至于includeConshdlr自写约束处理器,那是另一个量级的工作——你要自己实现分离、强制、检查、传播四套回调,收益只在“现成的线性化写法太松”时才明显,普通业务模型用addConsIndicator加合适的界就够了。
说个教训。早年我图省事,把所有临时割都设成removable=True又不开modifiable,跑小实例一切正常,规模一上去就开始出现“同一份数据两次跑出不同目标值”。查了两天才定位到是割在节点之间被删掉又重建,而我的约束里带着 Python 侧的状态。从那以后我每次改约束属性,都强制走一遍“小实例对照 + 参数快照 + 日志落盘”这三步,哪怕只改一个布尔值。希望这套流程也能帮你少走几个晚上。
本文还有配套的精品资源,点击获取