Cut generation not working lazy constraints - lagrangian relaxation

[tlagos]

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.

[odow]

It’s a bit hard to figure out what is going on without a MWE.

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

[tlagos]

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.

[odow]

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

[tlagos]

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)

[tlagos]

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!

[odow]

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

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.