Constraining binary integer variables with indicator constraints

Specific Domains Optimization (Mathematical)
,
1

I want to constrain binary integer variables, x, with indicator constraints to express a change in a binary integer variable over a specific iterator, t.

I thought this could be nice since the different cases are obvious:

x[t-1] - x[t] = 0 - 0 = 0   # no change
x[t-1] - x[t] = 1 - 0 = 1   # shutting down
x[t-1] - x[t] = 0 - 1 = -1  # starting up
x[t-1] - x[t] = 1 - 1 = 0   # no change

However, it seems that there is something I don’t understand from the solver or the handling of the indicator constraints in JuMP. I wrote:

@constraint(m, x_up[n,t]  => {x[n,t-1] - x[n,t] == -1})

And I was returned a zero matrix.

If I negated the indicator (NB: the constraints are applied for every unit, n, in the system):

@constraint(m, !x_up[n,t]  => {x[n,t-1] - x[n,t] == -1})

I was returned something that would express x “going down” instead of “up”.

So, since I had something concrete other than a zero matrix, I was happy. I recalled that solvers like to work with tolerances and I then expressed my indicator constraints like:

@constraint(m, !x_up[n,t]  => {x[n,t-1] - x[n,t] >= 0-atol})
@constraint(m, !x_down[n,t] => {x[n,t-1] - x[n,t] <= 0+atol})

And it works. I get what I expected but it does not work without the negation (It does not work either when taking the tolerances into account). In other words

# This returns a zero matrix:
@constraint(m, x_up[n,t]  => {x[n,t-1] - x[n,t] <= 0+atol})    # up is 1 when -1
@constraint(m, x_down[n,t] => {x[n,t-1] - x[n,t] >= 0-atol})   # down is 1 when 1

# This returns the expected:
@constraint(m, !x_up[n,t]  => {x[n,t-1] - x[n,t] >= 0-atol})   # up is 0 when ">=0"
@constraint(m, !x_down[n,t] => {x[n,t-1] - x[n,t] <= 0+atol})  # down is 0 when "<=0"

Question: What is the reason behind this and is this expected behavior? The logic becomes more cryptic.

My operating code can be found in my github

NB: Cbc forced exited constraining binary integers with indicator constraints. CPLEX has no problems. I have not tried Gurobi yet.
This is the terminal when using Cbc (using atom as IDE):

command line - Cbc_C_Interface -solve -quit (default strategy 1)
Continuous objective value is 6.11352e+07 - 0.01 seconds
Cgl0003I 0 fixed, 648 tightened bounds, 203 strengthened rows, 1613 substitutions
Cgl0003I 0 fixed, 0 tightened bounds, 81 strengthened rows, 1642 substitutions
Cgl0003I 0 fixed, 0 tightened bounds, 47 strengthened rows, 1615 substitutions
Cgl0003I 0 fixed, 0 tightened bounds, 22 strengthened rows, 1615 substitutions
Cgl0003I 0 fixed, 0 tightened bounds, 7 strengthened rows, 1615 substitutions
Cgl0003I 0 fixed, 0 tightened bounds, 1 strengthened rows, 1615 substitutions
Cgl0003I 0 fixed, 0 tightened bounds, 0 strengthened rows, 1615 substitutions
Cgl0003I 0 fixed, 0 tightened bounds, 0 strengthened rows, 1615 substitutions
Cgl0003I 0 fixed, 0 tightened bounds, 0 strengthened rows, 1481 substitutions
Cgl0004I processed model has 6461 rows, 6348 columns (3560 integer (3560 of which binary)) and 19298 elements
Cbc0045I Fixing only non-zero variables.
Cbc0045I MIPStart solution provided values for 710 of 3560 integer variables, 65 variables are still fractional.
Cbc0038I Full problem 6461 rows 6348 columns, reduced to 1440 rows 1368 columns
Cbc0038I Probing was tried 0 times and created 0 cuts of which 0 were active after adding rounds of cuts (0.000 seconds)
Cbc0045I Mini branch and bound defined values for remaining variables in 0.68 seconds.
Cbc0045I MIPStart provided solution with cost 6.11679e+07

Julia has exited.
Press Enter to start a new session.
1 Like
2

Good question. I’ve just tested it: Your code works fine with Gurobi.

@constraint(m, x_up[n,t]  => {x[n,t-1] - x[n,t] == -1})

also works with Gurobi.

2 Likes
3

Aha! Interesting. I need to finalize the download procedure for Gurobi then! Thanks for testing.

1 Like
4

Maybe some solvers don’t allow equality constraints as part of indicators (e.g. SCIP). But that should raise an error during model building, not from the solver.

1 Like
5

This worked for me with Cbc.

using JuMP, Cbc
model = Model(Cbc.Optimizer)
set_silent(model)
@variable(model, z, Bin)
@variable(model, x[1:2], Bin)
@constraint(model, z => {x[1] - x[2] == -1})
@objective(model, Max, z)
optimize!(model)
@show value.(x) == [0.0, 1.0]

