First you need to add z
as a binary variable:
z = m.addVars(I, J, vtype=GRB.BINARY, name="z")
Then you need constraints to ensure that z[i, j] = 1
if and only if c[i, j] <= 150
.
One way to do this is by using indicator constraints:
z = 1 -> c <= 150
z = 0 -> c >= 150
This is equivalent to
c > 150 -> z = 0
c < 150 -> z = 1
You add these as follows:
m.addConstrs((z[i, j] == 1) >> (c[i][j] <= 150) for i in I for j in J)
m.addConstrs((z[i, j] == 0) >> (c[i][j] >= 150) for i in I for j in J)
You can also model this yourself explicitly:
If you have upper and lower bounds M
and m
on the value of c[i][j] - 150
(i.e., M >= c[i][j] - 150 >= m
for all i, j
), you can use the following constraints:
M * (1-z) >= c - 150
m * z <= c - 150
If c > 150
, the right-hand-sides of both inequalities will be positive. The first one then forces 1 - z = 1
and hence z = 0
. The second inequality will trivially be satisfied.
If c < 150
, the right-hand-sides are negative. The first inequality becomes trivial while the second one forces z = 1
.
For M
the maximum entry in c
will do, for m
you can choose -150
if all c[i][j]
are non-negative.
You add these constraints as follows:
m.addConstrs( M * (1 - z[i, j]) >= c[i][j] - 150 for i in I for j in J )
m.addConstrs( m * z[i,j] <= c[i][j] - 150 for i in I for j in J )
Note that I neglected the case where c = 150
. This is because for floating point numbers equalities are always only considered to be satisfied within tolerances and hence there is no easy way to distinguish between strict and non-strict inequalities. You can approximate this with an epsilon, e.g.:
z = 0 -> c >= 150 + epsilon