From 94e87573432d64559ba5dde5a6f2a622d3ffb4dd Mon Sep 17 00:00:00 2001 From: Guillaume Marques Date: Thu, 4 May 2023 14:35:47 +0200 Subject: [PATCH] Classical Benders algorithm (#860) * wip * loop * Update test/unit/Benders/benders_default.jl --- src/Algorithm/Algorithm.jl | 1 + src/Algorithm/benders.jl | 3 - src/Algorithm/benders/default.jl | 179 +++++++++++++++++++++++++-- src/Algorithm/benders/printer.jl | 144 +++++++++++++++++++++ src/Algorithm/benders/utils.jl | 4 +- src/Benders/Benders.jl | 64 ++++++++-- src/MathProg/MOIinterface.jl | 28 +++++ src/MathProg/solutions.jl | 35 +++++- test/unit/Benders/benders_default.jl | 46 +++++-- 9 files changed, 463 insertions(+), 41 deletions(-) create mode 100644 src/Algorithm/benders/printer.jl diff --git a/src/Algorithm/Algorithm.jl b/src/Algorithm/Algorithm.jl index ea3377ff0..6cc2fe454 100644 --- a/src/Algorithm/Algorithm.jl +++ b/src/Algorithm/Algorithm.jl @@ -47,6 +47,7 @@ include("colgen.jl") # Benders algorithm include("benders/utils.jl") include("benders/default.jl") +include("benders/printer.jl") include("benders.jl") # Presolve diff --git a/src/Algorithm/benders.jl b/src/Algorithm/benders.jl index ff88a2f0d..095095ed3 100644 --- a/src/Algorithm/benders.jl +++ b/src/Algorithm/benders.jl @@ -27,9 +27,6 @@ function get_units_usage(algo::BendersCutGeneration, reform::Reformulation) return units_usage end - - - mutable struct BendersCutGenRuntimeData optstate::OptimizationState spform_phase::Dict{FormId, FormulationPhase} diff --git a/src/Algorithm/benders/default.jl b/src/Algorithm/benders/default.jl index ba02758b1..b4d50c9e0 100644 --- a/src/Algorithm/benders/default.jl +++ b/src/Algorithm/benders/default.jl @@ -1,11 +1,11 @@ -struct BendersContext +struct BendersContext <: Benders.AbstractBendersContext reform::Reformulation optim_sense restr_master_solve_alg restr_master_optimizer_id::Int - nb_benders_ietartion_limits::Int + nb_benders_iteration_limits::Int rhs_helper::RhsCalculationHelper - + second_stage_cost_var_ids::Vector{VarId} separation_solve_alg function BendersContext(reform, alg) @@ -16,15 +16,26 @@ struct BendersContext alg.restr_master_optimizer_id, alg.max_nb_iterations, RhsCalculationHelper(reform), + _second_stage_cost_var_ids(reform), alg.separation_solve_alg ) end end +function _second_stage_cost_var_ids(reform) + var_ids = VarId[] + for (_, sp) in get_benders_sep_sps(reform) + id = sp.duty_data.second_stage_cost_var + @assert !isnothing(id) + push!(var_ids, id) + end + return var_ids +end + +Benders.is_minimization(ctx::BendersContext) = ctx.optim_sense == MinSense 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} result::OptimizationState{F,S} end @@ -35,11 +46,65 @@ function Benders.is_unbounded(master_res::BendersMasterResult) end Benders.get_primal_sol(master_res::BendersMasterResult) = get_best_lp_primal_sol(master_res.result) +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)) + + +function _reset_second_stage_cost_var_inc_vals(ctx::BendersContext) + for var_id in ctx.second_stage_cost_var_ids + setcurincval!(Benders.get_master(ctx), var_id, 0.0) + end + return +end + +function _update_second_stage_cost_var_inc_vals(ctx::BendersContext, master_res::BendersMasterResult) + isnothing(Benders.get_primal_sol(master_res)) && return + _reset_second_stage_cost_var_inc_vals(ctx) + for (var_id, val) in Benders.get_primal_sol(master_res) + setcurincval!(Benders.get_master(ctx), var_id, val) + end + return +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) - return BendersMasterResult(opt_state) + master_res = BendersMasterResult(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 + + # We can derive a cut from the extreme ray + certificates = MathProg.get_dual_infeasibility_certificate(master, getoptimizer(master, ctx.restr_master_optimizer_id)) + + if length(certificates) > 0 + opt_state = OptimizationState( + master; + termination_status = UNBOUNDED + ) + set_ip_primal_sol!(opt_state, first(certificates)) + set_lp_primal_sol!(opt_state, first(certificates)) + mast_result = BendersMasterResult(opt_state) + _update_second_stage_cost_var_inc_vals(ctx, mast_result) + certificate = true + else + # If there is no dual infeasibility certificate, we set the cost of the second stage + # cost variable to zero and solve the master. + # TODO: This trick can backfire on us if the optimizer finds the master unbounded + # and does not return any dual infeasibility certificate for several consecutive iterations. + # It this case, we can end up with the same first level solution over and over again + # and probably be trapped in an infinite loop. + # We can escape the infinite loop by implementing a cut duplication checker but the + # algorithm won't exit gracefully. + Benders.set_second_stage_var_costs_to_zero!(ctx) + mast_result = Benders.optimize_master_problem!(master, ctx, env) + Benders.reset_second_stage_var_costs!(ctx) + end + return mast_result, certificate end function Benders.set_second_stage_var_costs_to_zero!(ctx::BendersContext) @@ -75,6 +140,20 @@ 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) + spid = getuid(sp) + T = ctx.rhs_helper.T[spid] + + new_rhs = T * mast_primal_sol + + for (constr_id, constr) in getconstrs(sp) + if getduty(constr_id) <= BendSpTechnologicalConstr + setcurrhs!(sp, constr_id, - new_rhs[constr_id]) + end + end + return +end + struct GeneratedCut lhs::Dict{VarId, Float64} rhs::Float64 @@ -90,13 +169,17 @@ end Base.iterate(set::CutsSet) = iterate(set.cuts) Base.iterate(set::CutsSet, state) = iterate(set.cuts, state) -Benders.set_of_cuts(ctx::BendersContext) = CutsSet() +Benders.set_of_cuts(::BendersContext) = CutsSet() struct BendersSeparationResult{F,S} + second_stage_estimation::Float64 + second_stage_cost::Float64 result::OptimizationState{F,S} cut::GeneratedCut end +Benders.get_primal_sol(res::BendersSeparationResult) = get_best_lp_primal_sol(res.result) + function Benders.optimize_separation_problem!(ctx::BendersContext, sp::Formulation{BendersSp}, env) spid = getuid(sp) input = OptimizationState(sp) @@ -116,22 +199,35 @@ function Benders.optimize_separation_problem!(ctx::BendersContext, sp::Formulati cut_rhs = transpose(dual_sol) * ctx.rhs_helper.rhs[spid] cut = GeneratedCut(cut_lhs, cut_rhs) - return BendersSeparationResult(opt_state, cut) + + cost = getvalue(dual_sol) + 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) + + return BendersSeparationResult(estimated_cost, cost, opt_state, cut) end Benders.get_dual_sol(res::BendersSeparationResult) = get_best_lp_dual_sol(res.result) function Benders.push_in_set!(ctx::BendersContext, set::CutsSet, sep_result::BendersSeparationResult) - @show sep_result.cut - push!(set.cuts, sep_result.cut) - return true + sc = Benders.is_minimization(ctx) ? 1.0 : -1.0 + + eq = abs(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 + + # if cost of separation result > second cost variable in master result + if !eq && gt + push!(set.cuts, sep_result.cut) + return true + end + return false end function Benders.insert_cuts!(reform, ctx::BendersContext, cuts) master = Benders.get_master(ctx) - @show cuts - cuts_to_insert = cuts.cuts #GeneratedCut[] + cuts_to_insert = cuts.cuts # We add the new cuts cut_ids = ConstrId[] @@ -143,6 +239,63 @@ function Benders.insert_cuts!(reform, ctx::BendersContext, cuts) ) push!(cut_ids, getid(constr)) end - return cut_ids +end + +struct BendersIterationOutput <: Benders.AbstractBendersIterationOutput + min_sense::Bool + nb_new_cuts::Int + infeasible_master::Bool + infeasible_subproblem::Bool + time_limit_reached::Bool + master::Union{Nothing,Float64} +end + +Benders.benders_iteration_output_type(ctx::BendersContext) = BendersIterationOutput + +function Benders.new_iteration_output( + ::Type{BendersIterationOutput}, + is_min_sense, + nb_new_cuts, + infeasible_master, + infeasible_subproblem, + time_limit_reached, + master_value +) + return BendersIterationOutput( + is_min_sense, + nb_new_cuts, + infeasible_master, + infeasible_subproblem, + time_limit_reached, + master_value + ) +end + +struct BendersOutput <: Benders.AbstractBendersOutput + infeasible_master::Bool + infeasible_subproblem::Bool + time_limit_reached::Bool +end + +Benders.benders_output_type(::BendersContext) = BendersOutput + +function Benders.new_output( + ::Type{BendersOutput}, + benders_iter_output::BendersIterationOutput +) + return BendersOutput( + benders_iter_output.infeasible_master, + benders_iter_output.infeasible_subproblem, + benders_iter_output.time_limit_reached + ) +end + +Benders.stop_benders(::BendersContext, ::Nothing, benders_iteration) = false +function Benders.stop_benders(ctx::BendersContext, benders_iteration_output::BendersIterationOutput, benders_iteration) + return benders_iteration_output.infeasible_master || + benders_iteration_output.infeasible_subproblem || + benders_iteration_output.time_limit_reached || + benders_iteration_output.nb_new_cuts <= 0 || + ctx.nb_benders_iteration_limits <= benders_iteration end \ No newline at end of file diff --git a/src/Algorithm/benders/printer.jl b/src/Algorithm/benders/printer.jl new file mode 100644 index 000000000..4a8a59b38 --- /dev/null +++ b/src/Algorithm/benders/printer.jl @@ -0,0 +1,144 @@ +mutable struct BendersPrinterContext + inner::BendersContext + sp_elapsed_time::Float64 + mst_elapsed_time::Float64 + print::Bool + debug_mode::Bool + debug_print_master::Bool + debug_print_master_primal_solution::Bool + debug_print_master_dual_solution::Bool + debug_print_subproblem::Bool + debug_print_subproblem_primal_solution::Bool + debug_print_subproblem_dual_solution::Bool + debug_print_generated_cuts::Bool +end + +function BendersPrinterContext( + reform, alg; + print = false, + debug_print_master = false, + debug_print_master_primal_solution = false, + debug_print_master_dual_solution = false, + debug_print_subproblem = false, + debug_print_subproblem_primal_solution = false, + debug_print_subproblem_dual_solution = false, + debug_print_generated_cuts = false, +) + debug_mode = debug_print_master || + debug_print_master_primal_solution || + debug_print_subproblem || + debug_print_subproblem_primal_solution || + debug_print_subproblem_dual_solution || + debug_print_generated_cuts + return BendersPrinterContext( + BendersContext(reform, alg), + 0.0, + 0.0, + print, + debug_mode, + debug_print_master, + debug_print_master_primal_solution, + debug_print_master_dual_solution, + debug_print_subproblem, + debug_print_subproblem_primal_solution, + debug_print_subproblem_dual_solution, + debug_print_generated_cuts, + ) +end + +Benders.is_minimization(ctx::BendersPrinterContext) = Benders.is_minimization(ctx.inner) +Benders.get_reform(ctx::BendersPrinterContext) = Benders.get_reform(ctx.inner) +Benders.get_master(ctx::BendersPrinterContext) = Benders.get_master(ctx.inner) +Benders.get_benders_subprobs(ctx::BendersPrinterContext) = Benders.get_benders_subprobs(ctx.inner) + +function Benders.optimize_master_problem!(master, ctx::BendersPrinterContext, 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") + end + ctx.mst_elapsed_time = @elapsed begin + result = Benders.optimize_master_problem!(master, ctx.inner, env) + end + if ctx.debug_print_master + print(crayon"bold underline blue", "Master problem:", crayon"!bold !underline") + @show master + print(crayon"reset") + end + if ctx.debug_print_master_primal_solution + print(crayon"bold underline blue", "Master primal solution:", crayon"!bold !underline") + @show Benders.get_primal_sol(result) + print(crayon"reset") + end + if ctx.debug_print_master_dual_solution + print(crayon"bold underline blue", "Master dual solution:", crayon"!bold !underline") + @show Benders.get_dual_sol(result) + print(crayon"reset") + end + 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") + end + return result +end + +Benders.treat_unbounded_master_problem!(master, ctx::BendersPrinterContext, env) = Benders.treat_unbounded_master_problem!(master, ctx.inner, env) +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.set_of_cuts(ctx::BendersPrinterContext) = Benders.set_of_cuts(ctx.inner) + +function Benders.optimize_separation_problem!(ctx::BendersPrinterContext, sp::Formulation{BendersSp}, env) + 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 + 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) + cut_ids = Benders.insert_cuts!(reform, ctx.inner, cuts) +end + +Benders.benders_iteration_output_type(ctx::BendersPrinterContext) = Benders.benders_iteration_output_type(ctx.inner) + +Benders.stop_benders(ctx::BendersPrinterContext, benders_iter_output, benders_iteration) = Benders.stop_benders(ctx.inner, benders_iter_output, benders_iteration) + +Benders.benders_output_type(ctx::BendersPrinterContext) = Benders.benders_output_type(ctx.inner) + +function _benders_iter_str(iteration, benders_iter_output, sp_time::Float64, mst_time::Float64, optim_time::Float64) + master::Float64 = benders_iter_output.master + nb_new_cuts = benders_iter_output.nb_new_cuts + return @sprintf( + " ", + iteration, optim_time, mst_time, sp_time, nb_new_cuts, master + ) +end + +function Benders.after_benders_iteration(ctx::BendersPrinterContext, phase, env, iteration, benders_iter_output) + println(_benders_iter_str(iteration, benders_iter_output, ctx.sp_elapsed_time, ctx.mst_elapsed_time, elapsed_optim_time(env))) + if ctx.debug_mode + println(crayon"bold red", repeat('-', 30), " end of iteration ", iteration, " ", repeat('-', 30), crayon"reset") + end + return +end \ No newline at end of file diff --git a/src/Algorithm/benders/utils.jl b/src/Algorithm/benders/utils.jl index 989581646..bc6afb0e1 100644 --- a/src/Algorithm/benders/utils.jl +++ b/src/Algorithm/benders/utils.jl @@ -16,14 +16,14 @@ function _add_subproblem!(rhs, T, spid, sp) for (constr_id, constr) in getconstrs(sp) if getduty(constr_id) <= BendSpTechnologicalConstr && iscuractive(sp, constr) && isexplicit(sp, constr) push!(constr_ids, constr_id) - push!(constr_rhs, getcurrhs(sp, constr_id)) + push!(constr_rhs, getperenrhs(sp, constr_id)) end end rhs[spid] = sparsevec(constr_ids, constr_rhs, Coluna.MAX_NB_ELEMS) T[spid] = _submatrix( sp, constr_id -> getduty(constr_id) <= BendSpTechnologicalConstr, - var_id -> getduty(var_id) <= BendSpPosSlackFirstStageVar + var_id -> getduty(var_id) <= BendSpFirstStageRepVar ) return end diff --git a/src/Benders/Benders.jl b/src/Benders/Benders.jl index 93505ac31..a22a50bd6 100644 --- a/src/Benders/Benders.jl +++ b/src/Benders/Benders.jl @@ -5,6 +5,8 @@ using .MustImplement abstract type AbstractBendersContext end +@mustimplement "Benders" is_minimization(context::AbstractBendersContext) = nothing + @mustimplement "Benders" get_reform(context::AbstractBendersContext) = nothing @mustimplement "Benders" get_master(context::AbstractBendersContext) = nothing @@ -18,6 +20,8 @@ abstract type AbstractBendersContext end @mustimplement "Benders" get_primal_sol(res) = nothing +# If the master is unbounded +@mustimplement "Benders" treat_unbounded_master_problem!(master, context, env) = nothing # second stage variable costs @mustimplement "Benders" set_second_stage_var_costs_to_zero!(context) = nothing @@ -27,6 +31,8 @@ abstract type AbstractBendersContext end @mustimplement "Benders" update_sp_rhs!(context, sp, mast_primal_sol) = nothing +@mustimplement "Benders" set_sp_rhs_to_zero!(context, sp, mast_primal_sol) = nothing + @mustimplement "Benders" set_of_cuts(context) = nothing @mustimplement "Benders" optimize_separation_problem!(context, sp_to_solve, env) = nothing @@ -39,25 +45,61 @@ abstract type AbstractBendersContext end @mustimplement "Benders" insert_cuts!(reform, context, generated_cuts) = nothing +@mustimplement "Benders" benders_iteration_output_type(::AbstractBendersContext) = nothing + +abstract type AbstractBendersIterationOutput end + +@mustimplement "Benders" new_iteration_output(::Type{<:AbstractBendersIterationOutput}, is_min_sense, nb_cuts_inserted, infeasible_master, infeasible_subproblem, time_limit_reached, master_obj_val) = nothing + +@mustimplement "Benders" after_benders_iteration(::AbstractBendersContext, phase, env, iteration, benders_iter_output) = nothing + +@mustimplement "Benders" stop_benders(::AbstractBendersContext, benders_iter_output, iteration) = nothing + +@mustimplement "Benders" benders_output_type(::AbstractBendersContext) = nothing + +abstract type AbstractBendersOutput end + +@mustimplement "Benders" new_output(::Type{<:AbstractBendersOutput}, benders_iter_output) = nothing + +@mustimplement "BendersMasterResult" get_obj_val(master_res) = nothing + +function run_benders_loop!(context, env; iter = 1) + iteration = iter + phase = nothing + ip_primal_sol = nothing + benders_iter_output = nothing + while !stop_benders(context, benders_iter_output, iteration) + benders_iter_output = run_benders_iteration!(context, phase, env, ip_primal_sol) + after_benders_iteration(context, phase, env, iteration, benders_iter_output) + iteration += 1 + end + O = benders_output_type(context) + return new_output(O, benders_iter_output) +end + function run_benders_iteration!(context, phase, env, ip_primal_sol) master = get_master(context) mast_result = optimize_master_problem!(master, context, env) - @show mast_result + certificate = false # At first iteration, if the master does not contain any Benders cut, the master will be - # unbounded. we therefore solve the master by setting the cost of the second stage cost - # variable to 0 so that the problem won't be unbounded anymore. + # unbounded. The implementation must provide a routine to handle this case. if is_unbounded(mast_result) - set_second_stage_var_costs_to_zero!(context) - mast_result = optimize_master_problem!(master, context, env) - reset_second_stage_var_costs!(context) + mast_result, certificate = treat_unbounded_master_problem!(master, context, env) end # Master primal solution mast_primal_sol = get_primal_sol(mast_result) for (_, sp) in get_benders_subprobs(context) - update_sp_rhs!(context, sp, mast_primal_sol) + # Right-hand-side of linking constraints is not updated in the same way whether the + # master returns a dual infeasibility certificate or a primal solution. + # See Lemma 2 of "Implementing Automatic Benders Decomposition in a Modern MIP Solver" (Bonami et al., 2020) + if certificate + set_sp_rhs_to_zero!(context, sp, mast_primal_sol) + else + update_sp_rhs!(context, sp, mast_primal_sol) + end end generated_cuts = set_of_cuts(context) @@ -76,9 +118,11 @@ function run_benders_iteration!(context, phase, env, ip_primal_sol) end cut_ids = insert_cuts!(get_reform(context), context, generated_cuts) - @show master - @show cut_ids - return + nb_cuts_inserted = length(cut_ids) + O = benders_iteration_output_type(context) + is_min_sense = is_minimization(context) + master_obj_val = get_obj_val(mast_result) + return new_iteration_output(O, is_min_sense, nb_cuts_inserted, false, false, false, master_obj_val) end end \ No newline at end of file diff --git a/src/MathProg/MOIinterface.jl b/src/MathProg/MOIinterface.jl index eaa5f8b59..1f8e7843a 100644 --- a/src/MathProg/MOIinterface.jl +++ b/src/MathProg/MOIinterface.jl @@ -370,6 +370,34 @@ function get_dual_solutions(form::F, optimizer::MoiOptimizer) where {F <: Formul return solutions end +function get_dual_infeasibility_certificate(form::F, optimizer::MoiOptimizer) where {F <: Formulation} + inner = getinner(optimizer) + nb_certificates = MOI.get(inner, MOI.ResultCount()) + certificates = PrimalSolution{F}[] + + for res_idx in 1:nb_certificates + if MOI.get(inner, MOI.PrimalStatus(res_idx)) != MOI.INFEASIBILITY_CERTIFICATE + continue + end + + # The ray is stored in the primal status. + certificate_var_ids = VarId[] + certificate_var_vals = Float64[] + for (varid, var) in getvars(form) + moi_index = getindex(getmoirecord(var)) + MOI.is_valid(inner, moi_index) || continue + val = MOI.get(inner, MOI.VariablePrimal(res_idx), moi_index) + val = round(val, digits = Coluna.TOL_DIGITS) + if abs(val) > Coluna.TOL + push!(certificate_var_ids, varid) + push!(certificate_var_vals, val) + end + end + push!(certificates, PrimalSolution(form, certificate_var_ids, certificate_var_vals, 0.0, INFEASIBLE_SOL)) + end + return certificates +end + function _show_function(io::IO, moi_model::MOI.ModelLike, func::MOI.ScalarAffineFunction) for term in func.terms diff --git a/src/MathProg/solutions.jl b/src/MathProg/solutions.jl index 322837fed..f3d17accf 100644 --- a/src/MathProg/solutions.jl +++ b/src/MathProg/solutions.jl @@ -269,4 +269,37 @@ end Base.:(*)(m::DynamicSparseMatrix, s::AbstractSolution) = m * s.solution.sol Base.:(*)(m::DynamicSparseArrays.Transposed{<:DynamicSparseMatrix}, s::AbstractSolution) = m * s.solution.sol -LinearAlgebra.norm(s::AbstractSolution) = norm(s.solution.sol) \ No newline at end of file +LinearAlgebra.norm(s::AbstractSolution) = norm(s.solution.sol) + +############################################################################################ +# Infeasibility certificate +############################################################################################ +@enum PrimalOrDualInfeasibility PRIMAL_INFEASIBILITY DUAL_INFEASIBILITY +struct InfeasibilityCertificate{M} <: AbstractSolution + infeasibility::PrimalOrDualInfeasibility + vars::Solution{M,VarId,Float64} + constrs::Solution{M,ConstrId,Float64} +end + +function PrimalInfeasibilityCertificate( + form::M, constrids, constrvals, varids, varvals, cost; +) where {M<:AbstractFormulation} + @assert length(constrids) == length(constrvals) + @assert length(varids) == length(varvals) + status = INFEASIBLE_SOL + vars = Solution{M,VarId,Float64}(form, varids, varvals, cost, status) + constrs = Solution{M,ConstrId,Float64}(form, constrids, constrvals, cost, status) + return InfeasibilityCertificate{M}(PRIMAL_INFEASIBILITY, vars, constrs) +end + +function DualInfeasibilityCertificate( + form::M, constrids, constrvals, varids, varvals, cost; +) where {M<:AbstractFormulation} + @assert length(constrids) == length(constrvals) + @assert length(varids) == length(varvals) + status = INFEASIBLE_SOL + vars = Solution{M,VarId,Float64}(form, varids, varvals, cost, status) + constrs = Solution{M,ConstrId,Float64}(form, constrids, constrvals, cost, status) + return InfeasibilityCertificate{M}(DUAL_INFEASIBILITY, vars, constrs) +end + diff --git a/test/unit/Benders/benders_default.jl b/test/unit/Benders/benders_default.jl index be592a152..5c8fa6b08 100644 --- a/test/unit/Benders/benders_default.jl +++ b/test/unit/Benders/benders_default.jl @@ -57,6 +57,15 @@ function benders_simple_example() end function benders_iteration_default() + # using JuMP, GLPK + # m = Model(GLPK.Optimizer) + # @variable(m, x[1:2]>= 0) + # @variable(m, y[1:2] >= 0) + # @constraint(m, -x[1] + 4x[2] + 2y[1] + 3y[2] >= 2) + # @constraint(m, x[1] + 3x[2] + y[1] + y[2] >= 3) + # @objective(m, Min, x[1] + 4x[2] + 2y[1] + 3y[2]) + # objectivevalue(m) + env, reform = benders_simple_example() master = Coluna.MathProg.getmaster(reform) @@ -68,18 +77,31 @@ function benders_iteration_default() ClMP.push_optimizer!(sp, () -> ClA.MoiOptimizer(GLPK.Optimizer())) end - println("\e[42m --------------- n --s------------- \e[00m") + alg = Coluna.Algorithm.BendersCutGeneration( + max_nb_iterations = 10 + ) + ctx = Coluna.Algorithm.BendersPrinterContext( + reform, alg; + 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) + + # println("\e[41m ------------1---------------- \e[00m") + # Coluna.Benders.run_benders_iteration!(ctx, nothing, env, nothing) + # println("\e[41m ------------2---------------- \e[00m") + # Coluna.Benders.run_benders_iteration!(ctx, nothing, env, nothing) + # println("\e[41m ------------3---------------- \e[00m") + # Coluna.Benders.run_benders_iteration!(ctx, nothing, env, nothing) + # println("\e[41m ------------4---------------- \e[00m") + # Coluna.Benders.run_benders_iteration!(ctx, nothing, env, nothing) - alg = Coluna.Algorithm.BendersCutGeneration() - ctx = Coluna.Algorithm.BendersContext(reform, alg) - println("\e[42m --------------- n --1------------- \e[00m") - Coluna.Benders.run_benders_iteration!(ctx, nothing, env, nothing) - println("\e[42m --------------- n --2------------- \e[00m") - Coluna.Benders.run_benders_iteration!(ctx, nothing, env, nothing) - println("\e[42m --------------- n --3------------- \e[00m") - Coluna.Benders.run_benders_iteration!(ctx, nothing, env, nothing) - println("\e[42m --------------- n --4------------- \e[00m") - Coluna.Benders.run_benders_iteration!(ctx, nothing, env, nothing) + Coluna.Benders.run_benders_loop!(ctx, env) end -register!(unit_tests, "benders_default", benders_iteration_default) \ No newline at end of file +register!(unit_tests, "benders_default", benders_iteration_default)