diff --git a/src/algorithms/colgen.jl b/src/algorithms/colgen.jl index d308057d9..8733743df 100644 --- a/src/algorithms/colgen.jl +++ b/src/algorithms/colgen.jl @@ -1,33 +1,61 @@ Base.@kwdef struct ColumnGeneration <: AbstractAlgorithm - option::Bool = false + max_nb_iterations::Int = 1000 + optimality_tol::Float64 = 1e-5 end -mutable struct ColumnGenTmpRecord +# Data stored while algorithm is running +mutable struct ColGenRuntimeData incumbents::Incumbents has_converged::Bool is_feasible::Bool + params::ColumnGeneration end -function ColumnGenTmpRecord(S::Type{<:AbstractObjSense}, node_inc::Incumbents) - i = Incumbents(S) - update_ip_primal_sol!(i, get_ip_primal_sol(node_inc)) - return ColumnGenTmpRecord(i, false, true) +function ColGenRuntimeData( + algparams::ColumnGeneration, form::Reformulation, node::Node +) + inc = Incumbents(form.master.obj_sense) + update_ip_primal_sol!(inc, get_ip_primal_sol(node.incumbents)) + return ColGenRuntimeData(inc, false, true, algparams) end # Data needed for another round of column generation -mutable struct ColumnGenerationRecord <: AbstractAlgorithmResult +mutable struct ColumnGenerationResult <: AbstractAlgorithmResult incumbents::Incumbents proven_infeasible::Bool end # Overload of the algorithm's prepare function -function prepare!(algo::ColumnGeneration, form, node) +function prepare!(alg::ColumnGeneration, form, node) @logmsg LogLevel(-1) "Prepare ColumnGeneration." return end -function should_do_ph_1(cg_rec::ColumnGenerationRecord) - primal_lp_sol = getsol(get_lp_primal_sol(cg_rec.incumbents)) +# Overload of the algorithm's run function +function run!(alg::ColumnGeneration, form::Reformulation, node::Node) + @logmsg LogLevel(-1) "Run ColumnGeneration." + algdata = ColGenRuntimeData(alg, form, node) + result = cg_main_loop(algdata, form, 2) + if should_do_ph_1(result) + record!(form, node) + set_ph_one(form.master) + result = cg_main_loop(algdata, form, 1) + end + if result.proven_infeasible + result.incumbents = Incumbents(getsense(result.incumbents)) + end + if result.proven_infeasible + @logmsg LogLevel(-1) "ColumnGeneration terminated with status INFEASIBLE." + else + @logmsg LogLevel(-1) "ColumnGeneration terminated with status FEASIBLE." + end + update!(node.incumbents, result.incumbents) + return result +end + +# Internal methods to the column generation +function should_do_ph_1(result::ColumnGenerationResult) + primal_lp_sol = getsol(get_lp_primal_sol(result.incumbents)) art_vars = filter(x->(getduty(x) isa ArtificialDuty), primal_lp_sol) if !isempty(art_vars) @logmsg LogLevel(-2) "Artificial variables in lp solution, need to do phase one" @@ -45,36 +73,11 @@ function set_ph_one(master::Formulation) return end -function run!(algo::ColumnGeneration, form, node) - @logmsg LogLevel(-1) "Run ColumnGeneration." - algdata = ColumnGenTmpRecord(form.master.obj_sense, node.incumbents) - cg_rec = cg_main_loop(algdata, form, 2) - if should_do_ph_1(cg_rec) - record!(form, node) - set_ph_one(form.master) - cg_rec = cg_main_loop(algdata, form, 1) - end - if cg_rec.proven_infeasible - cg_rec.incumbents = Incumbents(getsense(cg_rec.incumbents)) - end - if cg_rec.proven_infeasible - @logmsg LogLevel(-1) "ColumnGeneration terminated with status INFEASIBLE." - else - @logmsg LogLevel(-1) "ColumnGeneration terminated with status FEASIBLE." - end - update!(node.incumbents, cg_rec.incumbents) - return cg_rec -end - -# Internal methods to the column generation function update_pricing_problem!(spform::Formulation, dual_sol::DualSolution) - - masterform = spform.parent_formulation - + masterform = getmaster(spform) for (var_id, var) in Iterators.filter(_active_pricing_sp_var_ , getvars(spform)) setcurcost!(spform, var, computereducedcost(masterform, var_id, dual_sol)) end - return false end @@ -82,26 +85,25 @@ function update_pricing_target!(spform::Formulation) # println("pricing target will only be needed after automating convexity constraints") end -function insert_cols_in_master!(masterform::Formulation, - spform::Formulation, - sp_sols::Vector{PrimalSolution{S}}) where {S} - +function insert_cols_in_master!( + masterform::Formulation, spform::Formulation, + spsols::Vector{PrimalSolution{S}} +) where {S} sp_uid = getuid(spform) nb_of_gen_col = 0 - - for sp_sol in sp_sols - if contrib_improves_mlp(getbound(sp_sol)) + for spsol in spsols + if contrib_improves_mlp(getbound(spsol)) nb_of_gen_col += 1 ref = getvarcounter(masterform) + 1 name = string("MC", sp_uid, "_", ref) - resetsolvalue!(masterform, sp_sol) + resetsolvalue!(masterform, spsol) lb = 0.0 ub = Inf kind = Continuous duty = MasterCol sense = Positive mc = setprimaldwspsol!( - masterform, name, sp_sol, duty; lb = lb, ub = ub, + masterform, name, spsol, duty; lb = lb, ub = ub, kind = kind, sense = sense ) @logmsg LogLevel(-2) string("Generated column : ", name) @@ -122,17 +124,16 @@ function insert_cols_in_master!(masterform::Formulation, ==# end end - return nb_of_gen_col end contrib_improves_mlp(sp_primal_bound::PrimalBound{MinSense}) = (sp_primal_bound < 0.0 - 1e-8) contrib_improves_mlp(sp_primal_bound::PrimalBound{MaxSense}) = (sp_primal_bound > 0.0 + 1e-8) -function compute_pricing_db_contrib(spform::Formulation, - sp_sol_primal_bound::PrimalBound{S}, - sp_lb::Float64, - sp_ub::Float64) where {S} +function compute_pricing_db_contrib( + spform::Formulation, sp_sol_primal_bound::PrimalBound{S}, sp_lb::Float64, + sp_ub::Float64 +) where {S} # Since convexity constraints are not automated and there is no stab # the pricing_dual_bound_contrib is just the reduced cost * multiplicty if contrib_improves_mlp(sp_sol_primal_bound) @@ -143,12 +144,10 @@ function compute_pricing_db_contrib(spform::Formulation, return contrib end -function solve_sp_to_gencol!(masterform::Formulation, - spform::Formulation, - dual_sol::DualSolution, - sp_lb::Float64, - sp_ub::Float64) - +function solve_sp_to_gencol!( + masterform::Formulation, spform::Formulation, dual_sol::DualSolution, + sp_lb::Float64, sp_ub::Float64 +) #flag_need_not_generate_more_col = 0 # Not used flag_is_sp_infeasible = -1 #flag_cannot_generate_more_col = -2 # Not used @@ -198,19 +197,19 @@ function solve_sp_to_gencol!(masterform::Formulation, return insertion_status, pricing_db_contrib end -function solve_sps_to_gencols!(reformulation::Reformulation, - dual_sol::DualSolution{S}, - sp_lbs::Dict{FormId, Float64}, - sp_ubs::Dict{FormId, Float64}) where {S} - +function solve_sps_to_gencols!( + reform::Reformulation, dual_sol::DualSolution{S}, + sp_lbs::Dict{FormId, Float64}, sp_ubs::Dict{FormId, Float64} +) where {S} nb_new_cols = 0 dual_bound_contrib = DualBound{S}(0.0) - masterform = getmaster(reformulation) - sps = get_dw_pricing_sp(reformulation) + masterform = getmaster(reform) + sps = get_dw_pricing_sp(reform) for spform in sps sp_uid = getuid(spform) - gen_status, contrib = solve_sp_to_gencol!(masterform, spform, dual_sol, sp_lbs[sp_uid], sp_ubs[sp_uid]) - + gen_status, contrib = solve_sp_to_gencol!( + masterform, spform, dual_sol, sp_lbs[sp_uid], sp_ubs[sp_uid] + ) if gen_status > 0 nb_new_cols += gen_status dual_bound_contrib += float(contrib) @@ -221,15 +220,17 @@ function solve_sps_to_gencols!(reformulation::Reformulation, return (nb_new_cols, dual_bound_contrib) end -function compute_master_db_contrib(algdata::ColumnGenTmpRecord, - restricted_master_sol_value::PrimalBound{S}) where {S} +function compute_master_db_contrib( + algdata::ColGenRuntimeData, restricted_master_sol_value::PrimalBound{S} +) where {S} # TODO: will change with stabilization return DualBound{S}(restricted_master_sol_value) end -function update_lagrangian_db!(algdata::ColumnGenTmpRecord, - restricted_master_sol_value::PrimalBound{S}, - pricing_sp_dual_bound_contrib::DualBound{S}) where {S} +function update_lagrangian_db!( + algdata::ColGenRuntimeData, restricted_master_sol_value::PrimalBound{S}, + pricing_sp_dual_bound_contrib::DualBound{S} +) where {S} lagran_bnd = DualBound{S}(0.0) lagran_bnd += compute_master_db_contrib(algdata, restricted_master_sol_value) lagran_bnd += pricing_sp_dual_bound_contrib @@ -245,8 +246,10 @@ function solve_restricted_master!(master::Formulation) getprimalsols(opt_result), getdualsols(opt_result), elapsed_time) end -function generatecolumns!(algdata::ColumnGenTmpRecord, reform::Reformulation, - master_val, dual_sol, sp_lbs, sp_ubs) +function generatecolumns!( + algdata::ColGenRuntimeData, reform::Reformulation, master_val, + dual_sol, sp_lbs, sp_ubs +) nb_new_columns = 0 while true # TODO Replace this condition when starting implement stabilization nb_new_col, sp_db_contrib = solve_sps_to_gencols!(reform, dual_sol, sp_lbs, sp_ubs) @@ -264,21 +267,21 @@ end ph_one_infeasible_db(db::DualBound{MinSense}) = getvalue(db) > (0.0 + 1e-5) ph_one_infeasible_db(db::DualBound{MaxSense}) = getvalue(db) < (0.0 - 1e-5) -function cg_main_loop(algdata::ColumnGenTmpRecord, - reformulation::Reformulation, - phase::Int)::ColumnGenerationRecord +function cg_main_loop( + algdata::ColGenRuntimeData, reform::Reformulation, phase::Int +)::ColumnGenerationResult nb_cg_iterations = 0 # Phase II loop: Iterate while can generate new columns and # termination by bound does not apply - masterform = reformulation.master + masterform = reform.master sp_lbs = Dict{FormId, Float64}() sp_ubs = Dict{FormId, Float64}() # collect multiplicity current bounds for each sp - for spform in reformulation.dw_pricing_subprs + for spform in reform.dw_pricing_subprs sp_uid = getuid(spform) - lb_convexity_constr_id = reformulation.dw_pricing_sp_lb[sp_uid] - ub_convexity_constr_id = reformulation.dw_pricing_sp_ub[sp_uid] + lb_convexity_constr_id = reform.dw_pricing_sp_lb[sp_uid] + ub_convexity_constr_id = reform.dw_pricing_sp_ub[sp_uid] sp_lbs[sp_uid] = getcurrhs(getconstr(masterform, lb_convexity_constr_id)) sp_ubs[sp_uid] = getcurrhs(getconstr(masterform, ub_convexity_constr_id)) end @@ -290,7 +293,7 @@ function cg_main_loop(algdata::ColumnGenTmpRecord, if (phase != 1 && (master_status == MOI.INFEASIBLE || master_status == MOI.INFEASIBLE_OR_UNBOUNDED)) @error "Solver returned that restricted master LP is infeasible or unbounded (status = $master_status) during phase != 1." - return ColumnGenerationRecord(algdata.incumbents, true) + return ColumnGenerationResult(algdata.incumbents, true) end update_lp_primal_sol!(algdata.incumbents, primal_sols[1]) @@ -306,13 +309,13 @@ function cg_main_loop(algdata::ColumnGenTmpRecord, # generate new columns by solving the subproblems sp_time = @elapsed begin nb_new_col = generatecolumns!( - algdata, reformulation, master_val, dual_sols[1], sp_lbs, sp_ubs + algdata, reform, master_val, dual_sols[1], sp_lbs, sp_ubs ) end if nb_new_col < 0 @error "Infeasible subproblem." - return ColumnGenerationRecord(algdata.incumbents, true) + return ColumnGenerationResult(algdata.incumbents, true) end print_intermediate_statistics( @@ -323,42 +326,41 @@ function cg_main_loop(algdata::ColumnGenTmpRecord, dual_bound = get_ip_dual_bound(algdata.incumbents) primal_bound = get_lp_primal_bound(algdata.incumbents) - cur_gap = gap(primal_bound, dual_bound) - ip_primal_bound = get_ip_primal_bound(algdata.incumbents) - if diff(dual_bound, ip_primal_bound) < -1e-6 + if diff(dual_bound, ip_primal_bound) < algdata.params.optimality_tol algdata.has_converged = false # ??? @logmsg LogLevel(1) "Dual bound reached primal bound." - return ColumnGenerationRecord(algdata.incumbents, false) + return ColumnGenerationResult(algdata.incumbents, false) end if phase == 1 && ph_one_infeasible_db(dual_bound) algdata.is_feasible = false @logmsg LogLevel(1) "Phase one determines infeasibility." - return ColumnGenerationRecord(algdata.incumbents, true) + return ColumnGenerationResult(algdata.incumbents, true) end - if nb_new_col == 0 || cur_gap < 0.00001 #_params_.relative_optimality_tolerance - @logmsg LogLevel(1) "Column Generation Algorithm has converged." #nb_new_col cur_gap + if nb_new_col == 0 || gap(primal_bound, dual_bound) < algdata.params.optimality_tol + @logmsg LogLevel(1) "Column Generation Algorithm has converged." algdata.has_converged = true - return ColumnGenerationRecord(algdata.incumbents, false) + return ColumnGenerationResult(algdata.incumbents, false) end - if nb_cg_iterations > 1000 ##TDalgdata.max_nb_cg_iterations + if nb_cg_iterations > algdata.params.max_nb_iterations @warn "Maximum number of column generation iteration is reached." - return ColumnGenerationRecord(algdata.incumbents, false) + return ColumnGenerationResult(algdata.incumbents, false) end end - return ColumnGenerationRecord(algdata.incumbents, false) + return ColumnGenerationResult(algdata.incumbents, false) end -function print_intermediate_statistics(algdata::ColumnGenTmpRecord, - nb_new_col::Int, - nb_cg_iterations::Int, - mst_time::Float64, sp_time::Float64) +function print_intermediate_statistics( + algdata::ColGenRuntimeData, nb_new_col::Int, nb_cg_iterations::Int, + mst_time::Float64, sp_time::Float64 +) mlp = getvalue(get_lp_primal_bound(algdata.incumbents)) db = getvalue(get_ip_dual_bound(algdata.incumbents)) pb = getvalue(get_ip_primal_bound(algdata.incumbents)) @printf( - " \n", - nb_cg_iterations, _elapsed_solve_time(), mst_time, sp_time, nb_new_col, mlp, db, pb + " \n", + nb_cg_iterations, _elapsed_solve_time(), mst_time, sp_time, nb_new_col, mlp, db, pb ) + return end diff --git a/test/full_instances_tests.jl b/test/full_instances_tests.jl index f75d173fe..d582a6ade 100644 --- a/test/full_instances_tests.jl +++ b/test/full_instances_tests.jl @@ -43,6 +43,32 @@ function generalized_assignment_tests() @test CLD.GeneralizedAssignment.print_and_check_sol(data, problem, x) end + @testset "gap - ColGen max nb iterations" begin + data = CLD.GeneralizedAssignment.data("smallgap3.txt") + + coluna = JuMP.with_optimizer( + Coluna.Optimizer, params = CL.Params( + global_strategy = CL.GlobalStrategy( + CL.SimpleBnP( + colgen = CL.ColumnGeneration( + max_nb_iterations = 8 + ) + ), + CL.SimpleBranching(), + CL.DepthFirst() + ) + ), + default_optimizer = with_optimizer(GLPK.Optimizer) + ) + + problem, x, dec = CLD.GeneralizedAssignment.model(data, coluna) + + JuMP.optimize!(problem) + @test abs(JuMP.objective_value(problem) - 438.0) <= 0.00001 + @test MOI.get(problem.moi_backend.optimizer, MOI.TerminationStatus()) == MOI.OPTIMAL + @test CLD.GeneralizedAssignment.print_and_check_sol(data, problem, x) + end + @testset "gap with penalties - pure master variables" begin data = CLD.GeneralizedAssignment.data("smallgap3.txt")