# Cut generation not working lazy constraints - lagrangian relaxation

**URL:** <https://discourse.julialang.org/t/cut-generation-not-working-lazy-constraints-lagrangian-relaxation/9948>\
**Category:** Optimization (Mathematical)\
**Created:** [March 24, 2018, 6:24pm UTC](https://discourse.julialang.org/t/cut-generation-not-working-lazy-constraints-lagrangian-relaxation/9948 "2018-03-24T18:24:40Z")\
**Posts on this page:** 7\
**Page:** 1

<div class="post-metadata">

**Author:** ![tlagos](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tlagos/32/4490_2.png) [@tlagos](https://discourse.julialang.org/u/tlagos)\
**Post date:** [March 24, 2018, 6:24pm UTC](https://discourse.julialang.org/t/cut-generation-not-working-lazy-constraints-lagrangian-relaxation/9948/1 "2018-03-24T18:24:40Z")

</div>

Hello,

I have the following implementation to generate cuts for the Lagrangian relaxation problem for the shortest path with time constraints:

# \max\_{\lambda\ge0} { L( \lambda ) }

# X is a set of extreme points

```
function solveLag(c,t,T,X)
  MLR,λ,z=masterLgr(transpose(c)*X,transpose(t)*X,T)

  function cut(cb)
     println("----\nInside cut callback")
     masterval=getvalue(z);
     𝝺=getvalue(λ)
     col,val=ColGen(c,t,𝝺,0,E,V,s,r) #Generates an extreme point (path) 
     X=hcat(X,col)
     val=val - 𝝺*T
     # println("col=",col)
     if masterval+0.0001<=val
       return
     end
     @lazyconstraint(cb, dot(c,col) + λ*(dot(t,col) - T) >= z)
  end

  addlazycallback(MLR, cut)
  solve(MLR)
  obj=getobjectivevalue(MLR)
  𝝺=getvalue(λ)
  return obj,𝝺,X
end

```

Doesn’t enter in the function cut and ends with a solution that has masterval+0.0001\>val.

---

<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:** [March 25, 2018, 5:14am UTC](https://discourse.julialang.org/t/cut-generation-not-working-lazy-constraints-lagrangian-relaxation/9948/2 "2018-03-25T05:14:33Z")

</div>

It’s a bit hard to figure out what is going on without a [MWE](https://stackoverflow.com/help/mcve).

Here is an example using lazy constraints  
[https://github.com/JuliaOpt/JuMP.jl/blob/v0.18.0/examples/simplelazy.jl](https://github.com/JuliaOpt/JuMP.jl/blob/v0.18.0/examples/simplelazy.jl)

---

<div class="post-metadata">

**Author:** ![tlagos](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tlagos/32/4490_2.png) [@tlagos](https://discourse.julialang.org/u/tlagos)\
**Post date:** [March 27, 2018, 12:52pm UTC](https://discourse.julialang.org/t/cut-generation-not-working-lazy-constraints-lagrangian-relaxation/9948/3 "2018-03-27T12:52:55Z")

</div>

I think I was able to reproduce the same kind of problem in this example:

```
using JuMP
using Gurobi
slvr=GurobiSolver(Method=0)

# \min \{y_1 + 4y_2 + 2x_1 + 3x_2 \}
# s.t. y_1 - 3y_2 + x_1 -2x_2 \leq -2
# y_1 + 3y_2 + x_1 + x_2 \geq 3
# x_1,x_2\geq 0 , y_1,y_2 positive integers
#Initially we choose y_1^1 = y_2^1 = 0 , z^1=-\infty

function definemaster(d,D,b,λ)
    PM = Model(solver=slvr)
    @variable(PM, z)
    @variable(PM, y[i in 1:length(d)] >= 0, Int)
    λD= [sum(D[j,i]*λ[j] for j=1:length(b)) for i=1:length(d)]
    @constraint(PM,Ep[k in 1:1],z>= sum(y[i]*(d[i]-λD[i]) for i =1:length(d)) + sum(λ[j]*b[j] for j =1:length(b)) )
    # println(Ep[1])
    @objective(PM, :Min, z)
    return PM,y,z
end

function solvesubproblem(b,Dy,A,c)
    sub=Model(solver=slvr)
    @variable(sub,λ[j in 1:length(b)]>=0)
    @constraint(sub,x[i =1:length(c)], sum( A[j,i] * λ[j] for j =1:length(b) ) <= c[i])
    @objective(sub,:Max,sum( (b[j]-Dy[j])*λ[j] for j =1:length(b) ))
    status=solve(sub)
    if status==:Optimal
        return :Optimal,getvalue(λ),getobjectivevalue(sub),getdual(x)
    elseif status==:Unbounded
        return :Unbounded,getvalue(λ)
    else
        println("I'm in a really bad shape, status is ",status)
        return status
    end
end

function solvebenders(c,d,A,D,b)
    out=solvesubproblem(b,zeros(b),A,c)
    PM,y,z=definemaster(d,D,b,out[2])
    function cut(cb)
        println("----\nInside cut callback")
        masterval=getvalue(z);
        println("Current objective value of master problem is: ", masterval)
        𝐲=getvalue(y)
        println("Solution for 𝐲 is: ", 𝐲)
        Dy=[sum(D[j,i]*𝐲[i] for i=1:length(d)) for j=1:length(b)]
        out=solvesubproblem(b,Dy,A,c)
        λ=out[2]
        λD= [sum(D[j,i]*λ[j] for j=1:length(b)) for i=1:length(d)]
        expr=zero(AffExpr)
        if out[1]==:Unbounded
            expr+= sum( b[j]*λ[j] for j =1:length(b) )
            for j in 1:length(d)
                push!(expr,-λD[j],y[j])
            end
        elseif out[1]==:Optimal
            if -0.0001+ sum( b[j]*λ[j] for j =1:length(b)) + sum( (d[j]-λD[j])*𝐲[j] for j = 1:length(d) )<=masterval<=0.0001+ sum( b[j]*λ[j] for j =1:length(b)) + sum( (d[j]-λD[j])*𝐲[j] for j = 1:length(d) )
                println("masterval=",masterval)
                println( sum( b[j]*λ[j] for j =1:length(b)) + sum( (d[j]-λD[j])*𝐲[j] for j = 1:length(d)))
                return
            end
            expr+= sum( b[j]*λ[j] for j =1:length(b) )
            for j in 1:length(d)
                push!(expr,d[j]-λD[j],y[j])
            end
            push!(expr,-1.0,z)
        else
            println("please break")
        end
        println(expr)
        @lazyconstraint(cb,expr<= 0)
        println("The master problem to solve:")
        print(PM)
    end
    addlazycallback(PM, cut)
    solve(PM)
    obj=getobjectivevalue(PM)
    𝐲=getvalue(y)
end

b=[2,3]
D=[-1 3 ; 1 3]
A=[-1 2 ; 1 1]
c=[2,3]
d=[1,4]
obj,y,out,x=solvebenders(c,d,A,D,b)

```

It should finish with the solution \lambda=(0,0), x=(0,0), y=(0,1),z=4, but returns  
WARNING: Not solved to optimality, status: Unbounded  
WARNING: Unbounded ray not available  
WARNING: Variable value not defined for component of y. Check that the model was properly solved.

---

<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:** [March 27, 2018, 8:05pm UTC](https://discourse.julialang.org/t/cut-generation-not-working-lazy-constraints-lagrangian-relaxation/9948/4 "2018-03-27T20:05:02Z")

</div>

When I run your example, it enters the cut callback.

To prevent unboundedness on the initial iterations, you need to place a lower bound on the cost-to-go variable (i.e. z).

You may also want to look up `push!` and `append!` ([Expressions and Constraints — JuMP -- Julia for Mathematical Optimization 0.18 documentation](http://www.juliaopt.org/JuMP.jl/0.18/refexpr.html?highlight=append#methods)) as more efficient ways to build your cut expression than using `+=`.

Moreover, since you are using Gurobi, it might be good to have a read of [https://github.com/JuliaOpt/Gurobi.jl#reusing-the-same-gurobi-environment-for-multiple-solves](https://github.com/JuliaOpt/Gurobi.jl#reusing-the-same-gurobi-environment-for-multiple-solves)

---

<div class="post-metadata">

**Author:** ![tlagos](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tlagos/32/4490_2.png) [@tlagos](https://discourse.julialang.org/u/tlagos)\
**Post date:** [March 27, 2018, 8:10pm UTC](https://discourse.julialang.org/t/cut-generation-not-working-lazy-constraints-lagrangian-relaxation/9948/5 "2018-03-27T20:10:01Z")

</div>

Finally the problem was solver dependent.  
When I use CPLEX this works fine, but when I use Gurobi it stops before adding any new lazy constraint.

```
using JuMP
using CPLEX
slvr=CplexSolver()

# \min \{y_1 + 4y_2 + 2x_1 + 3x_2 \}
# s.t. y_1 - 3y_2 + x_1 -2x_2 \leq -2
# y_1 + 3y_2 + x_1 + x_2 \geq 3
# x_1,x_2\geq 0 , y_1,y_2 positive integers
#Initially we choose y_1^1 = y_2^1 = 0 , z^1=-\infty

function definemaster(d,D,b,λ)
    PM = Model(solver=slvr)
    @variable(PM, z)
    @variable(PM, y[i in 1:length(d)] >= 0, Int)
    λD= [sum(D[j,i]*λ[j] for j=1:length(b)) for i=1:length(d)]
    @constraint(PM,Ep[k in 1:1],z>= sum(y[i]*(d[i]-λD[i]) for i =1:length(d)) + sum(λ[j]*b[j] for j =1:length(b)) )
    # println(Ep[1])
    @objective(PM, :Min, z)
    return PM,y,z
end

function solvesubproblem(b,Dy,A,c)
    sub=Model(solver=slvr)
    @variable(sub,λ[j in 1:length(b)]>=0)
    @constraint(sub,x[i =1:length(c)], sum( A[j,i] * λ[j] for j =1:length(b) ) <= c[i])
    @objective(sub,:Max,sum( (b[j]-Dy[j])*λ[j] for j =1:length(b) ))
    status=solve(sub)
    if status==:Optimal
        return :Optimal,getvalue(λ),getobjectivevalue(sub),getdual(x)
    elseif status==:Unbounded
        return :Unbounded,getvalue(λ)
    else
        println("I'm in a really bad shape, status is ",status)
        return status
    end
end

function solvebenders(c,d,A,D,b)
    out=solvesubproblem(b,zeros(b),A,c)
    PM,y,z=definemaster(d,D,b,out[2])
    function cut(cb)
        println("----\nInside cut callback")
        masterval=getvalue(z);
        println("Current objective value of master problem is: ", masterval)
        𝐲=getvalue(y)
        println("Solution for 𝐲 is: ", 𝐲)
        Dy=[sum(D[j,i]*𝐲[i] for i=1:length(d)) for j=1:length(b)]
        out=solvesubproblem(b,Dy,A,c)
        λ=out[2]
        λD= [sum(D[j,i]*λ[j] for j=1:length(b)) for i=1:length(d)]
        expr=zero(AffExpr)
        if out[1]==:Unbounded
            expr+= sum( b[j]*λ[j] for j =1:length(b) )
            for j in 1:length(d)
                push!(expr,-λD[j],y[j])
            end
        elseif out[1]==:Optimal
            if -0.0001+ sum( b[j]*λ[j] for j =1:length(b)) + sum( (d[j]-λD[j])*𝐲[j] for j = 1:length(d) )<=masterval<=0.0001+ sum( b[j]*λ[j] for j =1:length(b)) + sum( (d[j]-λD[j])*𝐲[j] for j = 1:length(d) )
                println("masterval=",masterval)
                println( sum( b[j]*λ[j] for j =1:length(b)) + sum( (d[j]-λD[j])*𝐲[j] for j = 1:length(d)))
                return
            end
            expr+= sum( b[j]*λ[j] for j =1:length(b) )
            for j in 1:length(d)
                push!(expr,d[j]-λD[j],y[j])
            end
            push!(expr,-1.0,z)
        else
            println("please break")
        end
        println(expr)
        @lazyconstraint(cb,expr<= 0)
        println("The master problem to solve:")
        print(PM)
    end
    addlazycallback(PM, cut)
    solve(PM)
    obj=getobjectivevalue(PM)
    𝐲=getvalue(y)
    Dy=[sum(D[j,i]*𝐲[i] for i=1:length(d)) for j=1:length(b)]
    out=solvesubproblem(b,Dy,A,c)
    return obj,𝐲,out
end

b=[2,3]
D=[-1 3 ; 1 3]
A=[-1 2 ; 1 1]
c=[2,3]
d=[1,4]
obj,y,out=solvebenders(c,d,A,D,b)

```

---

<div class="post-metadata">

**Author:** ![tlagos](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tlagos/32/4490_2.png) [@tlagos](https://discourse.julialang.org/u/tlagos)\
**Post date:** [March 27, 2018, 8:41pm UTC](https://discourse.julialang.org/t/cut-generation-not-working-lazy-constraints-lagrangian-relaxation/9948/6 "2018-03-27T20:41:13Z")

</div>

With Gurobi it was fixed by adding the lower bound to the z:

```
function definemaster(d,D,b,λ)
    PM = Model(solver=slvr)
    @variable(PM, z>=-1.0)
    @variable(PM, y[i in 1:length(d)] >= 0, Int)
    λD= [sum(D[j,i]*λ[j] for j=1:length(b)) for i=1:length(d)]
    @constraint(PM,Ep[k in 1:1],z>= sum(y[i]*(d[i]-λD[i]) for i =1:length(d)) + sum(λ[j]*b[j] for j =1:length(b)) )
    # println(Ep[1])
    @objective(PM, :Min, z)
    return PM,y,z
end

```

Thank you very much Oscar!

---

<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:** [March 27, 2018, 9:54pm UTC](https://discourse.julialang.org/t/cut-generation-not-working-lazy-constraints-lagrangian-relaxation/9948/7 "2018-03-27T21:54:03Z")

</div>

This is one reason that JuMP 0.19 will no longer support solver-independent callbacks – they aren’t solver independent. See:

> [@Cutting plane algorithm performance varies with Gurobi and CPLEX](https://discourse.julialang.org/t/cutting-plane-algorithm-performance-varies-with-gurobi-and-cplex/8597/2):
>
> Callbacks are tied in very deeply to how a MIP solver works, and they’re just not the same across solvers. (This took some experience to discover.) In JuMP 0.19 we will no longer be providing solver-independent callbacks for this reason; they create a false expectation that you can abstract away the solver when you implement the callbacks. The replacement will be solver-specific callbacks that expose the same interface as the low-level API.