From ae28e27a5902f3cee03d75d2bce68888542e0c01 Mon Sep 17 00:00:00 2001 From: Guillaume Marques Date: Thu, 4 May 2023 22:42:55 +0200 Subject: [PATCH] fix benders when the master has integral variables (#866) --- src/Algorithm/basic/solveipform.jl | 2 +- src/Algorithm/benders/default.jl | 26 ++++++++++++++++++++++---- src/Algorithm/benders/printer.jl | 21 +++++++++++++++++---- test/unit/Benders/benders_default.jl | 24 +++++++++++++++--------- 4 files changed, 55 insertions(+), 18 deletions(-) diff --git a/src/Algorithm/basic/solveipform.jl b/src/Algorithm/basic/solveipform.jl index 9bb651934..15f781e3e 100644 --- a/src/Algorithm/basic/solveipform.jl +++ b/src/Algorithm/basic/solveipform.jl @@ -49,7 +49,7 @@ _optimizer_params(::Formulation, algo::SolveIpForm, ::UserOptimizer) = algo.user _optimizer_params(form::Formulation, algo::SolveIpForm, ::CustomOptimizer) = getinner(getoptimizer(form, algo.optimizer_id)) _optimizer_params(::Formulation, ::SolveIpForm, ::NoOptimizer) = nothing -function run!(algo::SolveIpForm, env::Env, form::Formulation, input::OptimizationState) +function run!(algo::SolveIpForm, env::Env, form::Formulation, input::OptimizationState, optimizer_id = 1) opt = getoptimizer(form, algo.optimizer_id) params = _optimizer_params(form, algo, opt) if params !== nothing diff --git a/src/Algorithm/benders/default.jl b/src/Algorithm/benders/default.jl index 2f200c2e6..fae251076 100644 --- a/src/Algorithm/benders/default.jl +++ b/src/Algorithm/benders/default.jl @@ -37,6 +37,7 @@ Benders.get_reform(ctx::BendersContext) = ctx.reform Benders.get_master(ctx::BendersContext) = getmaster(ctx.reform) Benders.get_benders_subprobs(ctx::BendersContext) = get_benders_sep_sps(ctx.reform) struct BendersMasterResult{F,S} + ip_solver::Bool result::OptimizationState{F,S} end @@ -45,9 +46,15 @@ function Benders.is_unbounded(master_res::BendersMasterResult) return status == ClB.UNBOUNDED end -Benders.get_primal_sol(master_res::BendersMasterResult) = get_best_lp_primal_sol(master_res.result) +function Benders.get_primal_sol(master_res::BendersMasterResult) + if master_res.ip_solver + return get_best_ip_primal_sol(master_res.result) + end + return get_best_lp_primal_sol(master_res.result) +end + Benders.get_dual_sol(master_res::BendersMasterResult) = get_best_lp_dual_sol(master_res.result) -Benders.get_obj_val(master_res::BendersMasterResult) = getvalue(get_best_lp_primal_sol(master_res.result)) +Benders.get_obj_val(master_res::BendersMasterResult) = getvalue(Benders.get_primal_sol(master_res)) function _reset_second_stage_cost_var_inc_vals(ctx::BendersContext) @@ -69,7 +76,8 @@ end function Benders.optimize_master_problem!(master, ctx::BendersContext, env) rm_input = OptimizationState(master) opt_state = run!(ctx.restr_master_solve_alg, env, master, rm_input, ctx.restr_master_optimizer_id) - master_res = BendersMasterResult(opt_state) + ip_solver = typeof(ctx.restr_master_solve_alg) <: SolveIpForm + master_res = BendersMasterResult(ip_solver, opt_state) _update_second_stage_cost_var_inc_vals(ctx, master_res) return master_res end @@ -77,6 +85,16 @@ end function Benders.treat_unbounded_master_problem!(master, ctx::BendersContext, env) mast_result = nothing certificate = false + ip_solver = typeof(ctx.restr_master_solve_alg) <: SolveIpForm + + # In the unbounded case, to get a dual infeasibility certificate, we need to relax the + # integrality and solve the master again. (at least with GLPK & the current implementation of SolveIpForm) + if ip_solver + relax_integrality!(master) + rm_input = OptimizationState(master) + opt_state = run!(SolveLpForm(get_dual_solution = true), env, master, rm_input, ctx.restr_master_optimizer_id) + enforce_integrality!(master) + end # We can derive a cut from the extreme ray certificates = MathProg.get_dual_infeasibility_certificate(master, getoptimizer(master, ctx.restr_master_optimizer_id)) @@ -88,7 +106,7 @@ function Benders.treat_unbounded_master_problem!(master, ctx::BendersContext, en ) set_ip_primal_sol!(opt_state, first(certificates)) set_lp_primal_sol!(opt_state, first(certificates)) - mast_result = BendersMasterResult(opt_state) + mast_result = BendersMasterResult(ip_solver, opt_state) _update_second_stage_cost_var_inc_vals(ctx, mast_result) certificate = true else diff --git a/src/Algorithm/benders/printer.jl b/src/Algorithm/benders/printer.jl index 4a8a59b38..e4020eaf3 100644 --- a/src/Algorithm/benders/printer.jl +++ b/src/Algorithm/benders/printer.jl @@ -79,7 +79,20 @@ function Benders.optimize_master_problem!(master, ctx::BendersPrinterContext, en return result end -Benders.treat_unbounded_master_problem!(master, ctx::BendersPrinterContext, env) = Benders.treat_unbounded_master_problem!(master, ctx.inner, env) +function Benders.treat_unbounded_master_problem!(master, ctx::BendersPrinterContext, env) + result, c = Benders.treat_unbounded_master_problem!(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") + @show master + print(crayon"bold underline blue", "Master primal solution:", crayon"!bold !underline") + @show Benders.get_primal_sol(result) + print(crayon"reset") + println(crayon"bold blue", repeat('-', 80), crayon"reset") + end + return result, c +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) @@ -90,14 +103,14 @@ function Benders.optimize_separation_problem!(ctx::BendersPrinterContext, sp::Fo 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 - ctx.sp_elapsed_time = @elapsed begin - result = Benders.optimize_separation_problem!(ctx.inner, sp, env) - end if ctx.debug_print_subproblem print(crayon"bold underline green", "Separation problem:", crayon"!bold !underline") @show sp print(crayon"reset") end + ctx.sp_elapsed_time = @elapsed begin + result = Benders.optimize_separation_problem!(ctx.inner, sp, env) + 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) diff --git a/test/unit/Benders/benders_default.jl b/test/unit/Benders/benders_default.jl index 1c766c704..eee8d976d 100644 --- a/test/unit/Benders/benders_default.jl +++ b/test/unit/Benders/benders_default.jl @@ -206,18 +206,24 @@ function benders_iter_default_A_integer() alg = Coluna.Algorithm.BendersCutGeneration( max_nb_iterations = 10, - separation_solve_alg = Coluna.Algorithm.SolveIpForm() + restr_master_solve_alg = Coluna.Algorithm.SolveIpForm() ) ctx = Coluna.Algorithm.BendersPrinterContext( reform, alg; - print = true + print = true, + # debug_print_master = true, + # debug_print_master_primal_solution = true, + # debug_print_master_dual_solution = true, + # debug_print_subproblem = true, + # debug_print_subproblem_primal_solution = true, + # debug_print_subproblem_dual_solution = true ) Coluna.set_optim_start_time!(env) result = Coluna.Benders.run_benders_loop!(ctx, env) @test result.mlp ≈ 4.0 end -register!(unit_tests, "benders_default", benders_iter_default_A_integer; x = true) +register!(unit_tests, "benders_default", benders_iter_default_A_integer) # B with continuous first stage finds optimal solution @@ -254,9 +260,9 @@ register!(unit_tests, "benders_default", benders_iter_default_B_continuous) # B with integer first stage finds optimal solution # Error occurs during test TODO fix # expected output: -# mlp = 7.083333333333333 -# x1 = 0.0, x2 = 1.0 -# y1 = 3.1666666666666665, y2 = 0.3333333333333333 +# mlp = 7 +# x1 = 0.0, x2 = 2.0 +# y1 = 2, y2 = 0 function benders_iter_default_B_integer() #env, reform = benders_simple_example() env, reform = benders_form_B() @@ -271,7 +277,7 @@ function benders_iter_default_B_integer() alg = Coluna.Algorithm.BendersCutGeneration( max_nb_iterations = 10, - separation_solve_alg = Coluna.Algorithm.SolveIpForm() + restr_master_solve_alg = Coluna.Algorithm.SolveIpForm() ) ctx = Coluna.Algorithm.BendersPrinterContext( reform, alg; @@ -280,6 +286,6 @@ function benders_iter_default_B_integer() Coluna.set_optim_start_time!(env) result = Coluna.Benders.run_benders_loop!(ctx, env) - @test result.mlp ≈ 7.083333333333333 + @test result.mlp ≈ 7 end -register!(unit_tests, "benders_default", benders_iter_default_B_integer; x = true) \ No newline at end of file +register!(unit_tests, "benders_default", benders_iter_default_B_integer) \ No newline at end of file