Indicator modeling of TempConstr of MLinExpr
Answered
I want to realize such a problem, whose objective function depends on which polytope the variable X locates in. X is a 1-D MVar, A and b are 2-D and 1-D NumPy.array respectively. The obj1, obj2, and obj3 are MLinExpr expressed by X. Those three constraints construct non-overlapping polytopes. (Actually the problem is derived from multiparametric programming, where these polytopes are critical regions.)
And I want to use indicator constraints to model it. So I add another binary var z to indicate which MLinExpr constraint is satisfied. But when I use the following small example, I found that I can not use a single binary to indicate a whole constraint expressed by matrix multiplication, i.e. if z==1, all rows in AX<=b are satisfied.
from gurobipy import *
import numpy as np
M=Model()
X = M.addMVar(12)
A = np.random.rand(5,12)
z = M.addVar(vtype=GRB.BINARY)
M.addConstr((z==1) >> (A@X<=0))
The above code raises ValueError: can only convert an array of size 1 to a Python scalar
But I can use groups of z to indicate each row of a MLinExpr constraint. If A is m*n, e.g. (z[i]==1) >> ((A@X)[i] <=0). Until all of these rows are indicated, I still need to use another binary variable y to indicate which objective function to choose. e.g.
(y[1]==1) >> ((z[1]+...z[m])==m)
(y[1]==1) >> (obj == obj1)
I don't know whether the above description is correct. I think if I need to model like this, the number of binary variables will be too large. Can someone help me to model the problem I post at the top? Thanks!!!!!!
-
Hi Jiang,
I believe you need the following code to model the feasible set as the union of non-overlapping polytopes:
import gurobipy as gp
from gurobipy import GRB
import numpy as np
# Some random data defining three polygons
A1 = np.random.rand(5,12).round(2)
b1 = np.random.rand(5).round(2)
A2 = np.random.rand(6,12).round(2)
b2 = np.random.rand(6).round(2)
A3 = np.random.rand(3,12).round(2)
b3 = np.random.rand(3).round(2)
model = gp.Model()
# Spatial variables
X = model.addMVar(12, name="x")
# Polytope selection variables
z = model.addVars([1, 2, 3], vtype=GRB.BINARY, name="z")
# Indicator constraints must be added one row at a time
for i in range(A1.shape[0]):
model.addConstr((z[1] == 1) >> (A1[i, :] @ X <= b1[i]))
for i in range(A2.shape[0]):
model.addConstr((z[2] == 1) >> (A2[i, :] @ X <= b2[i]))
for i in range(A3.shape[0]):
model.addConstr((z[3] == 1) >> (A3[i, :] @ X <= b3[i]))
# Enforce that exactly one polytope constraint set is satisfied
model.addConstr(z.sum() == 1)In this form, the binary variables z[1] ... z[3] select which polytope the point X lies in. z[1] must be applied to the entire polygon A1 @ X <= b1, but as you pointed out, the code:
m.addConstr((z[1] == 1) >> (A1 @ X <= b1)) # raises ValueError
does not work. This is because the indicator constraint syntax is not integrated with the matrix-oriented API (as of Gurobi v10). As a result, you can only add a single constraint at a time by iterating through the rows of A1 and entries of b1. This is done by slicing the matrix A1 as I have shown in the example above.
Defining the objective function
With the feasible region described using the model above, you can then define your objective using the same variables z and indicator constraints, similar to what you described:
obj = model.addVar()
model.addConstr((z[1] == 1) >> (obj == obj1))
model.addConstr((z[2] == 1) >> (obj == obj2))
...I would also note that if obj1, obj2, etc are constants and not some function of other variables, then you can define your objective directly over the binary variables z, without additional constraints. In that case, the objective would be simply written as:
obj1 * z[1] + obj2 * z[2] + ...
1 -
Hi Simon:
Thank you so much for this concrete answer!
Still, there is something unclear for me about your code. That is,
# Indicator constraints must be added one row at a time
for i in range(A1.shape[0]):
model.addConstr((z[1] == 1) >> (A1[i, :] @ X <= b1[i]))I understand what you want to model is using z to indicate which polytope the variable X is in. But maybe you want to model like only all rows of matrix A1@X are satisfied, then z[1] should be 1. In your code, I think there may exist such a condition that: A has 5 rows, in this loop, if A[1,:]@X, A[2,:]@X, A[3,:]@X is smaller than b1[1], b1[2], b1[3] respectively, then z[1]==1 until now. But if A[3,:]@X >= b1[3], I think z[1] will absolute turn into 0. The problem is, in the last iteration, if A[5,:]@X<= 0, then z[1] will turn back to 1 again. So at last, z[1]==1 but not all rows in A are satisfied the constraints. Is it true?
Actually, to deal with this problem, I introduce a new variable y to indicate each row of A. For example,
# Use y to indicate each row constraint
y = model.addVars([0,1,2,3,4], vtype=GRB.BINARY)
for i in range(A1.shape[0]):
model.addConstr((y[i] == 1) >> (A1[i, :] @ X <= b1[i]))
# Use z to indicate whether all rows constraints are satisfied (if z==1, X in this polytope)
model.addConstr(z[1] == and_([y[0],y[1],y[2],y[3],y[4]]))
model.addConstr((z[1]==1) >> (obj = obj1))I don't know whether my code is correct. And I am glad if you can answer my confusion about your code. Thanks a lot !!!!!
0 -
In your code, I think there may exist such a condition that: A has 5 rows, in this loop, if A[1,:]@X, A[2,:]@X, A[3,:]@X is smaller than b1[1], b1[2], b1[3] respectively, then z[1]==1 until now. But if A[3,:]@X >= b1[3], I think z[1] will absolute turn into 0. The problem is, in the last iteration, if A[5,:]@X<= 0, then z[1] will turn back to 1 again. So at last, z[1]==1 but not all rows in A are satisfied the constraints. Is it true?
No, this is not correct. There are two key points here:
1. z[1] can only take one value, and
2. the indicator constraint only applies to the z[1] = 1 condition, i.e. if z[1] = 1, then A1[i, :] @ X <= b1[i] must be true, but a feasible solution can have z[1] = 0 and A1[i, :] @ X <= b1[i].Thus, in the situation you describe (where A1[3,:]@X >= b1[3]) z[1] must be 0. This is the correct result, since if A1[3,:]@X >= b1[3], then X lies outside the polytope defined by A1X <= b1.
Your code produces the same result in terms of the feasible set of X, but the y variables you have introduced are redundant.
0
Please sign in to leave a comment.
Comments
3 comments