# Updating expressions for indicator constraints in JuMP

**URL:** <https://discourse.julialang.org/t/updating-expressions-for-indicator-constraints-in-jump/59059>\
**Category:** Optimization (Mathematical)\
**Tags:** jump, optimization\
**Created:** [April 11, 2021, 10:29pm UTC](https://discourse.julialang.org/t/updating-expressions-for-indicator-constraints-in-jump/59059 "2021-04-11T22:29:26Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![Maria95](https://avatars.discourse-cdn.com/v4/letter/m/e47c2d/32.png) [@Maria95](https://discourse.julialang.org/u/Maria95)\
**Post date:** [April 11, 2021, 10:29pm UTC](https://discourse.julialang.org/t/updating-expressions-for-indicator-constraints-in-jump/59059/1 "2021-04-11T22:29:26Z")

</div>

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`

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

 ![JuMP equation](https://global.discourse-cdn.com/julialang/original/3X/4/1/41219e9ce4ef34fdea3f28684d086d779f6ac072.png)

My MWE:

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

```

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [April 11, 2021, 11:53pm UTC](https://discourse.julialang.org/t/updating-expressions-for-indicator-constraints-in-jump/59059/2 "2021-04-11T23:53:48Z")

</div>

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.3.0](https://docs.mosek.com/modeling-cookbook/mio.html#exact-absolute-value)

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

```julia
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[1], 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[3] == 0)

```

---

<div class="post-metadata">

**Author:** ![Maria95](https://avatars.discourse-cdn.com/v4/letter/m/e47c2d/32.png) [@Maria95](https://discourse.julialang.org/u/Maria95)\
**Post date:** [April 13, 2021, 12:55am UTC](https://discourse.julialang.org/t/updating-expressions-for-indicator-constraints-in-jump/59059/3 "2021-04-13T00:55:38Z")

</div>

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:  
 ![expressions](https://global.discourse-cdn.com/julialang/original/3X/6/3/63f1c84b3f3a8df52e8ed6fce3f478672aa6c5a7.png)  
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?

> [@odow](#):
>
> ```julia
> A = rand(3, 2)
> 
> ```

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

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [April 13, 2021, 1:02am UTC](https://discourse.julialang.org/t/updating-expressions-for-indicator-constraints-in-jump/59059/4 "2021-04-13T01:02:59Z")

</div>

Don’t you just want

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

```

---

<div class="post-metadata">

**Author:** ![Maria95](https://avatars.discourse-cdn.com/v4/letter/m/e47c2d/32.png) [@Maria95](https://discourse.julialang.org/u/Maria95)\
**Post date:** [April 13, 2021, 1:10am UTC](https://discourse.julialang.org/t/updating-expressions-for-indicator-constraints-in-jump/59059/5 "2021-04-13T01:10:34Z")

</div>

Thanks again and also for sending Mosek Modelling cookbook.

> [@odow](#):
>
> ```julia
> y[i] == (1 + Factor1[i]) * y_p[i] + (1 + Factor2[i]) * y_n[i]
> 
> ```

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 😃

---

<div class="post-metadata">

**Author:** ![Maria95](https://avatars.discourse-cdn.com/v4/letter/m/e47c2d/32.png) [@Maria95](https://discourse.julialang.org/u/Maria95)\
**Post date:** [April 13, 2021, 6:52pm UTC](https://discourse.julialang.org/t/updating-expressions-for-indicator-constraints-in-jump/59059/6 "2021-04-13T18:52:08Z")

</div>

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:

```julia
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[1], 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[1] = A[:,1]'*x_r - b[1]
	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

```

---

<div class="post-metadata">

**Author:** ![odow](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/odow/32/28685_2.png) [@odow](https://discourse.julialang.org/u/odow)\
**Post date:** [April 13, 2021, 8:15pm UTC](https://discourse.julialang.org/t/updating-expressions-for-indicator-constraints-in-jump/59059/7 "2021-04-13T20:15:55Z")

</div>

I made a typo. It should be

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

```

---

<div class="post-metadata">

**Author:** ![Maria95](https://avatars.discourse-cdn.com/v4/letter/m/e47c2d/32.png) [@Maria95](https://discourse.julialang.org/u/Maria95)\
**Post date:** [April 14, 2021, 9:28pm UTC](https://discourse.julialang.org/t/updating-expressions-for-indicator-constraints-in-jump/59059/8 "2021-04-14T21:28:20Z")

</div>

Thank you, much appreciated