# Cut generation not working lazy constraints - lagrangian relaxation

Hello,

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

# 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

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.

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

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
end
println(expr)
@lazyconstraint(cb,expr<= 0)
println("The master problem to solve:")
print(PM)
end
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.

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

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
end
println(expr)
@lazyconstraint(cb,expr<= 0)
println("The master problem to solve:")
print(PM)
end
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)


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!

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

1 Like