 # Updating expressions for indicator constraints in JuMP

Hello, I am trying to implement a constraint with a conditional statement that I am making piecewise affine by introducing binary variables and using the big M method. I know that JuMP supports indicator constraints but I have not been able to find any relevant examples for more complicated use cases such as this one. My question may seem long, but I wanted to give enough information so that it is clear.

When I run my code and look at how the expressions are being built, it does not look like that it is doing what it should be doing. I will really appreciate if some can guide me on:
What changes do I need to make to the code? If indicator variables can be used it will make the code more readable. Could this constraint be expressed using indicator constraints in JuMP?

The constraint requires selecting between two different parameters, say factor 1 and factor 2, depending upon whether the value is positive or negative at the previous period:
The constraint requires defining an expression, let’s call it Linear_exp
`Linear_exp = @expression(model, A*x-b)` (where A and b are constants, and x are variables)
Now define another expression: `Affine_exp`

``````Affine_exp = Linear_exp
# Without (piecewise) linearizing, the constraint will look like:
for k=2:length(Affine_exp)
Affine_exp[k] = Linear_exp[i] + (Affine_exp[k-1] > 0 Affine_exp[k-1]*factor1[i] : Affine_exp[k-1] *factor2[i])
end
@constraint(model, zero_balance, Affine_exp[end] == 0)
``````

(I don’t think I can write latex here so I am pasting an image with some equations)

My MWE:

``````using JuMP, Cbc
function some_data()
cost = [500,400,300,600,450]
A= [100 120 130 140 150; 100 150 200 300 90; 100 200 150 140 150; 100 120 130 140 150; 200 300 400 500 200; 0 0 0 0 0; 0 0 0 0 0; 0 0 0 0 0]
b = [200, 300, 400, 1000, 500, 20, 30, 10]
factor1 = [0.001, 0.002, 0.003, 0.004, 0.0045, 0.005,0.0055, 0.006]
factor2 = [0.0015, 0.0025, 0.0035, 0.0048,0.0057, 0.007, 0.0065,0.0075]
t = length(b)
n = size(A,2)

return cost, A, b, factor1, factor2,n , t
end

function bigM(A,b,factor1,factor2, n,t)
# Calculate suitable values of big M
aux_M = A*ones(n) - b
M = Vector{Float64}(undef,length(aux_M))
M = aux_M
for i  = 2:length(M)
M[i] = aux_M[i] + (M[i-1] > 0 ? M[i-1]*factor1[i] : M[i-1] * factor2[i])
end
M = round.(abs.(M)) .+ 100
return M
end

function indicator_model(cost,A,b,factor1,factor2,n,t)
model = Model(Cbc.Optimizer)
@variable(model,0<= x[i=1:n] <= 1)
@variable(model,y[i=1:t], Bin)
@objective(model, Min, sum(cost[i]*x[i] for i in 1:n))

Linear_exp = @expression(model, A*x-b)
Affine_exp =[AffExpr(0.0) for i = 1:length(Linear_exp)]
aux_Affine_exp = [AffExpr(0.0) for i = 1:length(Linear_exp)]
Affine_exp = Linear_exp
aux_Affine_exp = 0.0

for k = 2:length(Linear_exp)
@constraint(model, Affine_exp[k-1] <= (1-y[k]) *M[k])
@constraint(model,  0 <= aux_Affine_exp[k] - factor1[k] * Affine_exp[k-1] + (1-y[k]) *M[k] - (1-y[k]) *M[k] <= 0)
@constraint(model, Affine_exp[k-1] >= -y[k] *M[k])
@constraint(model, -y[k] *M[k] <= aux_Affine_exp[k] - factor2[k])
@constraint(model, 0 <= aux_Affine_exp[k] - factor2[k] * Affine_exp[k-1] + (1-y[k]) *M[k] - y[k] *M[k] <= 0)
Affine_exp[k] = Linear_exp[k]+aux_Affine_exp[k]
end

@constraint(model,con_zero, Affine_exp[t] == 0.0)

return x, model

end

cost, A, b, factor1, factor2,n , t = some_data()
M = bigM(A,b,factor1,factor2, n,t)

x, model = indicator_model(cost,A,b,factor1,factor2,n,t)
``````

You’re close, but your big-M reformulation isn’t quite right.

Mosek has a good section on this trick: 9 Mixed integer optimization — MOSEK Modeling Cookbook 3.2.2

Here’s something to point you in the right direction (I haven’t tested it; there might be typos):

``````A = rand(3, 2)
b = rand(3)
Factor1 = rand(3)
Factor2 = rand(3)

