手把手教你用Python调用SCIP求解最优化模型

作者：刘兴禄，清华大学，清华-伯克利深圳学院博士在读

上一期我们介绍了Python调用SCIP时，SCIP的安装与配置，这一期我们来介绍使用Python调用SCIP完成建模和求解，以及解的输出等。

一个简单的例子

我们用SCIP给出的官方案例作为引入，先让大家大致看一下Python调用SCIP求解最优化模型的大题流程。

我们以最经典的对称旅行商问题(asymmetric traveling salesman problem, ATSP)为例。其具体模型为

min ⁡   ∑ i ∑ j c i j x i j ∑ i ∈ V x i j = 1,  ∀ j ∈ V , i ≠ j ∑ j ∈ V x i j = 1,  ∀ i ∈ V , i ≠ j μ i − μ j + N x i j ⩽ N − 1,  ∀ i , j ∈ V , i , j ≠ 0, i ≠ j x i j ∈ { 0, 1 } , μ i ⩾ 0,  μ i ∈ R 1 ∀ i ∈ V \begin{aligned} \min\text{ }\sum_i{\sum_j{c_{ij}x_{ij}}} \\ \sum_{i\in V}{x_{ij}}=\text{1, } & \hspace{1cm}\forall j\in V,i\ne j \\ \sum_{j\in V}{x_{ij}}=\text{1, } & \hspace{1cm}\forall i\in V,i\ne j \\ \mu _i-\mu _j+Nx_{ij}\leqslant N-\text{1, } & \hspace{1cm}\forall i,j\in V,i,j\ne \text{0,}i\ne j \\ x_{ij}\in \left\{ \text{0,}1 \right\} , \mu _i\geqslant \text{0, }\mu _i\in \mathbf{R}^1 &\hspace{1cm}\forall i\in V \end{aligned} min i∑​j∑​cij​xij​i∈V∑​xij​=1, j∈V∑​xij​=1, μi​−μj​+Nxij​⩽N−1, xij​∈{0,1},μi​⩾0, μi​∈R1​∀j∈V,i=j∀i∈V,i=j∀i,j∈V,i,j=0,i=j∀i∈V​

关于TSP更为细致的介绍，请参照我的另外一篇文章

https://blog.csdn.net/HsinglukLiu/article/details/107848461
运筹学修炼日记：TSP中两种不同消除子环路的方法及callback实现（Python调用Gurobi求解，附以王者荣耀视角解读callback的工作逻辑）

Python调用SCIP求解ATSP的具体代码如下。

##@file atsp.py
#@brief solve the asymmetric traveling salesman problem

"""
formulations implemented:
    - mtz -- Miller-Tucker-Zemlin's potential formulation
Copyright (c) by Joao Pedro PEDROSO and Mikio KUBO, 2012
"""
from pyscipopt import Model, quicksum, multidict

def mtz(n,c):
    """mtz: Miller-Tucker-Zemlin's model for the (asymmetric) traveling salesman problem
    (potential formulation)
    Parameters:
        - n: number of nodes
        - c[i,j]: cost for traversing arc (i,j)
    Returns a model, ready to be solved.
    """

    model = Model("atsp - mtz")

    x,u = {},{}
    for i in range(1,n+1):
        u[i] = model.addVar(lb=0, ub=n-1, vtype="C", name="u(%s)"%i)
        for j in range(1,n+1):
            if i != j:
                x[i,j] = model.addVar(vtype="B", name="x(%s,%s)"%(i,j))

    for i in range(1,n+1):
        model.addCons(quicksum(x[i,j] for j in range(1,n+1) if j != i) == 1, "Out(%s)"%i)
        model.addCons(quicksum(x[j,i] for j in range(1,n+1) if j != i) == 1, "In(%s)"%i)

    for i in range(1,n+1):
        for j in range(2,n+1):
            if i != j:
                model.addCons(u[i] - u[j] + (n-1)*x[i,j] <= n-2, "MTZ(%s,%s)"%(i,j))

    model.setObjective(quicksum(c[i,j]*x[i,j] for (i,j) in x), "minimize")
    
    model.data = x,u
    return model


def sequence(arcs):
    """sequence: make a list of cities to visit, from set of arcs"""
    succ = {}
    for (i,j) in arcs:
        succ[i] = j
    curr = 1    # first node being visited
    sol = [curr]
    for i in range(len(arcs)-2):
        curr = succ[curr]
        sol.append(curr)
    return sol


