Skip to content

Commit

Permalink
fix Benders algorithm logic (#911)
Browse files Browse the repository at this point in the history 
* fix Benders algorithm logic

* forgot to commit a file
  • Loading branch information
@guimarqu
guimarqu authored May 27, 2023
1 parent 00f3bda commit 70aa230
Show file tree
Hide file tree
Showing 4 changed files with 194 additions and 83 deletions.
146 changes: 95 additions & 51 deletions src/Algorithm/benders/default.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -96,7 +96,7 @@ function Benders.optimize_master_problem!(master, ctx::BendersContext, env)
return master_res
end

function Benders.treat_unbounded_master_problem!(master, ctx::BendersContext, env)
function Benders.treat_unbounded_master_problem_case!(master, ctx::BendersContext, env)
mast_result = nothing
ip_solver = typeof(ctx.restr_master_solve_alg) <: SolveIpForm

Expand Down Expand Up @@ -174,7 +174,7 @@ function Benders.update_sp_rhs!(ctx::BendersContext, sp, mast_primal_sol)
return
end

function Benders.set_sp_rhs_to_zero!(ctx::BendersContext, sp, mast_primal_sol)
function Benders.setup_separation_for_unbounded_master_case!(ctx::BendersContext, sp, mast_primal_sol)
spid = getuid(sp)
T = ctx.rhs_helper.T[spid]

Expand Down Expand Up @@ -238,29 +238,40 @@ struct BendersSeparationResult{F,S}
unbounded_master::Bool
end

Benders.get_obj_val(res::BendersSeparationResult) = res.second_stage_cost
Benders.get_primal_sol(res::BendersSeparationResult) = get_best_lp_primal_sol(res.result)
Benders.is_infeasible(res::BendersSeparationResult) = res.infeasible
Benders.is_unbounded(res::BendersSeparationResult) = res.unbounded

# This is a kind of phase 1 for the separation problem when it is infeasible.
function _optimize_feasibility_separation_problem!(ctx, sp::Formulation{BendersSp}, env)
for (varid, _) in getvars(sp)
if !isanArtificialDuty(getduty(varid))
setcurcost!(sp, varid, 0.0)
end
function _compute_cut_lhs(ctx, sp, dual_sol, feasibility_cut)
cut_lhs = Dict{VarId, Float64}()

coeffs = transpose(ctx.rhs_helper.T[getuid(sp)]) * dual_sol
for (varid, coeff) in zip(findnz(coeffs)...)
cut_lhs[varid] = coeff
end
_activate_art_vars(sp)

input = OptimizationState(sp)
opt_state = run!(ctx.separation_solve_alg, env, sp, input)
if feasibility_cut
cut_lhs[sp.duty_data.second_stage_cost_var] = 0.0
else
cut_lhs[sp.duty_data.second_stage_cost_var] = 1.0
end
return cut_lhs
end

for (varid, _) in getvars(sp)
if !isanArtificialDuty(getduty(varid))
setcurcost!(sp, varid, getperencost(sp, varid))
function _compute_cut_rhs_contrib(ctx, sp, dual_sol)
spid = getuid(sp)
bounds_contrib_to_rhs = 0.0
for (varid, (val, active_bound)) in get_var_redcosts(dual_sol)
if active_bound == MathProg.LOWER || active_bound == MathProg.LOWER_AND_UPPER
bounds_contrib_to_rhs += val * getcurlb(sp, varid)
elseif active_bound == MathProg.UPPER
bounds_contrib_to_rhs += val * getcurub(sp, varid)
end
end
_deactivate_art_vars(sp)
return opt_state

cut_rhs = transpose(dual_sol) * ctx.rhs_helper.rhs[spid] + bounds_contrib_to_rhs
return cut_rhs
end

function Benders.optimize_separation_problem!(ctx::BendersContext, sp::Formulation{BendersSp}, env, unbounded_master)
Expand All @@ -274,59 +285,88 @@ function Benders.optimize_separation_problem!(ctx::BendersContext, sp::Formulati
opt_state = run!(ctx.separation_solve_alg, env, sp, input)

if getterminationstatus(opt_state) == UNBOUNDED
return BendersSeparationResult(estimated_cost, nothing, opt_state, true, false, nothing, nothing, unbounded_master)
end

feasibility_cut = false
infeasible = getterminationstatus(opt_state) == INFEASIBLE
if infeasible
feasibility_cut = true
opt_state = _optimize_feasibility_separation_problem!(ctx, sp, env)
return BendersSeparationResult(estimated_cost, nothing, opt_state, false, true, nothing, nothing, unbounded_master)
end

infeasible = getterminationstatus(opt_state) == INFEASIBLE
if infeasible
return BendersSeparationResult(estimated_cost, nothing, opt_state, false, true, nothing, nothing, unbounded_master)
if getterminationstatus(opt_state) == INFEASIBLE
return BendersSeparationResult(estimated_cost, nothing, opt_state, true, false, nothing, nothing, unbounded_master)
end

dual_sol = get_best_lp_dual_sol(opt_state)
cost = getvalue(dual_sol)
min_sense = Benders.is_minimization(ctx)
sc = min_sense ? 1.0 : - 1.0
if sc * cost < sc * estimated_cost
# Unbounded error if in Master unbounded case
if !feasibility_cut && unbounded_master
return BendersSeparationResult(estimated_cost, nothing, opt_state, false, true, nothing, nothing, true)

# println(">>>>> $cost < $estimated_cost ")
# if sc * cost < sc * estimated_cost + 1e-5
# # Unbounded error if in Master unbounded case
# if unbounded_master
# return BendersSeparationResult(estimated_cost, nothing, opt_state, false, true, nothing, nothing, true)
# else
# # Optimal solution in the othercase
# return BendersSeparationResult(estimated_cost, cost, opt_state, false, false, dual_sol, nothing, unbounded_master)
# end
# end

cut_lhs = _compute_cut_lhs(ctx, sp, dual_sol, false)
cut_rhs = _compute_cut_rhs_contrib(ctx, sp, dual_sol)

cut = GeneratedCut(min_sense, cut_lhs, cut_rhs)
return BendersSeparationResult(estimated_cost, cost, opt_state, false, false, dual_sol, cut, unbounded_master)
end

function Benders.master_is_unbounded(ctx::BendersContext, second_stage_cost, unbounded_master_case)
estimated_cost = 0
for (spid, sp) in Benders.get_benders_subprobs(ctx)
second_stage_cost_var = sp.duty_data.second_stage_cost_var
estimated_cost += getcurincval(Benders.get_master(ctx), second_stage_cost_var)
end

min_sense = Benders.is_minimization(ctx)
sc = min_sense ? 1.0 : - 1.0
return sc * second_stage_cost < sc * estimated_cost + 1e-5 && unbounded_master_case
end

function Benders.treat_infeasible_separation_problem_case!(ctx::BendersContext, sp::Formulation{BendersSp}, env, unbounded_master_case)
second_stage_cost_var = sp.duty_data.second_stage_cost_var
@assert !isnothing(second_stage_cost_var)
estimated_cost = getcurincval(Benders.get_master(ctx), second_stage_cost_var)

for (varid, _) in getvars(sp)
if !isanArtificialDuty(getduty(varid))
setcurcost!(sp, varid, 0.0)
end
# Optimal solution in the othercase
end
_activate_art_vars(sp)

cut_lhs = Dict{VarId, Float64}()
input = OptimizationState(sp)
opt_state = run!(ctx.separation_solve_alg, env, sp, input)

coeffs = transpose(ctx.rhs_helper.T[getuid(sp)]) * dual_sol
for (varid, coeff) in zip(findnz(coeffs)...)
cut_lhs[varid] = coeff
for (varid, _) in getvars(sp)
if !isanArtificialDuty(getduty(varid))
setcurcost!(sp, varid, getperencost(sp, varid))
end
end
_deactivate_art_vars(sp)

if feasibility_cut
cut_lhs[sp.duty_data.second_stage_cost_var] = 0.0
else
cut_lhs[sp.duty_data.second_stage_cost_var] = 1.0
if getterminationstatus(opt_state) == INFEASIBLE # should not happen
error("A")
return BendersSeparationResult(estimated_cost, nothing, opt_state, false, true, nothing, nothing, unbounded_master_case)
end

bounds_contrib_to_rhs = 0.0
for (varid, (val, active_bound)) in get_var_redcosts(dual_sol)
if active_bound == MathProg.LOWER || active_bound == MathProg.LOWER_AND_UPPER
bounds_contrib_to_rhs += val * getcurlb(sp, varid)
elseif active_bound == MathProg.UPPER
bounds_contrib_to_rhs += val * getcurub(sp, varid)
end
dual_sol = get_best_lp_dual_sol(opt_state)
cost = getvalue(dual_sol)
min_sense = Benders.is_minimization(ctx)
sc = min_sense ? 1.0 : - 1.0

if sc * cost <= 0
error("B") # should not happen
end

cut_rhs = transpose(dual_sol) * ctx.rhs_helper.rhs[spid] + bounds_contrib_to_rhs
cut_lhs = _compute_cut_lhs(ctx, sp, dual_sol, true)
cut_rhs = _compute_cut_rhs_contrib(ctx, sp, dual_sol)
cut = GeneratedCut(min_sense, cut_lhs, cut_rhs)

return BendersSeparationResult(estimated_cost, cost, opt_state, false, false, dual_sol, cut, unbounded_master)
return BendersSeparationResult(estimated_cost, cost, opt_state, false, false, dual_sol, cut, unbounded_master_case)
end

function Benders.get_dual_sol(res::BendersSeparationResult)
Expand All @@ -337,9 +377,13 @@ function Benders.get_dual_sol(res::BendersSeparationResult)
end

function Benders.push_in_set!(ctx::BendersContext, set::CutsSet, sep_result::BendersSeparationResult)
if isnothing(sep_result.cut)
return false
end

sc = Benders.is_minimization(ctx) ? 1.0 : -1.0
eq = false #bs(sep_result.second_stage_cost - sep_result.second_stage_estimation) < 1e-5
gt = sc * sep_result.second_stage_cost > sc * sep_result.second_stage_estimation
gt = sc * sep_result.second_stage_cost + 1e-5 > sc * sep_result.second_stage_estimation

# if cost of separation result > second cost variable in master result
if !eq && gt
Expand Down
39 changes: 35 additions & 4 deletions src/Algorithm/benders/printer.jl
View file
Original file line number Diff line number Diff line change
Expand Up @@ -79,8 +79,8 @@ function Benders.optimize_master_problem!(master, ctx::BendersPrinterContext, en
return result
end

function Benders.treat_unbounded_master_problem!(master, ctx::BendersPrinterContext, env)
result = Benders.treat_unbounded_master_problem!(master, ctx.inner, env)
function Benders.treat_unbounded_master_problem_case!(master, ctx::BendersPrinterContext, env)
result = Benders.treat_unbounded_master_problem_case!(master, ctx.inner, env)
if ctx.debug_print_master || ctx.debug_print_master_primal_solution || ctx.debug_print_master_dual_solution
println(crayon"bold blue", repeat('-', 80), crayon"reset")
println(crayon"bold underline blue", "Treat unbounded master", crayon"reset")
Expand All @@ -96,7 +96,7 @@ end
Benders.set_second_stage_var_costs_to_zero!(ctx::BendersPrinterContext) = Benders.set_second_stage_var_costs_to_zero!(ctx.inner)
Benders.reset_second_stage_var_costs!(ctx::BendersPrinterContext) = Benders.reset_second_stage_var_costs!(ctx.inner)
Benders.update_sp_rhs!(ctx::BendersPrinterContext, sp, primal_sol) = Benders.update_sp_rhs!(ctx.inner, sp, primal_sol)
Benders.set_sp_rhs_to_zero!(ctx::BendersPrinterContext, sp, primal_sol) = Benders.set_sp_rhs_to_zero!(ctx.inner, sp, primal_sol)
Benders.setup_separation_for_unbounded_master_case!(ctx::BendersPrinterContext, sp, primal_sol) = Benders.setup_separation_for_unbounded_master_case!(ctx.inner, sp, primal_sol)
Benders.set_of_cuts(ctx::BendersPrinterContext) = Benders.set_of_cuts(ctx.inner)
Benders.set_of_sep_sols(ctx::BendersPrinterContext) = Benders.set_of_sep_sols(ctx.inner)

Expand All @@ -105,7 +105,7 @@ function Benders.optimize_separation_problem!(ctx::BendersPrinterContext, sp::Fo
println(crayon"bold green", repeat('-', 80), crayon"reset")
end
if ctx.debug_print_subproblem
print(crayon"bold underline green", "Separation problem:", crayon"!bold !underline")
print(crayon"bold underline green", "Separation problem (unbounded master = $unbounded_master):", crayon"!bold !underline")
@show sp
print(crayon"reset")
end
Expand All @@ -128,6 +128,35 @@ function Benders.optimize_separation_problem!(ctx::BendersPrinterContext, sp::Fo
return result
end

function Benders.treat_infeasible_separation_problem_case!(ctx::BendersPrinterContext, sp, env, unbounded_master_case)
result = Benders.treat_infeasible_separation_problem_case!(ctx.inner, sp, env, unbounded_master_case)
if ctx.debug_print_subproblem || ctx.debug_print_subproblem_primal_solution || ctx.debug_print_subproblem_dual_solution
println(crayon"bold green", repeat('-', 80), crayon"reset")
end
if ctx.debug_print_subproblem
print(crayon"bold underline green", "Phase 1 Separation problem (unbounded_master = $unbounded_master_case):", crayon"!bold !underline")
@show sp
print(crayon"reset")
end
ctx.sp_elapsed_time = @elapsed begin
result = Benders.treat_infeasible_separation_problem_case!(ctx.inner, sp, env, unbounded_master_case)
end
if ctx.debug_print_subproblem_primal_solution
print(crayon"bold underline green", "Separation problem primal solution:", crayon"!bold !underline")
@show Benders.get_primal_sol(result)
print(crayon"reset")
end
if ctx.debug_print_subproblem_dual_solution
print(crayon"bold underline green", "Separation problem dual solution:", crayon"!bold !underline")
@show Benders.get_dual_sol(result)
print(crayon"reset")
end
if ctx.debug_print_subproblem || ctx.debug_print_subproblem_primal_solution || ctx.debug_print_subproblem_dual_solution
println(crayon"bold green", repeat('-', 80), crayon"reset")
end
return result
end

Benders.push_in_set!(ctx::BendersPrinterContext, set, sep_result) = Benders.push_in_set!(ctx.inner, set, sep_result)

function Benders.insert_cuts!(reform, ctx::BendersPrinterContext, cuts)
Expand Down Expand Up @@ -160,3 +189,5 @@ end
function Benders.build_primal_solution(context::BendersPrinterContext, mast_primal_sol, sep_sp_sols)
return Benders.build_primal_solution(context.inner, mast_primal_sol, sep_sp_sols)
end

Benders.master_is_unbounded(ctx::BendersPrinterContext, second_stage_cost, unbounded_master_case) = Benders.master_is_unbounded(ctx.inner, second_stage_cost, unbounded_master_case)
Loading

0 comments on commit 70aa230

Please sign in to comment.