Benders Subproblem Infeasible

Specific Domains Optimization (Mathematical)
1

Hey @odow Here is the full code. The optimality cuts work well (in other instances). But the feasibility cuts seem not to cut the solution and keep appearing (the same cut) in the next iteration.

using JuMP, Gurobi, Printf

# Z* = 435
# instance with Benders subproblem infeasible 
p =[58	22	67	32;
    15	49	87	52;
    49	19	48	71;
    45	26	40	74;
    95	24	81	73] 

global (n, m) = size(p)   

# control parameters/inicialization 
global MAXIMUM_TIME = 3600 # 3600.0
global ABSOLUTE_OPTIMALITY_GAP = 1e-6
global upper_bound = 10^10
global inicio = time()

global H = Int(sum(p)) 
global Δ = minimum(p) 
global B = Int(round(H/Δ)) 
global P = Array{Int64}(undef, n, m) 
[P[j,k] = floor(p[j,k]/Δ)+1 for j=1:n for k=1:m]
global π = Array{Float64}(undef, n, m)
[π[j,k] = P[j,k] - p[j,k]/Δ for j=1:n for k=1:m] 

# MASTER PROBLEM 
modelo = Model(Gurobi.Optimizer)
set_optimizer_attribute(modelo, "TimeLimit", 3600) # Gurobi 

@variables modelo begin
    x[j=1:n,k=1:m,b=1:B,l=0:1], Bin 
    u[j=1:n,k=1:m,b=1:B,l=0:1] >= 0 
    θ >= 0
end

@objective(modelo, Min, θ)  

@constraints modelo begin
    c1[j=1:n,k=1:m], sum(sum(x[j,k,b,l] for l=0:1) for b=1:B) == 1
    c2[k=1:m,b=1:B], sum(sum(sum(x[j,k,a,l] for a=max(b-P[j,k]-l+2,1):b) for l=0:1) for j=1:n) <= 1
    c4[j=1:n,k=2:m], sum(sum((b*Δ*x[j,k-1,b,l]-Δ*u[j,k-1,b,l]+x[j,k-1,b,l]*p[j,k-1]) for b=1:B) for l=0:1) <= sum(sum(b*Δ*x[j,k,b,l]-Δ*u[j,k,b,l] for b=1:B) for l=0:1)
    c6[j=1:n,k=1:m,b=1:B,l=0:1], ((1-l)*(1-π[j,k])+(1/Δ))*x[j,k,b,l] <= u[j,k,b,l]
    c7[j=1:n,k=1:m,b=1:B,l=0:1], u[j,k,b,l] <= (1-l*π[j,k])*x[j,k,b,l]
end


# SUBPROBLEM 
function subproblem(x, u)
    modelosub = Model(Gurobi.Optimizer)
    set_optimizer_attribute(modelosub, "TimeLimit", 3600) # Gurobi
    set_optimizer_attribute(modelosub, "InfUnbdInfo", 1) # to get the UnbdRay attribute (senão dá erro)

    @variables modelosub begin
        v[j=1:n,k=1:m,b=1:B,l=0:1] >= 0 
        Cmax >= 0
    end

    @objective(modelosub, Min, Cmax)

    con3 = @constraint(modelosub, [k=1:m,b=1:B], sum(sum(v[j,k,b,l] for j=1:n) for l=0:1) <= 1 - sum(sum(sum(x[j,k,a,l] for a=max(b-P[j,k]-l+2,1):b-1) - u[j,k,b,l] for j=1:n) for l=0:1) ) 
    con5 = @constraint(modelosub, [j=1:n,k=1:m,b=1:B,l=0:1; b+P[j,k]+l-1<=B], v[j,k,b+P[j,k]+l-1,l] == (2-l-π[j,k])*x[j,k,b,l] - u[j,k,b,l])
    con8 = @constraint(modelosub, [j=1:n,k=1:m,b=1:B,l=0:1; b+P[j,k]+l-1<=B], - v[j,k,b+P[j,k]+l-1,l] <= (l-1)*(x[j,k,b,l] - π[j,k]*x[j,k,b,l]))
    con9 = @constraint(modelosub, [j=1:n,k=1:m,b=1:B,l=0:1; b+P[j,k]+l-1<=B], v[j,k,b+P[j,k]+l-1,l] <= ((1/Δ)- l*π[j,k])*x[j,k,b,l] ) 
    con10 = @constraint(modelosub, [j=1:n], Cmax >= sum(sum(b*Δ*x[j,m,b,l]+x[j,m,b,l]*p[j,m] - Δ*u[j,m,b,l] for b=1:B) for l=0:1))

    optimize!(modelosub)

    if termination_status(modelosub) == OPTIMAL || termination_status(modelosub) == INFEASIBLE
        return (obj = dual_objective_value(modelosub), status = termination_status(modelosub), λ1 = dual.(con3), λ3 = dual.(con5), λ6 = dual.(con8), λ7 = dual.(con9), λ8 = dual.(con10))
    end 
end

