PySCIPOpt中的指示符变量

Question

在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

但我不知道如何设置约束。我试着读大M，但似乎仍在挣扎

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

因此，一旦指标变量起作用，我最终会想到以下几点：

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

如果有更好的方法，也可以提供任何建议

Answer 1

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

大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}

在pyscipot中，可以使用model.addConsIndicator来指定这一点。注意，当x没有好的上界时，可以使用这种方法。您可能还想看看model.addConsCardinality以获得更简单的方法

编码

这里是琐碎的部分。将数学公式转换成代码。SCIP Python接口的作者已经尽了最大努力使这一步骤变得非常简单。因此，对于专门从事复制粘贴的程序员来说：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]都是非零的

我很困惑，为什么需要提供代码。答案中的详细描述应该足以让您走上正确的道路。在任何时候，我都喜欢在一堆代码之上建立一个精心制作、陈述良好的数学模型