• Gurobi Staff

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

# 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

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:

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] + ...

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
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)

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 !!!!!

• Gurobi Staff

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.