if __name__ == "__main__":
    n = 5
    c = { (1,1):0,  (1,2):1989,  (1,3):102, (1,4):102, (1,5):103,
          (2,1):104, (2,2):0,  (2,3):11,  (2,4):104, (2,5):108,
          (3,1):107, (3,2):108, (3,3):0,  (3,4):19,  (3,5):102,
          (4,1):109, (4,2):102, (4,3):107, (4,4):0,  (4,5):15,
          (5,1):13,  (5,2):103, (5,3):104, (5,4):101, (5,5):0,
         }

    model = mtz(n,c)
    model.hideOutput() # silent mode
    model.optimize()
    cost = model.getObjVal()
    print()
    print("Miller-Tucker-Zemlin's model:")
    print("Optimal value:", cost)
    #model.printAttr("X")
    for v in model.getVars():
        if model.getVal(v) > 0.001:
            print(v.name, "=", model.getVal(v))

    x,u = model.data
    sol = [i for (p,i) in sorted([(int(model.getVal(u[i])+.5),i) for i in range(1,n+1)])]
    print(sol)
    arcs = [(i,j) for (i,j) in x if model.getVal(x[i,j]) > .5]
    sol = sequence(arcs)
    print(sol)

运行结果为

Miller-Tucker-Zemlin's model:
Optimal value: 330.0
x(1,4) = 1.0
x(2,3) = 1.0
x(3,5) = 1.0
x(4,2) = 1.0
x(5,1) = 1.0
u(5) = 4.0
u(2) = 2.0
u(3) = 3.0
u(4) = 1.0

[1, 4, 2, 3, 5]
[1, 4, 2, 3]

接下来我们来详细介绍具体的用法。

Python调用SCIP求解最优化模型的一般步骤

使用Python调用SCIP求解最优化模型的一般步骤一般分为