model = Model()
@variables(model, begin
x[1:2]
y[1:3]
y_p[1:3] >= 0
y_n[1:3] >= 0
z[1:3], Bin
end)
@expression(model, expr[i=1:3], A[:, i]' * x - b[i] + y[i])
fix(y, 0.0)
for i in 2:3
@constraints(model, begin
y_p[i] - y_n[i] == expr[i - 1]
y[i] == Factor1[i] * y_p[i] + Factor2[i] * y_n[i]
y_p[i] <= M * z[i]
y_n[i] <= M * (1 - z[i])
end)
end
@constraint(model, expr == 0)
``````

Thank you @odow much appreciated. You are right my big M formulation is not correct.

I would like to check few points with you, I hope that is ok

The code with your sample data works, but anything else I try the solver returns infeasible problem (even the problems I know are feasible). This could be because I should have defined the expression: If I change this line of code: `y[i] == Factor1[i] * y_p[i] + Factor2[i] * y_n[i]` to
`y[i] == Factor1[i] * y_p[i] + Factor2[i] * y_n[i] + y[i-1]` will this effectively change this constraint? I must be missing another line of code that needs to be changed as well to implement this change.

In your code, you have used one big M rather than 2. Is that deliberate or it is just for simplicity?

I think you intended `A = rand(2,3)`

Don’t you just want

``````y[i] == (1 + Factor1[i]) * y_p[i] + (1 + Factor2[i]) * y_n[i]
``````

Thanks again and also for sending Mosek Modelling cookbook.

I also thought about this one but for some reason decided to add `y[i-1]`
I will spend some time to understand why the solver is returning the problem infeasible for different data. I will mark this as solved but I may have to come back with questions later Hi @odow I have tried your suggestion, but it doesn’t seem to update the constraints correctly:
` y[i] == (1 + Factor1[i]) * y_p[i] + (1 + Factor2[i]) * y_n[i]`

However, if I change this to: `y[i] == Factor1[i-1] * y_p[i] + Factor2[i-1] * y_n[i] + expr[i-1]` then I do get feasible result sometimes if the solver doesn’t take for ever.

I do not understand why there should be any difference between these two lines of code. Could this be a bug in JuMP or there is something wrong in the formulation?

Cbc is taking forever when I give it a larger problem, so I tried Bonmin and Couenne but they both finish with a message “The LP is infeasible or too expensive”.

A full working example is here:

``````using JuMP, Cbc

function some_data()
cost = [500,400,300,600,450]
A= [100 120 130 140 150; 100 150 200 300 90; 100 200 150 140 150; 100 120 130 140 150; 200 300 400 500 200]
b = [200, 300, 400, 1000, 500, 20, 30, 3]
# Append matrix with zeros to ensure dimension matches with b for any operations on the expression
if size(A,2) < length(b)
A = hcat(A,zeros(size(A,1),length(b)-size(A,2)))
elseif length(b) < size(A,2)
b = vcat(b,zeros(size(A,2)-length(b)))
end
Factor1 = [0.001, 0.002, 0.003, 0.004, 0.0045, 0.005,0.0055, 0.006]
Factor2 = [0.0015, 0.0025, 0.0035, 0.0048,0.0057, 0.007, 0.0065,0.0075]
t = length(b)
n = size(A,1)

return cost, A, b, Factor1, Factor2, n , t
end

function indicator_cons(A,b,Factor1,Factor2,M, n, t)
model = Model(Cbc.Optimizer)

@variables(model, begin
0 <= x[1:n] <= 1
y[1:t]
y_p[1:t] >= 0
y_n[1:t] >= 0
z[1:t], Bin
end)

@expression(model, expr[i=1:t], A[:, i]' * x - b[i] + y[i])
fix(y, 0.0)

for i in 2:t
@constraints(model, begin
y_p[i] - y_n[i] == expr[i - 1]
y[i] == (1+Factor1[i-1]) * y_p[i] + (1+Factor2[i-1]) * y_n[i]
# change the previous line to
#y[i] == Factor1[i-1] * y_p[i] + Factor2[i-1] * y_n[i] + expr[i-1]
y_p[i] <= M * z[i]
y_n[i] <= M * (1 - z[i])
end)
end
@constraint(model, con_zero, expr[t] == 0)

return x, y, expr, con_zero, model
end

function audit(x_r,A,b,Factor1,Factor2,t)
# Function to test results
lastvalue = Vector{Float64}(undef,length(b))

lastvalue = A[:,1]'*x_r - b
for i in 2:t
if lastvalue[i-1] > 0
lastvalue[i] = A[:,i]'*x_r - b[i] +  lastvalue[i-1]*(1+Factor1[i-1])
else
lastvalue[i] = A[:,i]'*x_r - b[i] +  lastvalue[i-1]*(1+Factor2[i-1])
end
end
return lastvalue[t]
end

# Test code
cost, A, b, Factor1, Factor2,n , t = some_data()

M=1000000
x, y, expr, con_zero, model = indicator_cons(A,b,Factor1,Factor2,M, n, t)
optimize!(model)

x_r = value.(x)
audit(x_r,A,b,Factor1,Factor2,t)
# This return -377.71088024605297
# Retry by changing to y[i] == Factor1[i-1] * y_p[i] + Factor2[i-1] * y_n[i] + expr[i-1]
# audit(x_r,A,b,Factor1,Factor2,t)  returns -1.3500311979441904e-13
``````

I made a typo. It should be

``````y[i] == (1 + Factor1[i]) * y_p[i] - (1 + Factor2[i]) * y_n[I]
^
-, not +
``````
1 Like

Thank you, much appreciated