model = Model(Cbc.Optimizer)
set_silent(model)
@variable(model, z, Bin)
@variable(model, x[1:2], Bin)
@constraint(model, !z => {x[1] - x[2] == -1})
@objective(model, Min, z)
optimize!(model)
@show value.(x) == [0.0, 1.0]
1 Like
6

Simplicity is good! Thanks for boiling the trouble down.

I see now that my question is regarding some of the underlying inference that ought to be understood from indicator constraints:

This inference works as it should while having z decision variable in the objective function like in your cases. I.e. as from your case:

z = 1 if obj. func. is to maximize z, and thus is x[1] = 0 and x[2] = 1
z = 0 if obj. func. is to minimize z, and thus is x[1] = 0 and x[2] = 0

But by adding some terms, where the intuition is to force a value in x[1]=0 and x[2]=1 so that z=1, then the inference is gone so that z actually is returned as 0.

Added is a variable y[1:2] with the constraints that it will require sum(y)==2 and a constraint that will define the feasible region for y with the help of values stored in d. These values forces what we want. The objective function is to minimize sum(y)+z:

using JuMP, Cbc
d = [0.9 0.9; 1.8 2.1]
model = Model(Cbc.Optimizer)
set_silent(model)
@variable(model, z, Bin)
@variable(model, x[1:2], Bin)
@variable(model, y[1:2])            
@constraint(model, sum(y[:]) == 2)                 
for i in 1:2                      
    @constraint(model, d[i,1]*x[i] <= y[i])
    @constraint(model, y[i] <= d[i,2]*x[i])
end
@constraint(model, z => {x[1] - x[2] == -1})
@objective(model, Min, sum(y)+z)                   
optimize!(model)
@show value.(x) == [0.0, 1.0]                       # true
@show value.(x)                                     # [0.0, 1.0]
@show value.(y)                                     # [0.0, 2.0]
@show value.(z)                                     # 0.0

And here with the cryptic logical expression but which works and keeps the inference:

using JuMP, Cbc
d = [0.9 0.9; 1.8 2.1]
model = Model(Cbc.Optimizer)
set_silent(model)
@variable(model, z, Bin)
@variable(model, x[1:2], Bin)
@variable(model, y[1:2])            
@constraint(model, sum(y[:]) == 2)                 
for i in 1:2                      
    @constraint(model, d[i,1]*x[i] <= y[i])
    @constraint(model, y[i] <= d[i,2]*x[i])
end
@constraint(model, !z => {x[1] - x[2] >= 0-10e-6})
@objective(model, Min, sum(y)+z)                   
optimize!(model)
@show value.(x) == [0.0, 1.0]                       # true
@show value.(x)                                     # [0.0, 1.0]
@show value.(y)                                     # [0.0, 2.0]
@show value.(z)                                     # 1.0

And just for a common sanity check. Do we agree that:

!z => {x[1] - x[2] >= 0-10e-6}    # z=0 if 0 - 1 = -1 >= 0, and z=0 <=> !z=1

# Is equal to:
z => {x[1] - x[2] < 0-10e-6}      # z=1 if 0 - 1 = -1 < 0

# Which is similar to (disregarding tolerances and assuming integers)
z => {x[1] - x[2] == -1}          # z=1 if 0 - 1 = -1 == -1
7

Your confusion is likely that these are one-way implications, not iff <=>.

If LHS, then the RHS holds. But the reverse is not true. If the RHS holds, we can say nothing about the LHS.

3 Likes
8

That might be my confusion :slight_smile: Is there a neat way to implement such a relation?

If not, I guess I cannot turn my constraint around and demand that z=1 if (!x[1] && x[2]) or similar:

@constraint(model, (!x[1] && x[2])  => {z == 1})

I hope you see what I am trying to achieve.

Even though there is no bi-implication intended in the indicator constraint, the behavior of the rather obscure logic seems to work quite well (despite it is horrible to read).
Here is a plot from my code (linked to in the main post). In the middle row, you see the “unit activation” values, where red and green respectively represent the unit being turned off or on. The dark green points indicate the time where the unit is turned on.

image
image900×600 159 KB

I would of course really want to refrain from using “inadvisable” methods but I am not quite sure how to achieve my goal in another way.

9

I guess I cannot turn my constraint around

No.

The trick is to find constraint(s) that have an equivalent truth-table. For example
z <=> x[1] - x[2] == -1, where x, z are binary, is equivalent to

@constraint(model, x[1] - x[2] - 2 * (1 - z) <= -1)

and (!x[1] && x[2]) => {z == 1} is equivalent to

@constraint(model, z >= -x[1] + x[2])

The MOSEK Modeling Cookbook (9 Mixed integer optimization — MOSEK Modeling Cookbook 3.3.0) has a number of helpful examples.

2 Likes
10

Thank you very much for your elaboration. The reference to the modeling cookbook is very useful. It will hopefully decrease the number of new topics on basics from my side :sweat_smile:

I updated the equivalent stackoverflow with the details from your answer above.

Many thanks once again and all the best!

1 Like