I am trying to implement Benders Decomposition for a graph problem, I have a master problem and subproblem and in each iteration, I have to solve the master problem change the RHS of the subproblem constraint and solve the subproblem and add lazy constraint until both problems converge, my problem that I can’t access the objective of the master problem with the following error
LoadError: MathOptInterface.OptimizeInProgress{MathOptInterface.TerminationStatus}: Cannot get the result as the the MOI.optimize! has not finished.
Moreover, to change the constraint I am deleting the constraint and creating new ones with the new values, there is any way to do the following in a more efficient way ??
using JuMP
using Random
import Gurobi
######################################
initial = 0
n = 32
s = 33
E = [(0, 1), (7, 9), (15, 27), (11, 1), (14, 9), (17, 16), (2, 33), (33, 17),
(1, 10), (7, 27), (15, 28), (19, 1), (15, 9), (18, 16), (33, 2), (18, 33),
(1, 11), (9, 14), (16, 17), (20, 1), (16, 9), (19, 16), (3, 33), (33, 18),
(1, 19), (9, 15), (16, 18), (3, 2), (17, 9), (20, 16), (33, 3), (19, 33),
(1, 20), (9, 16), (16, 19), (4, 2), (20, 9), (22, 16), (4, 33), (33, 19),
(2, 3), (9, 17), (16, 20), (5, 2), (21, 9), (25, 16), (33, 4), (20, 33),
(2, 4), (9, 20), (16, 22), (6, 2), (22, 9), (26, 16), (5, 33), (33, 20),
(2, 5), (9, 21), (16, 25), (7, 2), (23, 9), (27, 16), (33, 5), (21, 33),
(2, 6), (9, 22), (16, 26), (8, 2), (26, 9), (32, 16), (6, 33), (33, 21),
(2, 7), (9, 23), (16, 27), (9, 2), (27, 9), (23, 17), (33, 6), (22, 33),
(2, 8), (9, 26), (16, 32), (5, 3), (32, 9), (19, 18), (7, 33), (33, 22),
(2, 9), (9, 27), (17, 23), (9, 3), (11, 10), (20, 19), (33, 7), (23, 33),
(3, 5), (9, 32), (18, 19), (13, 3), (12, 10), (25, 19), (8, 33), (33, 23),
(3, 9), (10, 11), (19, 20), (14, 3), (13, 10), (26, 19), (33, 8), (24, 33),
(3, 13), (10, 12), (19, 25), (15, 3), (13, 12), (29, 19), (9, 33), (33, 24),
(3, 14), (10, 13), (19, 26), (16, 3), (26, 12), (30, 19), (33, 9), (25, 33),
(3, 15), (12, 13), (19, 29), (17, 3), (28, 12), (21, 20), (10, 33), (33, 25),
(3, 16), (12, 26), (19, 30), (9, 5), (19, 13), (28, 20), (33, 10), (26, 33),
(3, 17), (12, 28), (20, 21), (15, 5), (25, 13), (23, 22), (11, 33), (33, 26),
(5, 9), (13, 19), (20, 28), (23, 5), (26, 13), (24, 22), (33, 11), (27, 33),
(5, 15), (13, 25), (22, 23), (25, 5), (28, 13), (25, 22), (12, 33), (33, 27),
(5, 23), (13, 26), (22, 24), (29, 5), (29, 13), (26, 22), (33, 12), (28, 33),
(5, 25), (13, 28), (22, 25), (7, 6), (16, 14), (27, 22), (13, 33), (33, 28),
(5, 29), (13, 29), (22, 26), (9, 6), (16, 15), (27, 24), (33, 13), (29, 33),
(6, 7), (14, 16), (22, 27), (14, 6), (20, 15), (27, 26), (14, 33), (33, 29),
(6, 9), (15, 16), (24, 27), (21, 6), (21, 15), (30, 29), (33, 14), (30, 33),
(6, 14), (15, 20), (26, 27), (22, 6), (22, 15), (31, 29), (15, 33), (33, 30),
(6, 21), (15, 21), (29, 30), (26, 6), (24, 15), (0, 33), (33, 15), (31, 33),
(6, 22), (15, 22), (29, 31), (27, 6), (26, 15), (33, 0), (16, 33), (33, 31),
(6, 26), (15, 24), (1, 0), (9, 7), (27, 15), (1, 33), (33, 16), (32, 33),
(6, 27), (15, 26), (10, 1), (27, 7), (28, 15), (33, 1), (17, 33), (33, 32)]
E = sort(E)
###################################################
sub_model = direct_model(Gurobi.Optimizer())
set_optimizer_attribute(sub_model, "Threads", 1)
set_optimizer_attribute(sub_model, "presolve", 0)
@variable(sub_model, y[v=initial:n] >= 0)
@variable(sub_model,u[v=initial:s] >= 0)
x_H = Dict((i, j) => 0 for i in initial:s for j in initial:s if (i,j) in E)
x_H[(28, 15)] = 1
x_H[(15, s)] = 1
x_H[(s, 28)] = 1
for i = initial:n
sub_c1=@constraint(sub_model, y[i] == sum(x_H[(i,j)] for j in initial:s if (i,j) in E))
set_name(sub_c1, "sub_cons1_($(i))")
end
for (i, j) in E
if j != s
sub_c2=@constraint(sub_model,u[i] - u[j] <= n * (1 - x_H[(i, j)]) - 1)
set_name(sub_c2, "sub_cons2_($(i),$(j))")
end
end
sub_c3=@constraint(sub_model,u[s] == 1)
set_name(sub_c3, "sub_cons3")
@objective(sub_model, Min, -1*(sum(y[w] for w in initial:n)))
############################################
master_model = direct_model(Gurobi.Optimizer())
set_optimizer_attribute(master_model, "Threads", 1)
set_optimizer_attribute(master_model, "presolve", 0)
@variable(master_model, Z <= 0)
@variable(master_model,x[E],binary = true)
for (i, j) in E
m_c1=@constraint(master_model,x[(i, j)] + x[(j, i)] <= 1)
set_name(m_c1, "master_cons1_($(i),$(j))")
end
m_c2=@constraint(master_model,sum(x[(s, v)] for v in initial:n) == 1)
set_name(m_c2, "master_cons2")
m_c3=@constraint(master_model,sum(x[(v,s)] for v in initial:n) == 1)
set_name(m_c3, "master_cons3")
for (i, j) in E
if i != s && j != s
m_c4=@constraint(master_model,(sum(x[(i, g)] for g in initial:s
if (g, i) in E && g != i && g != j)+sum(x[(j, g)] for g in initial:s
if (g, j) in E && g != i && g != j)) <= 1)
set_name(m_c4, "master_cons4_($(i),$(j))")
end
end
for i in initial:n
m_c5=@constraint(master_model,sum(x[(i, j)] for j in initial:s if (i, j) in E) <= 1)
set_name(m_c5, "master_cons5_($(i))")
m_c6=@constraint(master_model,sum(x[(j,i)] for j in initial:s if (j,i) in E) <= 1)
set_name(m_c6, "master_cons6_($(i))")
m_c7=@constraint(master_model,sum(x[(j, i)] for j in initial:s if (j, i) in E) ==
sum(x[(i, j)] for j in initial:s if (i, j) in E))
set_name(m_c7, "master_con7_($(i))")
end
@objective(master_model, Min, Z)
################################################
gamma = Dict((1, v) => 0 for v in initial:n)
zeta = Dict((1,i, j) => 0 for i in initial:s for j in initial:n if (i,j) in E)
mu = Dict(1 => 0)
#################################################
optimize!(sub_model)
if termination_status(sub_model) == MOI.OPTIMAL
k = 1
for v in initial:n
gamma[k, v] = JuMP.dual(constraint_by_name(sub_model,"sub_cons1_($(v))"))
end
for (i, j) in E
if j != s
zeta[k, i, j] = JuMP.dual(constraint_by_name(sub_model,"sub_cons2_($(i),$(j))"))
end
end
mu[k] = JuMP.dual(constraint_by_name(sub_model,"sub_cons3"))
cut=@constraint(master_model, (sum(x[(i, j)] * gamma[k, i] for (i, j) in E if i != s)
+sum((n * (1 - x[(i, j)]) - 1) * zeta[k, i, j] for (i, j) in E if j != s)
+ mu[k]) <= Z)
set_name(cut, "cut_($(k))")
else
k = 1
for v in initial:n
gamma[k, v] = JuMP.dual(constraint_by_name(sub_model,"sub_cons1_($(v))"))
end
for (i, j) in E
if j != s
zeta[k, i, j] = JuMP.dual(constraint_by_name(sub_model,"sub_cons2_($(i),$(j))"))
end
end
mu[k] = JuMP.dual(constraint_by_name(sub_model,"sub_cons3"))
cut=@constraint(master_model, (sum(x[(i, j)] * gamma[k, i] for (i, j) in E if i != s)
+sum((n * (1 - x[(i, j)]) - 1) * zeta[k, i, j] for (i, j) in E if j != s)
+ mu[k]) <= Z)
set_name(cut, "cut_($(k))")
end
###################################################
function mycallback2(cb_data)
model_status = callback_node_status(cb_data, master_model)
if model_status == MOI.CALLBACK_NODE_STATUS_INTEGER
Z_val = callback_value.(Ref(cb_data), Z)
x_val = callback_value.(Ref(cb_data), x)
end
mp_solution =JuMP.objective_bound(master_model)
status = getSubSolution(x_val,mp_solution)
if status == 3
add_feasibility_cut()
elseif status == 2
add_optimality_cut()
else
println("terminate")
return
end
end
function getSubSolution(x,mas_solution)
for (a,b) in E
if b != s
delete(sub_model,constraint_by_name(sub_model,"sub_cons2_($(a),$(b))"))
sub_c2=@constraint(sub_model,u[a] - u[b] <= n * (1 - x[(a, b)]) - 1)
set_name(sub_c2, "sub_cons2_($(a),$(b))")
end
end
println("change constraint")
for a in initial:n
temp = 0
for b in initial:s
if (a,b) in E
temp+=x[(a,b)]
end
end
delete(sub_model,constraint_by_name(sub_model,"sub_cons1_($(a))"))
sub_c1=@constraint(sub_model, y[a] == temp)
set_name(sub_c1, "sub_cons1_($(a))")
end
optimize!(sub_model)
sub_solution = JuMP.objective_value(sub_model)
if termination_status(sub_model) == MOI.OPTIMAL
if abs(mas_solution) - abs(sub_solution) <= 0.0001
println("True")
return 1
else
return 2
end
elseif termination_status(sub_model) != MOI.OPTIMAL
return 3
else
println("non of the above")
end
end
function add_feasibility_cut()
println("feas_cut")
k = rand(Int, 1)
for v in initial:n
gamma[k, v] = JuMP.dual(constraint_by_name(sub_model,"sub_cons1_($(v))"))
end
for (i, j) in E
if j != s
zeta[k, i, j] = JuMP.dual(constraint_by_name(sub_model,"sub_cons2_($(i),$(j))"))
end
end
mu[k] = JuMP.dual(constraint_by_name(sub_model,"sub_cons3"))
con = @build_constraint((sum(x[(i, j)] * gamma[k, i] for (i, j) in E if i != s)
+ sum((n * (1 - x[i, j]) - 1) * zeta[k, i, j] for (i, j) in E if j != s)
+ mu[k]) <= 0)
MOI.submit(cec_model, MOI.LazyConstraint(cb_data), con)
end
function add_optimality_cut()
println("opt_cut")
k = rand(Int, 1)
for v in initial:n
gamma[k, v] = JuMP.dual(constraint_by_name(sub_model,"sub_cons1_($(v))"))
end
for (i, j) in E
if j != s
zeta[k, i, j] = JuMP.dual(constraint_by_name(sub_model,"sub_cons2_($(i),$(j))"))
end
end
mu[k] = JuMP.dual(constraint_by_name(sub_model,"sub_cons3"))
con = @build_constraint((sum(x[(i, j)] * gamma[k, i] for (i, j) in E if i != s)
+ sum((n * (1 - x[i, j]) - 1) * zeta[k, i, j] for (i, j) in E if j != s)
+ mu[k]) <= Z)
MOI.submit(cec_model, MOI.LazyConstraint(cb_data), con)
end
MOI.set(master_model, MOI.RawParameter("LazyConstraints"), 1)
MOI.set(master_model, MOI.LazyConstraintCallback(),mycallback2)
optimize!(master_model)