# BENDERS DECOMPOSITION METHOD 
iter = 0
cut_opt = 0
cut_fact = 0 
tempo = time() - inicio
MAXIMUM_ITER = 10
while iter < MAXIMUM_ITER 
    global tempo = time() - inicio
    global iter += 1

    optimize!(modelo)

    lower_bound = objective_value(modelo)
    x_k = value.(x)
    u_k = value.(u)
    ret = subproblem(x_k, u_k)
    upper_bound_iter = ret.obj
    global upper_bound = min(upper_bound_iter, upper_bound)
    gap = (upper_bound - lower_bound) / upper_bound

    println()
    println("===========")
    println("Iteration  Lower Bound  Upper Bound          Gap")
    println(string(iter), "         ", string(lower_bound), "\t", string(upper_bound), "                ", "$(@sprintf("%.6f",gap))")
    println("===========")

    if gap < ABSOLUTE_OPTIMALITY_GAP
        println()
        println("Terminating with the optimal solution")
        break
    end

    if ret.status == OPTIMAL 
        cut = @constraint(modelo, θ >= 
            sum(sum((1 - sum(sum(sum(x[j,k,a,l] for a=max(b-P[j,k]-l+2,1):b-1) - u[j,k,b,l] for j=1:n) for l=0:1))*ret.λ1[k,b] for k=1:m) for b=1:B) 
            + sum(sum(sum(sum(((2-l-π[j,k])*x[j,k,b,l] - u[j,k,b,l])*ret.λ3[j,k,b,l] for j=1:n if b+P[j,k]+l-1<=B) for k=1:m) for b=1:B) for l=0:1)
            + sum(sum(sum(sum(((l-1)*(x[j,k,b,l] - π[j,k]*x[j,k,b,l]))*ret.λ6[j,k,b,l] for j=1:n if b+P[j,k]+l-1<=B) for k=1:m) for b=1:B) for l=0:1)
            + sum(sum(sum(sum((((1/Δ)- l*π[j,k])*x[j,k,b,l])*ret.λ7[j,k,b,l] for j=1:n if b+P[j,k]+l-1<=B) for k=1:m) for b=1:B) for l=0:1)
            + sum(sum(sum(b*Δ*x[j,m,b,l]+x[j,m,b,l]*p[j,m] - Δ*u[j,m,b,l] for b=1:B) for l=0:1)*ret.λ8[j] for j=1:n) )
        
        global cut_opt += 1
        
        #@info "Adding the optimality cut $(cut)"

    elseif ret.status == INFEASIBLE  
        cut = @constraint(modelo, 0 <= 
            sum(sum((1 - sum(sum(sum(x[j,k,a,l] for a=max(b-P[j,k]-l+2,1):b-1) - u[j,k,b,l] for j=1:n) for l=0:1))*ret.λ1[k,b] for k=1:m) for b=1:B) 
            + sum(sum(sum(sum(((2-l-π[j,k])*x[j,k,b,l] - u[j,k,b,l])*ret.λ3[j,k,b,l] for j=1:n if b+P[j,k]+l-1<=B) for k=1:m) for b=1:B) for l=0:1)
            + sum(sum(sum(sum(((l-1)*(x[j,k,b,l] - π[j,k]*x[j,k,b,l]))*ret.λ6[j,k,b,l] for j=1:n if b+P[j,k]+l-1<=B) for k=1:m) for b=1:B) for l=0:1)
            + sum(sum(sum(sum((((1/Δ)- l*π[j,k])*x[j,k,b,l])*ret.λ7[j,k,b,l] for j=1:n if b+P[j,k]+l-1<=B) for k=1:m) for b=1:B) for l=0:1)
            + sum(sum(sum(b*Δ*x[j,m,b,l]+x[j,m,b,l]*p[j,m] - Δ*u[j,m,b,l] for b=1:B) for l=0:1)*ret.λ8[j] for j=1:n) )

        global cut_fact += 1

        @info "Adding the feasibility cut $(cut)"
    end 
    
end

#println("Z = ", objective_value(modelo))
tempo = time() - inicio
println()
println("TEMPO TOTAL = ", tempo)
println("Núm iterações = ", iter)
println("Núm cuts opt = ", cut_opt)
println("Núm cuts fact = ", cut_fact)
2

The optimality cuts work well (in other instances)

Do they? I don’t think your cut formulae are correct. For example, you don’t appear to be using the subproblem objective in your cuts.

I’ve just made a PR to update the JuMP docs: [docs] simplify Benders tutorial with copy of state variables by odow · Pull Request #3832 · jump-dev/JuMP.jl · GitHub

You can simplify the cut computation considerably if you introduce a copy of the state variable into the subproblem, and add a x_copy == x_in constraint.

Also note that is far more efficient if you do in-place re-solves of the subproblem, instead of re-building it from scratch each time. See this part of the tutorial: Benders decomposition · JuMP

1 Like
3

Thanks, @odow I will check and test in detail.

4

I’ve opened a PR to add feasibility cuts to the documentation: [docs] add feasibility cuts to Benders decomposition by odow · Pull Request #3833 · jump-dev/JuMP.jl · GitHub

5

Just closing the loop on this:

I have significantly revised the Benders decomposition tutorial: Benders decomposition · JuMP

It uses a more realistic example that doesn’t use the Matrix based input, and it now features a section on feasibility cuts: Benders decomposition · JuMP

Let me know if you have questions after reading :smile:

1 Like