Using pyscipopt in python I have set up a continuous variable and binary variable:
For j in J:
x[j] = model.addVar(vtype="C", name="x(%s)"%j)
y[j] = model.addVar(vtype="B", name="y(%s)"%j)
What I am trying to do is get a simple indicator variable of the form:
when x[j] = 0 -> y[j] = 0, and when x[j] > 0 -> y[j] = 1
But I cannot figure out how to set up the constraint. I tried to read up on Big M but still seem to be struggling.
And ultimately, the goal is to get the optimization model to use just one of x[3] or x[7], or neither, but not both in the solution.
So eventually once the indicator variable is working, I was thinking something along the lines of:
model.addCons(y[3] + y[7] <= 1)
Feel free to offer any suggestions on this too if there is a better way to do it.
just one of x[3] or x[7], or neither, but not both in the solution
A big-M approach would be:
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}
An approach based on indicator constraints would look like:
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}
In PySCIPOpt this can be specified using model.addConsIndicator. Note that this approach can be used when we don't have good upper bounds on x. You also may want to have a look at model.addConsCardinality for a simpler approach.
Coding
Here is the trivial part. Convert the mathematical formulation into code. The authors of the Python interface for SCIP have done their best to make this step very easy. So for the programmers who specialize in copy-paste: the big-M approach can be coded as:
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])))
The solution looks like:
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
Indeed not both x[3] and x[7] are nonzero.
The same result can be obtained by using 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])))
The results are:
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
Again not both x[3] and x[7] are non-zero.
PS. I am confused, why providing the code was needed. The detailed description in the answer should be more than enough to get you on the correct path. I prefer a carefully crafted and well-stated mathematical model above a bunch of code any time.