PySCIPOpt 中的指示变量

Indicator Variable in PySCIPOpt

在python中使用pyscipopt 我设置了一个连续变量和二进制变量:

For j in J:
   x[j] = model.addVar(vtype="C", name="x(%s)"%j)   
   y[j] = model.addVar(vtype="B", name="y(%s)"%j)

我想要做的是获得一个简单的指标变量,其形式为:

when x[j] = 0 -> y[j] = 0, and when x[j] > 0 -> y[j] = 1

但我不知道如何设置约束。我试图阅读有关 Big M 的内容,但似乎仍在苦苦挣扎。

最终,目标是让优化模型在解决方案中仅使用 x[3] 或 x[7] 之一,或两者都不使用,但不能同时使用两者。

所以最终一旦指标变量起作用,我就在想一些类似的事情:

model.addCons(y[3] + y[7] <= 1)

如果有更好的方法,也欢迎提出任何建议。

just one of x[3] or x[7], or neither, but not both in the solution

big-M 方法 将是:

x[3] <= b[3]*U[3]      # b is a binary variable
x[7] <= b[7]*U[7]      # U is an upperbound on x (constant)
b[3]+b[7] <= 1
0 <= x[3] <= U[3]
0 <= x[7] <= U[7]
b[3],b[7] ∈ {0,1}

基于指标约束的方法如下所示:

b[3]==0 ==> x[3]=0   # or x[3]<=0 (this is used in addConsIndicator)
b[7]==0 ==> x[7]=0
b[3]+b[7] <= 1
b[3],b[7] ∈ {0,1}

在 PySCIPOpt 中,这可以使用 model.addConsIndicator 指定。请注意,当我们在 x 上没有很好的上限时,可以使用这种方法。您可能还想查看 model.addConsCardinality 以获得更简单的方法。

编码

这是琐碎的部分。将数学公式转换为代码。 SCIP 的 Python 接口的作者已尽最大努力使这一步变得非常简单。所以对于擅长copy-paste的程序员来说:big-M方式可以编码为:

from pyscipopt import Model

J = range(10)
U = [10*j for j in J]  # upper bound on x


model = Model("Formulation1")

#
# variables
#
x = [model.addVar(vtype="C", lb=0, ub=U[j]) for j in J]
b = [model.addVar(vtype="B") for j in J]

#
# objective
#
model.setObjective(sum(x), sense="maximize")

#
# constraints
#
model.addCons(x[3] <= b[3]*U[3])
model.addCons(x[7] <= b[7]*U[7])
model.addCons(b[3]+b[7]<=1)


#
# solve
#
model.optimize()

#
# print solution
#
for j in J:
    print("j=%2s, x[j]=%5s, b[j]=%5s" % (j,model.getVal(x[j]),model.getVal(b[j])))

解决方案如下:

j= 0, x[j]= -0.0, b[j]= -0.0
j= 1, x[j]= 10.0, b[j]= -0.0
j= 2, x[j]= 20.0, b[j]= -0.0
j= 3, x[j]=  0.0, b[j]=  0.0
j= 4, x[j]= 40.0, b[j]= -0.0
j= 5, x[j]= 50.0, b[j]= -0.0
j= 6, x[j]= 60.0, b[j]= -0.0
j= 7, x[j]= 70.0, b[j]=  1.0
j= 8, x[j]= 80.0, b[j]= -0.0
j= 9, x[j]= 90.0, b[j]= -0.0

确实不是 x[3] 和 x[7] 都非零。

使用model.addConsCardinality:

可以得到相同的结果
from pyscipopt import Model

J = range(10)
U = [10*j for j in J]  # upper bound on x


model = Model("Formulation3")

#
# variables
#
x = [model.addVar(vtype="C", lb=0, ub=U[j]) for j in J]

#
# objective
#
model.setObjective(sum(x), sense="maximize")

#
# constraints
#
# We use model.addConsCardinality
#    Add a cardinality constraint that allows at most 'cardval' many nonzero variables.
# Only the following two parameters are used:
#       :param consvars: list of variables to be included
#       :param cardval: nonnegative integer
#
model.addConsCardinality([x[3],x[7]],1)


#
# solve
#
model.optimize()

#
# print solution
#
for j in J:
    print("j=%2s, x[j]=%5s" % (j,model.getVal(x[j])))

结果是:

j= 0, x[j]= -0.0
j= 1, x[j]= 10.0
j= 2, x[j]= 20.0
j= 3, x[j]=  0.0
j= 4, x[j]= 40.0
j= 5, x[j]= 50.0
j= 6, x[j]= 60.0
j= 7, x[j]= 70.0
j= 8, x[j]= 80.0
j= 9, x[j]= 90.0

同样,并非 x[3] 和 x[7] 都非零。

PS。我很困惑,为什么需要提供代码。答案中的详细描述应该足以让您走上正确的道路。比起一堆代码,我更喜欢一个精心设计和表述良好的数学模型。