1. 创建模型对象: pyscipopt.scip.Model(modelname)
2. 创建决策变量：addVar(name='', vtype='C', lb=0.0, ub=None, obj=0.0, priceVar=False)
3. 创建目标函数：pyscipopt.scip.Model.setObjective(coeffs, sense='minimize', clear='true)
4. 创建约束条件：addCons()、addConsAnd()、addConsOr()、addConsXor() 等
5. 求解模型：model.optimize()
6. 获得解的信息并输出：model.getBestSol()、model.getObjVal()、model.getVal(z)等。

接下来我们来逐步介绍。

创建模型对象

调用过程中，所有的组成部分，包括决策变量，约束等，都是存储在一个模型对象中的，也就是一个pyscipopt.scip.Model类的对象。需要用到的函数为

• pyscipopt.scip.Model(modelname) : creat and return a pyscipopt.scip.Model object

下面举一个例子

from pyscipopt import Model 
model = Model("puzzle") 

创建决策变量

创建决策变量，一般用函数

• addVar(name='', vtype='C', lb=0.0, ub=None, obj=0.0, priceVar=False)

其中各个参数的详细解读如下：

下面我们给出几个例子

from pyscipopt import Model 
model = Model() 

# Create decision variables : example 
x = model.addVar(vtype="I", name="octopusses")
y = model.addVar(vtype="I", name="turtles")
z = model.addVar(vtype="I", name="cranes")

如果是较多的变量，我们可以用字典进行存储，比如上面的ATSP的例子。

x,u = {},{}
for i in range(1,n+1):
    u[i] = model.addVar(lb=0, ub=n-1, vtype="C", name="u(%s)"%i)
    for j in range(1,n+1):
        if i != j:
            x[i,j] = model.addVar(vtype="B", name="x(%s,%s)"%(i,j))

之后，如果我们想提取变量的名字等信息，可以使用

• var.name
  

设置目标函数

设置目标函数，一般使用函数

• pyscipopt.scip.Model.setObjective(coeffs, sense='minimize', clear='true)

具体的参数解释如下

可以看下面的例子

model.setObjective(x + y) 

如果是比较复杂的目标函数，就可以使用创建表达式的全局函数quicksum()来完成。

比如下面的目标函数

$$\begin{array}{r}\min \sum _{i=1}^5{c}_i{x}_i\end{array}$$

具体代码如下

model = Model() 
c = [1, 2, 3, 4, 5]
x = {} 
for i in range(len(c)):
    x[i] = model.addVar(vtype="I", name="x" + str(i))
    
# set objective function by quicksum function 
model.setObjective(
        quicksum(c[i]*x[i] for i in range(len(c))),
        "minimize")   

创建完目标函数以后，可以查看目标函数的一些属性，比如

• 查看目标函数的sense，也就是$\max$还是$\min$
• model.getObjectiveSense()
• model.setMaximize： 直接设置模型的sense为$\max$

创建约束

创建约束的函数主要是：addCons()、addConsAnd()、addConsOr()、addConsXor() 等

创建一般约束

比如，想要加入下面的约束

$$\begin{array}{r}\sum _{i=1}^5{c}_i{x}_i⩽10\end{array}$$

则代码如下

model.addCons(quicksum(x[i] for i in range(len(c))) <= 10, "Width")  

更为复杂的情况，可以见下面的例子

min ⁡   ∑ i ∑ j c i j x i j ∑ i ∈ V x i j = 1,  ∀ j ∈ V , i ≠ j ∑ j ∈ V x i j = 1,  ∀ i ∈ V , i ≠ j μ i − μ j + N x i j ⩽ N − 1,  ∀ i , j ∈ V , i , j ≠ 0, i ≠ j x i j ∈ { 0, 1 } , μ i ⩾ 0,  μ i ∈ R 1 ∀ i ∈ V \begin{aligned} \min\text{ }\sum_i{\sum_j{c_{ij}x_{ij}}} \\ \sum_{i\in V}{x_{ij}}=\text{1, } & \hspace{1cm}\forall j\in V,i\ne j \\ \sum_{j\in V}{x_{ij}}=\text{1, } & \hspace{1cm}\forall i\in V,i\ne j \\ \mu _i-\mu _j+Nx_{ij}\leqslant N-\text{1, } & \hspace{1cm}\forall i,j\in V,i,j\ne \text{0,}i\ne j \\ x_{ij}\in \left\{ \text{0,}1 \right\} , \mu _i\geqslant \text{0, }\mu _i\in \mathbf{R}^1 &\hspace{1cm}\forall i\in V \end{aligned} min i∑​j∑​cij​xij​i∈V∑​xij​=1, j∈V∑​xij​=1, μi​−μj​+Nxij​⩽N−1, xij​∈{0,1},μi​⩾0, μi​∈R1​∀j∈V,i=j∀i∈V,i=j∀i,j∈V,i,j=0,i=j∀i∈V​

# 一些例子 
model = Model("atsp - mtz")

x,u = {},{}
for i in range(1,n+1):
    u[i] = model.addVar(lb=0, ub=n-1, vtype="C", name="u(%s)"%i)
    for j in range(1,n+1):
        if i != j:
            x[i,j] = model.addVar(vtype="B", name="x(%s,%s)"%(i,j))
                
for i in range(1,n+1):
    model.addCons(quicksum(x[i,j] for j in range(1,n+1) if j != i) == 1, "Out(%s)"%i)
    model.addCons(quicksum(x[j,i] for j in range(1,n+1) if j != i) == 1, "In(%s)"%i)

for i in range(1,n+1):
    for j in range(2,n+1):
        if i != j:
            model.addCons(u[i] - u[j] + (n-1)*x[i,j] <= n-2, "MTZ(%s,%s)"%(i,j))

创建广义约束

这里的广义约束就指逻辑约束等。比如And， Or等。
具体函数为

• model.addConsAnd([x,y,z], r): 表示， 当$x=y=z=1$时，$r=1$; 否则$r=0$；
• model.addConsOr([x,y,z], r) 表示， 当$x=1\text{ or }y=1\text{ or }z=1$时，$r=1$; 否则$r=0$；
• model.addConsXor([x,y,z], r):表示异或
  -异或
  

下面是SCIP给出的例子

##@file tutorial/logical.py
#@brief Tutorial example on how to use AND/OR/XOR constraints.
"""
N.B.: standard SCIP XOR constraint works differently from AND/OR by design.
The constraint is set with a boolean rhs instead of an integer resultant.
cf. http://listserv.zib.de/pipermail/scip/2018-May/003392.html
A workaround to get the resultant as variable is here proposed.
Public Domain, WTFNMFPL Public Licence
"""
from pyscipopt import Model
from pyscipopt import quicksum

def _init():
    model = Model()
    model.hideOutput()
    x = model.addVar("x","B")
    y = model.addVar("y","B")
    z = model.addVar("z","B")
    return model, x, y, z

def _optimize(name, m):
    m.optimize()
    print("* %s constraint *" % name)
    objSet = bool(m.getObjective().terms.keys())
    print("* Is objective set? %s" % objSet)
    if objSet:
        print("* Sense: %s" % m.getObjectiveSense())
    status = m.getStatus()
    print("* Model status: %s" % status)
    if status == 'optimal':
        for v in m.getVars():
            if v.name != "n":
                print("%s: %d" % (v, round(m.getVal(v))))
    else:
        print("* No variable is printed if model status is not optimal")
    print("")

def and_constraint(v=1, sense="minimize"):
    """ AND constraint """
    assert v in [0,1], "v must be 0 or 1 instead of %s" % v.__repr__()
    model, x, y, z = _init()
    r = model.addVar("r", "B")
    model.addConsAnd([x,y,z], r)
    model.addCons(x==v)
    model.setObjective(r, sense=sense)
    _optimize("AND", model)


def or_constraint(v=0, sense="maximize"):
    """ OR constraint"""
    assert v in [0,1], "v must be 0 or 1 instead of %s" % v.__repr__()
    model, x, y, z = _init()
    r = model.addVar("r", "B")
    model.addConsOr([x,y,z], r)
    model.addCons(x==v)
    model.setObjective(r, sense=sense)
    _optimize("OR", model)

def xors_constraint(v=1):
    """ XOR (r as boolean) standard constraint"""
    assert v in [0,1], "v must be 0 or 1 instead of %s" % v.__repr__()
    model, x, y, z = _init()
    r = True    
    model.addConsXor([x,y,z], r)
    model.addCons(x==v)
#     print('v=', v)  
    _optimize("Standard XOR (as boolean)", model)

def xorc_constraint(v=0, sense="maximize"):
    """ XOR (r as variable) custom constraint"""
    assert v in [0,1], "v must be 0 or 1 instead of %s" % v.__repr__()
    model, x, y, z = _init()
    r = model.addVar("r", "B")
    n = model.addVar("n", "I") # auxiliary
    model.addCons(r+quicksum([x,y,z]) == 2*n)
    model.addCons(x==v)
    model.setObjective(r, sense=sense)
    _optimize("Custom XOR (as variable)", model)

if __name__ == "__main__":
    and_constraint()  
    or_constraint()
    xors_constraint()
    xorc_constraint()

运行结果如下

* AND constraint *
* Is objective set? True
* Sense: minimize
* Model status: optimal
x: 1
y: 0
z: 0
r: 0

* OR constraint *
* Is objective set? True
* Sense: maximize
* Model status: optimal
x: 0
y: 1
z: 1
r: 1

* Standard XOR (as boolean) constraint *
* Is objective set? False
* Model status: optimal
x: 1
y: 1
z: 1

* Custom XOR (as variable) constraint *
* Is objective set? True
* Sense: maximize
* Model status: optimal
x: 0
y: 1
z: 0
r: 1

还有addGenConsIndicator()也是类似的用法，例如加入指示约束，

• 如果$y=1$，则$x+2y=1$

这样的约束就可以使用addGenConsIndicator()完成。

注意：SCIP似乎没有experssion的对象，也就是没有linear expression等的表达式对象，因此在添加约束的时候，也不能使用addTerms等的操作，只能通过quicksum拼凑表达式。

但是在Column Generation的时候，按列建模的时候，是可以使用addTerms的。

例如

        # add new column to the master problem
        col = Column()
        for i in range(m):
            if t[K][i] > 0:
                col.addTerms(t[K][i], orders[i])
        x[K] = master.addVar(obj=1, vtype="I", name="x(%s)"%K, column=col)

说实话，这一点有一些些不太方便。

求解模型

求解模型非常简单，直接model.optimize()就可以，但是有几个涉及到的相关函数需要注意一下

• model.hideOutput(quiet = False) # silent mode : 关闭log信息
• model.optimize():求解模型

model.hideOutput(quiet = False) # silent mode
model.optimize()

当然，求解之前，还可以对模型进行松弛等。

- relax = model.relax()：就是获得model的LP relaxation

获得解的信息

相关函数如下

• model.getObjVal(): 获得最优解对应的目标函数值
• model.getBestSol() ：获得最好解
• model.getVal(x) ：获得最优解中$x$的取值

例如

sol = model.getBestSol()
print("x: {}".format(sol[x]))
print("y: {}".format(sol[y]))  

运行结果

x: 0.0
y: 0.0

再比如

print("Optimal value:", model.getObjVal())
print((x.name, y.name), " = ", (model.getVal(x), model.getVal(y)))

运行结果

Optimal value: 0.0
('x', 'y')  =  (0.0, 0.0)

至此，Python(pyscipopt)调用SCIP求解最优化模型的一套流程就结束了，大家可以开开心心的调用了。

当然，很多小伙伴也许还会有些其他需求，见下文。

一些其他常用函数

• model.chgLhs():改变约束的左端项

model.chgLhs(Constraint cons, lhs)
-----------------------------
Docstring:
Change left hand side value of a constraint.

:param Constraint cons: linear or quadratic constraint
:param lhs: new left hand side (set to None for -infinity)
Type:      builtin_function_or_method

• model.chgRhs():改变约束的右端项

model.chgRhs(Constraint cons, rhs) 
-----------------------------
Docstring:
Change left hand side value of a constraint.

:param Constraint cons: linear or quadratic constraint
:param lhs: new left hand side (set to None for -infinity)
Type:      builtin_function_or_method

• model.chgVarLb():改变决策变量的下界

model.chgRhs(Variable var, lb) 
-----------------------------
Docstring:
Changes the lower bound of the specified variable.

:param Variable var: variable to change bound of
:param lb: new lower bound (set to None for -infinity)
Type:      builtin_function_or_method

• model.chgVarUb():改变决策变量的上界

model.chgVarUb(Variable var, ub) 
-----------------------------
Docstring:
Changes the upper bound of the specified variable.

:param Variable var: variable to change bound of
:param ub: new upper bound (set to None for +infinity)
Type:      builtin_function_or_method

• model.chgVarType():改变决策变量的类型

model.chgVarType(Variable var, vtype) 
-----------------------------
Docstring:
Changes the type of a variable

:param Variable var: variable to change type of
:param vtype: new variable type
Type:      builtin_function_or_method

• model.writeLP('Lp.lp'):将模型导出
• model.delCons(Constraint cons):删除约束
• model.delVar(var):删除一个决策变量
• model.getGap(): 获得当前模型的Gap, Docstring: Retrieve the gap, i.e. |(primalbound - dualbound)/min(|primalbound|,|dualbound|)|.
• model.getLhs(Constraint cons) ：获得约束的左端项
• model.getRhs(Constraint cons) ：获得约束的右端项
• model.getObjective() ：Retrieve objective function as Expr
• model.getObjVal() ：Retrieve the objective value of value of best solution.
• model.getStatus(): Retrieve solution status. 可能为optimal
• model.setParams(): Sets multiple parameters at once. :param params: dict mapping parameter names to their values.
• model.getNConss(): Retrieve the number of constraints.
• model.getNVars(): Retrieve number of variables in the problems
• model.getNVars(): Retrieve number of variables in the problems
• model.getDualMultiplier(Constraint cons): Retrieve the dual solution to a linear constraint.

这些函数，大概可以满足大部分的需求了。

今天就介绍到这里，下一期再见。

---

作者：刘兴禄，清华大学，清华-伯克利深圳学院博士在读

参考文献：
[1] Gerald Gamrath, Daniel Anderson, Ksenia Bestuzheva, Wei-Kun Chen, Leon Eifler, Maxime Gasse, Patrick Gemander, Ambros Gleixner, Leona Gottwald, Katrin Halbig, Gregor Hendel, Christopher Hojny, Thorsten Koch, Pierre Le Bodic, Stephen J. Maher, Frederic Matter, Matthias Miltenberger, Erik Mühmer, Benjamin Müller, Marc Pfetsch, Franziska Schlösser, Felipe Serrano, Yuji Shinano, Christine Tawfik, Stefan Vigerske, Fabian Wegscheider, Dieter Weninger, Jakob Witzig
Available at Optimization Online and as ZIB-Report 20-10, March 2020