fix benders when the master has integral variables (#866)

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

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,14 +76,25 @@ 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
 
 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

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)

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)
+register!(unit_tests, "benders_default", benders_iter_default_B_integer)