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[1] = Linear_exp[1]
# 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[1] = aux_M[1]
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[1] = Linear_exp[1]
aux_Affine_exp[1] = 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)
