From d5fd81d9a76a09ec3add960ac252edce3845587f Mon Sep 17 00:00:00 2001 From: Guillaume Marques Date: Wed, 21 Jul 2021 08:17:34 +0200 Subject: [PATCH 1/6] rm support of anonymous variables & constraints (#565) --- src/decomposition.jl | 45 +++++++++++++++++++++++++++++++++++++++++- test/issues_tests.jl | 47 ++++++++++++++++++++++++++++++++++++++++++++ 2 files changed, 91 insertions(+), 1 deletion(-) diff --git a/src/decomposition.jl b/src/decomposition.jl index 76e1c4700..ca2b6191e 100644 --- a/src/decomposition.jl +++ b/src/decomposition.jl @@ -492,18 +492,61 @@ function buildformulations!( return end +# Error messages for `check_annotations`. +_err_check_annotations(id::VarId) = error(""" +A variable (id = $id) is not annotated. +Make sure you do not use anonymous variables (variable with no name declared in JuMP macro variable). +Otherwise, open an issue at https://github.com/atoptima/Coluna.jl/issues +""") + +_err_check_annotations(id::ConstrId) = error(""" +A constraint (id = $id) is not annotated. +Make sure you do not use anonymous constraints (constraint with no name declared in JuMP macro variable). +Otherwise, open an issue at https://github.com/atoptima/Coluna.jl/issues +""") + +""" +Make sure that all variables and constraints of the original formulation are +annotated. Otherwise, it returns an error. +""" +function check_annotations(prob::Problem, annotations::Annotations) + origform = get_original_formulation(prob) + + for (varid, _) in getvars(origform) + if !haskey(annotations.ann_per_var, varid) + return _err_check_annotations(varid) + end + end + + for (constrid, _) in getconstrs(origform) + if !haskey(annotations.ann_per_constr, constrid) + return _err_check_annotations(constrid) + end + end + return true +end + +""" +Reformulate the original formulation of prob according to the annotations. +The environment maintains formulation ids. +""" function reformulate!(prob::Problem, annotations::Annotations, env::Env) + # Once the original formulation built, we close the "fill mode" of the + # coefficient matrix which is a super fast writing mode compared to the default + # writing mode of the dynamic sparse matrix. if getcoefmatrix(prob.original_formulation).fillmode closefillmode!(getcoefmatrix(prob.original_formulation)) end + decomposition_tree = annotations.tree if decomposition_tree !== nothing + check_annotations(prob, annotations) root = BD.getroot(decomposition_tree) reform = Reformulation() set_reformulation!(prob, reform) buildformulations!(prob, reform, env, annotations, reform, root) relax_integrality!(getmaster(reform)) - else + else # No decomposition provided by BlockDecomposition push_optimizer!( prob.original_formulation, prob.default_optimizer_builder diff --git a/test/issues_tests.jl b/test/issues_tests.jl index da3a5a7ea..a8fae4b0d 100644 --- a/test/issues_tests.jl +++ b/test/issues_tests.jl @@ -306,6 +306,49 @@ function continuous_vars_in_sp() return end +# issue https://github.com/atoptima/Coluna.jl/issues/553 +function unsupported_anonym_constrs_vars() + coluna = JuMP.optimizer_with_attributes( + Coluna.Optimizer, + "params" => CL.Params(solver=ClA.TreeSearchAlgorithm()), + "default_optimizer" => GLPK.Optimizer + ) + + function anonymous_var_model!(m) + y = @variable(m, binary = true) + @variable(m, 0 <= x[D] <= 1) + @constraint(m, sp[d in D], x[d] <= 0.85) + @objective(m, Min, sum(x) + y) + @dantzig_wolfe_decomposition(m, dec, D) + end + + function anonymous_constr_model!(m) + @variable(m, 0 <= x[D] <= 1) + sp = @constraint(m, [d in D], x[d] <= 0.85) + @objective(m, Min, sum(x)) + @dantzig_wolfe_decomposition(m, dec, D) + end + + @axis(D, 1:5) + m = BlockModel(coluna, direct_model=true) + anonymous_var_model!(m) + @test_throws ErrorException optimize!(m) + + m = BlockModel(coluna) + anonymous_var_model!(m) + # The variable is annotated in the master. + # @test_throws ErrorException optimize!(m) + + m = BlockModel(coluna, direct_model=true) + anonymous_constr_model!(m) + @test_throws ErrorException optimize!(m) + + m = BlockModel(coluna) + anonymous_constr_model!(m) + @test_throws ErrorException optimize!(m) + return +end + function test_issues_fixed() @testset "no_decomposition" begin solve_with_no_decomposition() @@ -338,6 +381,10 @@ function test_issues_fixed() @testset "continuous_vars_in_sp" begin continuous_vars_in_sp() end + + @testset "unsupported_anonym_constrs_vars" begin + unsupported_anonym_constrs_vars() + end end test_issues_fixed() \ No newline at end of file From fd5bead996fa440ee3ac2480ac6d8c790da73274 Mon Sep 17 00:00:00 2001 From: Lara di Cavalcanti Pontes <34425678+laradicp@users.noreply.github.com> Date: Thu, 22 Jul 2021 06:42:53 -0300 Subject: [PATCH 2/6] Add support for `AffExpr, MathOptInterface.Interval{Float64}` (#566) * add MOI.set for SplitIntervalBridge and add test * more general MOI.set method for SplitIntervalBridge and simpler test Co-authored-by: Guillaume Marques --- src/MOIwrapper.jl | 9 +++++++++ test/MathOptInterface/MOI_wrapper.jl | 23 +++++++++++++++++++++++ 2 files changed, 32 insertions(+) diff --git a/src/MOIwrapper.jl b/src/MOIwrapper.jl index fe6ce7854..dd6677099 100644 --- a/src/MOIwrapper.jl +++ b/src/MOIwrapper.jl @@ -425,6 +425,15 @@ end ############################################################################################ # Attributes of constraints ############################################################################################ +function MOI.set( + model::MOI.Bridges.LazyBridgeOptimizer{Coluna.Optimizer}, attr::MOI.AbstractConstraintAttribute, + bridge::MOI.Bridges.Constraint.SplitIntervalBridge, value +) + MOI.set(model.model, attr, bridge.lower, value) + MOI.set(model.model, attr, bridge.upper, value) + return +end + function MOI.set( model::Coluna.Optimizer, ::BD.ConstraintDecomposition, constrid::MOI.ConstraintIndex, annotation::BD.Annotation diff --git a/test/MathOptInterface/MOI_wrapper.jl b/test/MathOptInterface/MOI_wrapper.jl index c64b7c177..9534b8766 100644 --- a/test/MathOptInterface/MOI_wrapper.jl +++ b/test/MathOptInterface/MOI_wrapper.jl @@ -88,6 +88,29 @@ end @test JuMP.objective_value(model) == -JuMP.objective_value(model2) end +@testset "SplitIntervalBridge" begin + coluna = optimizer_with_attributes( + Coluna.Optimizer, + "params" => Coluna.Params( + solver=Coluna.Algorithm.TreeSearchAlgorithm() + ), + "default_optimizer" => GLPK.Optimizer + ) + + @axis(M, 1:1) + J = 1:1 + + model = BlockModel(coluna) + @variable(model, x[m in M, j in J]) + @constraint(model, mult[m in M], 1 <= sum(x[m,j] for j in J) <= 2) + @objective(model, Max, sum(x[m,j] for m in M, j in J)) + + @dantzig_wolfe_decomposition(model, decomposition, M) + + optimize!(model) + @test JuMP.objective_value(model) == 2.0 +end + const UNSUPPORTED_TESTS = [ "solve_qcp_edge_cases", # Quadratic constraints not supported "delete_nonnegative_variables", # variable deletion not supported From fdd948f2a953fbbc8fa0f2a1c9de3ff3291729c4 Mon Sep 17 00:00:00 2001 From: Guillaume Marques Date: Tue, 27 Jul 2021 17:23:51 +0200 Subject: [PATCH 3/6] Fix Benders reformulation & algorithm (#567) * rm support of anonymous variables & constraints * fix bounds of BendSpSlackSecondStageCostVar * wip * wip; bug in benders algo * WIP, bug fixed in toy instance, but test of old instance does not pass * wip * fixed * ok --- src/Algorithm/benders.jl | 46 +++++++++--------- src/MathProg/MOIinterface.jl | 1 - src/MathProg/clone.jl | 1 + src/MathProg/duties.jl | 8 +++- src/MathProg/formulation.jl | 44 ++++++++++++----- src/decomposition.jl | 93 +++++++++++++++++++++++------------- test/full_instances_tests.jl | 2 +- test/interfaces/model.jl | 1 - test/issues_tests.jl | 63 ++++++++++++++++++++++-- test/runtests.jl | 1 + 10 files changed, 185 insertions(+), 75 deletions(-) diff --git a/src/Algorithm/benders.jl b/src/Algorithm/benders.jl index 70e992cf7..416824506 100644 --- a/src/Algorithm/benders.jl +++ b/src/Algorithm/benders.jl @@ -24,7 +24,6 @@ function get_units_usage(algo::BendersCutGeneration, reform::Reformulation) end return units_usage end - mutable struct BendersCutGenRuntimeData optstate::OptimizationState spform_phase::Dict{FormId, FormulationPhase} @@ -33,13 +32,6 @@ mutable struct BendersCutGenRuntimeData #slack_cost_increase_applied::Bool end -function all_sp_in_phase2(algdata::BendersCutGenRuntimeData) - for (key, phase) in algdata.spform_phase - phase != PurePhase2 && return false - end - return true -end - function BendersCutGenRuntimeData(form::Reformulation, init_optstate::OptimizationState) optstate = OptimizationState(getmaster(form)) best_ip_primal_sol = get_best_ip_primal_sol(init_optstate) @@ -61,10 +53,12 @@ function run!( end function update_benders_sp_slackvar_cost_for_ph1!(spform::Formulation) + slack_cost = getobjsense(spform) === MinSense ? 1.0 : -1.0 for (varid, var) in getvars(spform) iscuractive(spform, varid) || continue - if getduty(varid) == BendSpSlackFirstStageVar - setcurcost!(spform, var, 1.0) + if getduty(varid) <= BendSpSlackFirstStageVar + setcurcost!(spform, var, slack_cost) + setcurub!(spform, var, getperenub(spform, var)) else setcurcost!(spform, var, 0.0) end @@ -76,7 +70,7 @@ end function update_benders_sp_slackvar_cost_for_ph2!(spform::Formulation) for (varid, var) in getvars(spform) iscuractive(spform, varid) || continue - if getduty(varid) == BendSpSlackFirstStageVar + if getduty(varid) <= BendSpSlackFirstStageVar setcurcost!(spform, var, 0.0) setcurub!(spform, var, 0.0) else @@ -169,9 +163,10 @@ function record_solutions!( )::Vector{ConstrId} recorded_dual_solution_ids = Vector{ConstrId}() + sense = getobjsense(spform) === MinSense ? 1.0 : -1.0 for dual_sol in get_lp_dual_sols(spresult) - if getvalue(dual_sol) > algo.feasibility_tol + if sense * getvalue(dual_sol) > algo.feasibility_tol (insertion_status, dual_sol_id) = setdualsol!(spform, dual_sol) if insertion_status push!(recorded_dual_solution_ids, dual_sol_id) @@ -187,9 +182,8 @@ end function insert_cuts_in_master!( masterform::Formulation, spform::Formulation, sp_dualsol_ids::Vector{ConstrId}, ) - sp_uid = getuid(spform) nb_of_gen_cuts = 0 - sense = (getobjsense(masterform) == MinSense ? Greater : Less) + sense = getobjsense(masterform) == MinSense ? Greater : Less for dual_sol_id in sp_dualsol_ids nb_of_gen_cuts += 1 @@ -265,10 +259,14 @@ function solve_sp_to_gencut!( # Solve sub-problem and insert generated cuts in master # @logmsg LogLevel(-3) "optimizing benders_sp prob" TO.@timeit Coluna._to "Bender Sep SubProblem" begin - optstate = run!(SolveLpForm(get_dual_solution = true), env, spform, OptimizationInput(OptimizationState(spform))) + optstate = run!( + SolveLpForm(get_dual_solution = true, relax_integrality = true), + env, spform, OptimizationInput(OptimizationState(spform)) + ) end optresult = getoptstate(optstate) + if getterminationstatus(optresult) == INFEASIBLE # if status != MOI.OPTIMAL sp_is_feasible = false # @logmsg LogLevel(-3) "benders_sp prob is infeasible" @@ -279,7 +277,7 @@ function solve_sp_to_gencut!( benders_sp_lagrangian_bound_contrib = compute_benders_sp_lagrangian_bound_contrib(algdata, spform, optresult) primalsol = get_best_lp_primal_sol(optresult) - spsol_relaxed = contains(primalsol, varid -> getduty(varid) == BendSpSlackFirstStageVar) + spsol_relaxed = contains(primalsol, varid -> getduty(varid) <= BendSpSlackFirstStageVar) benders_sp_primal_bound_contrib = 0.0 # compute benders_sp_primal_bound_contrib which stands for the sum of nu var, @@ -294,12 +292,12 @@ function solve_sp_to_gencut!( end end end - + if - algo.feasibility_tol <= get_lp_primal_bound(optresult) <= algo.feasibility_tol # no cuts are generated since there is no violation if spsol_relaxed if algdata.spform_phase[spform_uid] == PurePhase2 - error("In PurePhase2, art var were not supposed to be in sp forlumation ") + error("In PurePhase2, art var were not supposed to be in sp formulation ") end if algdata.spform_phase[spform_uid] == PurePhase1 error("In PurePhase1, if art var were in sol, the objective should be strictly positive.") @@ -326,7 +324,6 @@ function solve_sp_to_gencut!( continue end end - else # a cut can be generated since there is a violation recorded_dual_solution_ids = record_solutions!(algo, algdata, spform, optresult) if spsol_relaxed && algo.option_increase_cost_in_hybrid_phase @@ -418,9 +415,14 @@ end function solve_relaxed_master!(master::Formulation, env::Env) elapsed_time = @elapsed begin - optresult = TO.@timeit Coluna._to "relaxed master" run!(SolveLpForm(get_dual_solution = true), env, master, OptimizationInput(OptimizationState(master))) + optstate = TO.@timeit Coluna._to "relaxed master" begin + run!( + SolveLpForm(get_dual_solution = true, relax_integrality = true), + env, master, OptimizationInput(OptimizationState(master)) + ) + end end - return optresult, elapsed_time + return optstate, elapsed_time end function generatecuts!( @@ -470,8 +472,8 @@ function bend_cutting_plane_main_loop!( cur_gap = 0.0 optoutput, master_time = solve_relaxed_master!(masterform, env) - optresult = getoptstate(optoutput) + optresult = getoptstate(optoutput) if getterminationstatus(optresult) == INFEASIBLE db = - getvalue(DualBound(masterform)) pb = - getvalue(PrimalBound(masterform)) diff --git a/src/MathProg/MOIinterface.jl b/src/MathProg/MOIinterface.jl index 817cf2c24..d7313ce46 100644 --- a/src/MathProg/MOIinterface.jl +++ b/src/MathProg/MOIinterface.jl @@ -275,7 +275,6 @@ function get_dual_solutions(form::F, optimizer::MoiOptimizer) where {F <: Formul solcost += val * getcurrhs(form, id) val = round(val, digits = Coluna.TOL_DIGITS) if abs(val) > Coluna.TOL - @logmsg LogLevel(-4) string("Constr ", constr.name, " = ", val) push!(solconstrs, id) push!(solvals, val) end diff --git a/src/MathProg/clone.jl b/src/MathProg/clone.jl index 75941ec3f..ef3b8030e 100644 --- a/src/MathProg/clone.jl +++ b/src/MathProg/clone.jl @@ -1,3 +1,4 @@ +# TODO : these methods should not be part of MathProg. function clonevar!( originform::Formulation, destform::Formulation, diff --git a/src/MathProg/duties.jl b/src/MathProg/duties.jl index 185eef956..2b87c361c 100644 --- a/src/MathProg/duties.jl +++ b/src/MathProg/duties.jl @@ -23,7 +23,11 @@ mutable struct DwSp <: AbstractSpDuty end "A Benders subproblem of a formulation decomposed using Benders." -struct BendersSp <: AbstractSpDuty end +struct BendersSp <: AbstractSpDuty + slack_to_first_stage::Dict{VarId, VarId} +end + +BendersSp() = BendersSp(Dict{VarId, VarId}()) # # Duties tree for a Variable @@ -53,6 +57,8 @@ struct BendersSp <: AbstractSpDuty end AbstractBendSpVar <= Duty{Variable} AbstractBendSpSlackMastVar <= AbstractBendSpVar BendSpSlackFirstStageVar <= AbstractBendSpSlackMastVar + BendSpPosSlackFirstStageVar <= BendSpSlackFirstStageVar + BendSpNegSlackFirstStageVar <= BendSpSlackFirstStageVar BendSpSlackSecondStageCostVar <= AbstractBendSpSlackMastVar BendSpSepVar <= AbstractBendSpVar #BendSpPureVar <= AbstractBendSpVar diff --git a/src/MathProg/formulation.jl b/src/MathProg/formulation.jl index dd002f682..dc5b5d500 100644 --- a/src/MathProg/formulation.jl +++ b/src/MathProg/formulation.jl @@ -339,8 +339,8 @@ function setcol_from_sp_primalsol!( end function setcut_from_sp_dualsol!( - masterform::Formulation, - spform::Formulation, + masterform::Formulation{BendersMaster}, + spform::Formulation{BendersSp}, dual_sol_id::ConstrId, name::String, duty::Duty{Constraint}; @@ -351,32 +351,50 @@ function setcut_from_sp_dualsol!( is_explicit::Bool = true, moi_index::MoiConstrIndex = MoiConstrIndex() ) - rhs = getdualsolrhss(spform)[dual_sol_id] + objc = getobjsense(masterform) === MinSense ? 1.0 : -1.0 + + rhs = objc * getdualsolrhss(spform)[dual_sol_id] benders_cut_id = Id{Constraint}(duty, dual_sol_id) benders_cut_data = ConstrData( rhs, Essential, sense, inc_val, is_active, is_explicit ) + benders_cut = Constraint( benders_cut_id, name; constr_data = benders_cut_data, moi_index = moi_index ) + master_coef_matrix = getcoefmatrix(masterform) sp_coef_matrix = getcoefmatrix(spform) sp_dual_sol = getdualsolmatrix(spform)[:,dual_sol_id] - for (ds_constrid, ds_constr_val) in sp_dual_sol - ds_constr = getconstr(spform, ds_constrid) - if getduty(ds_constrid) <= AbstractBendSpMasterConstr - for (master_var_id, sp_constr_coef) in @view sp_coef_matrix[ds_constrid,:] - var = getvar(spform, master_var_id) - if getduty(master_var_id) <= AbstractBendSpSlackMastVar - master_coef_matrix[benders_cut_id, master_var_id] += ds_constr_val * sp_constr_coef + for (constr_id, constr_val) in sp_dual_sol + if getduty(constr_id) <= AbstractBendSpMasterConstr + for (var_id, coef) in @view sp_coef_matrix[constr_id,:] + master_var_id = nothing + if getduty(var_id) <= BendSpPosSlackFirstStageVar + orig_var_id = get(spform.duty_data.slack_to_first_stage, var_id, nothing) + if orig_var_id === nothing + error(""" + A subproblem first level slack variable is not mapped to a first level variable. + Please open an issue at https://github.com/atoptima/Coluna.jl/issues. + """) + end + master_var_id = getid(getvar(masterform, orig_var_id)) + elseif getduty(var_id) <= BendSpSlackSecondStageCostVar + master_var_id = var_id # identity + end + + if master_var_id !== nothing + master_coef_matrix[benders_cut_id, master_var_id] += objc * constr_val * coef end end end end + _addconstr!(masterform.manager, benders_cut) + if isexplicit(masterform, benders_cut) add!(masterform.buffer, getid(benders_cut)) end @@ -464,7 +482,7 @@ function _addlocalartvar!(form::Formulation, constr::Constraint, abs_cost::Float form, name1, MasterArtVar; cost = cost, lb = 0.0, ub = Inf, kind = Continuous ) var2 = setvar!( - form, name1, MasterArtVar; cost = cost, lb = 0.0, ub = Inf, kind = Continuous + form, name2, MasterArtVar; cost = cost, lb = 0.0, ub = Inf, kind = Continuous ) push!(constr.art_var_ids, getid(var1)) push!(constr.art_var_ids, getid(var2)) @@ -607,6 +625,7 @@ function computesolvalue(form::Formulation, sol_vec::AbstractDict{Id{Variable}, return val end +# TODO : remove (unefficient & specific to an algorithm) function computereducedcost(form::Formulation, varid::Id{Variable}, dualsol::DualSolution) redcost = getperencost(form, varid) coefficient_matrix = getcoefmatrix(form) @@ -621,7 +640,8 @@ function computereducedcost(form::Formulation, varid::Id{Variable}, dualsol::Dua return redcost end -function computereducedrhs(form::Formulation, constrid::Id{Constraint}, primalsol::PrimalSolution) +# TODO : remove (unefficient & specific to Benders) +function computereducedrhs(form::Formulation{BendersSp}, constrid::Id{Constraint}, primalsol::PrimalSolution) constrrhs = getperenrhs(form,constrid) coefficient_matrix = getcoefmatrix(form) for (varid, primal_val) in primalsol diff --git a/src/decomposition.jl b/src/decomposition.jl index ca2b6191e..a91245235 100644 --- a/src/decomposition.jl +++ b/src/decomposition.jl @@ -314,32 +314,12 @@ function instantiate_orig_vars!( annotations::Annotations, sp_ann ) - masterform = getmaster(spform) if haskey(annotations.vars_per_ann, sp_ann) vars = annotations.vars_per_ann[sp_ann] for (_, var) in vars clonevar!(origform, spform, spform, var, BendSpSepVar, cost = 0.0) end end - mast_ann = get(annotations, masterform) - if haskey(annotations.vars_per_ann, mast_ann) - vars = annotations.vars_per_ann[mast_ann] - for (id, var) in vars - duty, _ = _dutyexpofbendmastvar(var, annotations, origform) - if duty == MasterBendFirstStageVar - name = "μ[$(split(getname(origform, var), "[")[end])" - setvar!( - spform, name, BendSpSlackFirstStageVar; - cost = getcurcost(origform, var), - lb = getcurlb(origform, var), - ub = getcurub(origform, var), - kind = Continuous, - is_explicit = true, - id = Id{Variable}(BendSpSlackFirstStageVar, id, getuid(masterform)) - ) - end - end - end return end @@ -374,8 +354,54 @@ function create_side_vars_constrs!( spform::Formulation{BendersSp}, origform::Formulation{Original}, ::Env, - ::Annotations + annotations::Annotations ) + spcoef = getcoefmatrix(spform) + origcoef = getcoefmatrix(origform) + + # 1st level slack variables + masterform = getmaster(spform) + mast_ann = get(annotations, masterform) + if haskey(annotations.vars_per_ann, mast_ann) + vars = annotations.vars_per_ann[mast_ann] + for (varid, var) in vars + duty, _ = _dutyexpofbendmastvar(var, annotations, origform) + if duty == MasterBendFirstStageVar + name = string("μ⁺[", split(getname(origform, var), "[")[end], "]") + slack_pos = setvar!( + spform, name, BendSpPosSlackFirstStageVar; + cost = getcurcost(origform, var), + lb = getcurlb(origform, var), + ub = getcurub(origform, var), + kind = Continuous, + is_explicit = true, + id = Id{Variable}(BendSpPosSlackFirstStageVar, varid, getuid(masterform)) + ) + + name = string("μ⁻[", split(getname(origform, var), "[")[end], "]") + slack_neg = setvar!( + spform, name, BendSpNegSlackFirstStageVar; + cost = getcurcost(origform, var), + lb = getcurlb(origform, var), + ub = getcurub(origform, var), + kind = Continuous, + is_explicit = true + ) + + spform.duty_data.slack_to_first_stage[getid(slack_pos)] = varid + spform.duty_data.slack_to_first_stage[getid(slack_neg)] = varid + + for (constrid, coeff) in @view origcoef[:, varid] + spconstr = getconstr(spform, constrid) + if spconstr !== nothing + spcoef[getid(spconstr), getid(slack_pos)] = coeff + spcoef[getid(spconstr), getid(slack_neg)] = -coeff + end + end + end + end + end + sp_has_second_stage_cost = false global_costprofit_ub = 0.0 global_costprofit_lb = 0.0 @@ -383,46 +409,45 @@ function create_side_vars_constrs!( getduty(varid) == BendSpSepVar || continue orig_var = getvar(origform, varid) cost = getperencost(origform, orig_var) - if cost > 0.00001 + if cost > 0 global_costprofit_ub += cost * getcurub(origform, orig_var) global_costprofit_lb += cost * getcurlb(origform, orig_var) - elseif cost < - 0.00001 + elseif cost < 0 global_costprofit_ub += cost * getcurlb(origform, orig_var) global_costprofit_lb += cost * getcurub(origform, orig_var) end end - if global_costprofit_ub > 0.00001 || global_costprofit_lb < - 0.00001 + if global_costprofit_ub > 0 || global_costprofit_lb < 0 sp_has_second_stage_cost = true end if sp_has_second_stage_cost - sp_coef = getcoefmatrix(spform) sp_id = getuid(spform) # Cost constraint - nu = setvar!( - spform, "ν[$sp_id]", BendSpSlackSecondStageCostVar; - cost = 1.0, - lb = - global_costprofit_lb, + nu_pos = setvar!( + spform, "ν⁺[$sp_id]", BendSpSlackSecondStageCostVar; + cost = getobjsense(spform) == MinSense ? 1.0 : -1.0, + lb = global_costprofit_lb, ub = global_costprofit_ub, kind = Continuous, is_explicit = true ) - setcurlb!(spform, nu, 0.0) - setcurub!(spform, nu, Inf) + setcurlb!(spform, nu_pos, 0.0) + setcurub!(spform, nu_pos, Inf) cost = setconstr!( spform, "cost[$sp_id]", BendSpSecondStageCostConstr; rhs = 0.0, kind = Essential, - sense = Greater, + sense = getobjsense(spform) == MinSense ? Greater : Less, is_explicit = true ) - sp_coef[getid(cost), getid(nu)] = 1.0 + spcoef[getid(cost), getid(nu_pos)] = 1.0 for (varid, _) in getvars(spform) getduty(varid) == BendSpSepVar || continue - sp_coef[getid(cost), varid] = - getperencost(origform, varid) + spcoef[getid(cost), varid] = - getperencost(origform, varid) end end return diff --git a/test/full_instances_tests.jl b/test/full_instances_tests.jl index 61712c72d..ac96692cf 100644 --- a/test/full_instances_tests.jl +++ b/test/full_instances_tests.jl @@ -4,7 +4,7 @@ function full_instances_tests() generalized_assignment_tests() capacitated_lot_sizing_tests() lot_sizing_tests() - #facility_location_tests() + # facility_location_tests() cutting_stock_tests() cvrp_tests() end diff --git a/test/interfaces/model.jl b/test/interfaces/model.jl index 3308979c3..5aa86672b 100644 --- a/test/interfaces/model.jl +++ b/test/interfaces/model.jl @@ -25,7 +25,6 @@ mutable struct KnapsackLibOptimizer <: BlockDecomposition.AbstractCustomOptimize end function Coluna.Algorithm.get_units_usage(opt::KnapsackLibOptimizer, form) # form is Coluna Formulation - println("\e[41m get units usage \e[00m") units_usage = Tuple{AbstractModel, Coluna.ColunaBase.UnitType, Coluna.ColunaBase.UnitPermission}[] # TODO : the abstract model is KnapsackLibModel (opt.model) return units_usage diff --git a/test/issues_tests.jl b/test/issues_tests.jl index a8fae4b0d..209c69460 100644 --- a/test/issues_tests.jl +++ b/test/issues_tests.jl @@ -202,7 +202,7 @@ function branching_file_completion() for line in eachline(file) for pieceofdata in split(line) regex_match = match(r"n\d+", pieceofdata) - if regex_match != nothing + if regex_match !== nothing regex_match = regex_match.match push!(existing_nodes, parse(Int, regex_match[2:length(regex_match)])) end @@ -349,6 +349,61 @@ function unsupported_anonym_constrs_vars() return end +# issue https://github.com/atoptima/Coluna.jl/issues/554 +function simple_benders() + # Test in Min sense + coluna = optimizer_with_attributes( + Coluna.Optimizer, + "params" => Coluna.Params( + solver = Coluna.Algorithm.BendersCutGeneration() + ), + "default_optimizer" => GLPK.Optimizer + ) + + model = BlockModel(coluna, direct_model=true) + + @axis(S, 1:2) + + @variable(model, x, Bin) + @variable(model, y[i in S], Bin) + @constraint(model, purefirststage, x <= 1) + @constraint(model, tech1[S[1]], y[S[1]] <= x) + @constraint(model, tech2[S[2]], y[S[2]] <= 1-x) + @constraint(model, puresecondstage[s in S], y[s] <= 1) + @objective(model, Min, -sum(y)) + + @benders_decomposition(model, decomposition, S) + + optimize!(model) + @test objective_value(model) == -1.0 + + # Test in Max sense + coluna = optimizer_with_attributes( + Coluna.Optimizer, + "params" => Coluna.Params( + solver = Coluna.Algorithm.BendersCutGeneration() + ), + "default_optimizer" => GLPK.Optimizer + ) + + model = BlockModel(coluna, direct_model=true) + + @axis(S, 1:2) + + @variable(model, x, Bin) + @variable(model, y[i in S], Bin) + @constraint(model, purefirststage, x <= 1) + @constraint(model, tech1[S[1]], y[S[1]] <= x) + @constraint(model, tech2[S[2]], y[S[2]] <= 1-x) + @constraint(model, puresecondstage[s in S], y[s] <= 1) + @objective(model, Max, +sum(y)) + + @benders_decomposition(model, decomposition, S) + + optimize!(model) + @test objective_value(model) == 1.0 +end + function test_issues_fixed() @testset "no_decomposition" begin solve_with_no_decomposition() @@ -385,6 +440,8 @@ function test_issues_fixed() @testset "unsupported_anonym_constrs_vars" begin unsupported_anonym_constrs_vars() end -end -test_issues_fixed() \ No newline at end of file + @testset "simple benders decomposition" begin + simple_benders() + end +end diff --git a/test/runtests.jl b/test/runtests.jl index 5a0b2cc0c..d50afbba8 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -40,6 +40,7 @@ include("sol_disaggregation_tests.jl") rng = MersenneTwister(1234123) unit_tests() +test_issues_fixed() @testset "Full instances " begin full_instances_tests() From ec1fa6b6893721d9ae26f1db8a964de2cda27fcc Mon Sep 17 00:00:00 2001 From: Lara di Cavalcanti Pontes <34425678+laradicp@users.noreply.github.com> Date: Mon, 9 Aug 2021 17:33:35 -0300 Subject: [PATCH 4/6] MOI tests: LP unit tests and contlinear tests (#571) * first draft of MOI.ConstraintDual for Coluna.Optimizer * fix MOI.ConstraintDual for ScalarAffine and add draft for SingleVariable * fix lists of var and constr indices and fix var primal * add MOI.DualObjectiveValue and organize unsupported/failing tests * close fill mode before viewing matrix in get MOI.ConstraintFunction * update comments Co-authored-by: Guillaume Marques --- src/MOIwrapper.jl | 49 ++++++++++++++------- test/MathOptInterface/MOI_wrapper.jl | 66 +++++++++++++++++++++++----- 2 files changed, 88 insertions(+), 27 deletions(-) diff --git a/src/MOIwrapper.jl b/src/MOIwrapper.jl index dd6677099..2a0474bd7 100644 --- a/src/MOIwrapper.jl +++ b/src/MOIwrapper.jl @@ -216,10 +216,10 @@ end function MOI.get(model::Coluna.Optimizer, ::MOI.ListOfVariableIndices) indices = Vector{MathOptInterface.VariableIndex}() - for (key,value) in model.moi_varids + for (_, value) in model.moi_varids push!(indices, value) end - return sort!(indices) + return sort!(indices, by = x -> x.value) end ############################################################################################ @@ -286,11 +286,10 @@ end function MOI.get( model::Coluna.Optimizer, ::MOI.ListOfConstraintIndices{F, S} ) where {F<:MOI.SingleVariable, S} - orig_form = get_original_formulation(model.inner) indices = MOI.ConstraintIndex{F,S}[] - for (id, var) in model.vars - if S == MathProg.convert_coluna_kind_to_moi(getperenkind(orig_form, var)) - push!(indices, MOI.ConstraintIndex{F,S}(id.value)) + for (id, _) in model.constrs_on_single_var_to_vars + if S == typeof(MOI.get(model, MOI.ConstraintSet(), id)) + push!(indices, id) end end return sort!(indices, by = x -> x.value) @@ -302,8 +301,10 @@ function MOI.get( orig_form = get_original_formulation(model.inner) constrid = getid(model.constrs[index]) terms = MOI.ScalarAffineTerm{Float64}[] - for (varid, coeff) in @view getcoefmatrix(orig_form)[constrid, :] - push!(terms, MOI.ScalarAffineTerm(coeff, model.moi_varids[varid])) + coefmatrix = getcoefmatrix(orig_form) + coefmatrix.fillmode && closefillmode!(coefmatrix) + for (varid, coef) in @view coefmatrix[constrid, :] + push!(terms, MOI.ScalarAffineTerm(coef, model.moi_varids[varid])) end return MOI.ScalarAffineFunction(terms, 0.0) end @@ -619,6 +620,8 @@ function MOI.empty!(model::Coluna.Optimizer) model.env.varids = CleverDicts.CleverDict{MOI.VariableIndex, VarId}() model.moi_varids = Dict{VarId, MOI.VariableIndex}() model.constrs = Dict{MOI.ConstraintIndex, Constraint}() + model.constrs_on_single_var_to_vars = Dict{MOI.ConstraintIndex, VarId}() + model.constrs_on_single_var_to_names = Dict{MOI.ConstraintIndex, String}() if model.default_optimizer_builder !== nothing set_default_optimizer_builder!(model.inner, model.default_optimizer_builder) end @@ -703,18 +706,21 @@ function MOI.get(optimizer::Optimizer, ::MOI.ObjectiveValue) return getvalue(get_ip_primal_bound(optimizer.result)) end +function MOI.get(optimizer::Optimizer, ::MOI.DualObjectiveValue) + return getvalue(get_lp_dual_bound(optimizer.result)) +end + function MOI.get(optimizer::Optimizer, ::MOI.RelativeGap) return ip_gap(optimizer.result) end -function MOI.get(optimizer::Optimizer, ::MOI.VariablePrimal, ref::MOI.VariableIndex) +function MOI.get(optimizer::Optimizer, attr::MOI.VariablePrimal, ref::MOI.VariableIndex) id = getid(optimizer.vars[ref]) # This gets a coluna Id{Variable} - best_primal_sol = get_best_ip_primal_sol(optimizer.result) - if best_primal_sol === nothing - @warn "Coluna did not find a primal feasible solution." - return NaN + primalsols = get_ip_primal_sols(optimizer.result) + if 1 <= attr.N <= length(primalsols) + return get(primalsols[attr.N], id, 0.0) end - return get(best_primal_sol, id, 0.0) + return error("Invalid result index.") end function MOI.get(optimizer::Optimizer, ::MOI.VariablePrimal, refs::Vector{MOI.VariableIndex}) @@ -775,7 +781,20 @@ end MOI.get(optimizer::Optimizer, ::MOI.NodeCount) = optimizer.env.kpis.node_count MOI.get(optimizer::Optimizer, ::MOI.SolveTime) = optimizer.env.kpis.elapsed_optimization_time -# function MOI.get(optimizer::Optimizer, ::MOI.ConstraintDual, index::MOI.ConstraintIndex) +function MOI.get( + optimizer::Optimizer, attr::MOI.ConstraintDual, + index::MOI.ConstraintIndex{MOI.ScalarAffineFunction{Float64}} +) + dualsols = get_lp_dual_sols(optimizer.result) + if 1 <= attr.N <= length(dualsols) + return get(dualsols[attr.N], getid(optimizer.constrs[index]), 0.0) + end + return error("Invalid result index.") +end + +# function MOI.get( +# optimizer::Optimizer, attr::MOI.ConstraintDual, index::MOI.ConstraintIndex{F,S} +# ) where {F<:MOI.SingleVariable,S} # return 0.0 # end diff --git a/test/MathOptInterface/MOI_wrapper.jl b/test/MathOptInterface/MOI_wrapper.jl index 9534b8766..b63291ef3 100644 --- a/test/MathOptInterface/MOI_wrapper.jl +++ b/test/MathOptInterface/MOI_wrapper.jl @@ -3,7 +3,7 @@ const OPTIMIZER = Coluna.Optimizer() MOI.set(OPTIMIZER, MOI.RawParameter("default_optimizer"), GLPK.Optimizer) -const CONFIG = MOIT.TestConfig(atol=1e-6, rtol=1e-6) +const CONFIG = MOIT.TestConfig(atol=1e-6, rtol=1e-6, infeas_certificates = false) @testset "SolverName" begin @@ -130,8 +130,9 @@ const UNSUPPORTED_TESTS = [ "number_threads", # TODO : support of MOI.NumberOfThreads() "silent", # TODO : support of MOI.Silent() "time_limit_sec", # TODO : support of MOI.TimeLimitSec() - "solve_time", # TODO : support of MOI.SolveTime() - "solve_twice" # TODO : fix + "solve_unbounded_model", # default lower bound 0 + "solve_duplicate_terms_obj", # TODO: support duplicate terms + "solve_duplicate_terms_scalar_affine" # TODO: support duplicate terms ] MathOptInterface.Test.getconstraint @@ -178,6 +179,41 @@ const LP_TESTS = [ "solve_affine_lessthan" ] +const CONSTRAINTDUAL_SINGLEVAR = [ + "solve_with_lowerbound", + "solve_singlevariable_obj", + "solve_constant_obj", + "solve_single_variable_dual_max", + "solve_single_variable_dual_min", + "solve_duplicate_terms_obj", + "solve_blank_obj", + "solve_with_upperbound", + "linear1", + "linear2", + "linear10b", + "linear14" +] + +const MODIFY_DELETE = [ + # BUG + "linear1", # modify + "linear5", # modify + "linear11", # delete + "linear14" # delete +] + +const UNCOVERED_TERMINATION_STATUS = [ + "linear8b", # DUAL_INFEASIBLE or INFEASIBLE_OR_UNBOUNDED required + "linear8c" # DUAL_INFEASIBLE or INFEASIBLE_OR_UNBOUNDED required +] + +const SET_CONSTRAINTSET = [ + # BUG + "linear4", + "linear6", + "linear7" +] + @testset "Unit Basic/MIP" begin MOI.set(OPTIMIZER, MOI.RawParameter("params"), CL.Params(solver = ClA.SolveIpForm())) MOIT.unittest(OPTIMIZER, CONFIG, vcat(UNSUPPORTED_TESTS, LP_TESTS, MIP_TESTS)) @@ -197,13 +233,19 @@ MOI.set(BRIDGED, MOI.RawParameter("params"), CL.Params(solver = ClA.SolveIpForm( ]) end -# @testset "Unit LP" begin -# MOI.set(OPTIMIZER, MOI.RawParameter("params"), CL.Params(solver = ClA.SolveLpForm())) -# MOIT.unittest(OPTIMIZER, CONFIG, vcat(UNSUPPORTED_TESTS, MIP_TESTS, BASIC)) -# end +@testset "Unit LP" begin + MOI.set(BRIDGED, MOI.RawParameter("params"), CL.Params(solver = ClA.SolveLpForm( + update_ip_primal_solution=true, get_dual_solution=true, set_dual_bound=true + ))) + MOIT.unittest(BRIDGED, CONFIG, vcat(UNSUPPORTED_TESTS, MIP_TESTS, BASIC, CONSTRAINTDUAL_SINGLEVAR)) +end -# @testset "Continuous Linear" begin -# MOIT.contlineartest(OPTIMIZER, CONFIG, [ -# "partial_start" # VariablePrimalStart not supported -# ]) -# end +@testset "Continuous Linear" begin + MOIT.contlineartest(BRIDGED, CONFIG, vcat( + CONSTRAINTDUAL_SINGLEVAR, MODIFY_DELETE, UNCOVERED_TERMINATION_STATUS, SET_CONSTRAINTSET, [ + "partial_start", # VariablePrimalStart not supported + "linear1", # TODO: support duplicate terms + "linear10" # BUG: optimize twice changing sense from max to min fails + ] + )) +end From 6178b56cc61daa64839de4c1cdb1a33c9557c565 Mon Sep 17 00:00:00 2001 From: Guillaume Marques Date: Tue, 10 Aug 2021 08:21:38 +0200 Subject: [PATCH 5/6] Don't close fillmode in get terms of constraint (#576) --- .github/workflows/ci.yml | 1 + src/MOIwrapper.jl | 4 +--- 2 files changed, 2 insertions(+), 3 deletions(-) diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index ac21fd3a2..e9ecb9110 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -29,6 +29,7 @@ jobs: - run: julia --project=@. -e 'using Pkg; Pkg.add(PackageSpec(url="https://github.com/atoptima/ColunaDemos.jl.git"));' - run: julia --project=@. -e 'using Pkg; Pkg.add(PackageSpec(url="https://github.com/rafaelmartinelli/KnapsackLib.jl.git"));' - run: julia --project=@. -e 'using Pkg; Pkg.add(PackageSpec(url="https://github.com/atoptima/BlockDecomposition.jl.git"));' + - run: julia --project=@. -e 'using Pkg; Pkg.add(PackageSpec(url="https://github.com/atoptima/DynamicSparseArrays.jl.git"));' - uses: actions/cache@v1 env: cache-name: cache-artifacts diff --git a/src/MOIwrapper.jl b/src/MOIwrapper.jl index 2a0474bd7..37a06fe51 100644 --- a/src/MOIwrapper.jl +++ b/src/MOIwrapper.jl @@ -301,9 +301,7 @@ function MOI.get( orig_form = get_original_formulation(model.inner) constrid = getid(model.constrs[index]) terms = MOI.ScalarAffineTerm{Float64}[] - coefmatrix = getcoefmatrix(orig_form) - coefmatrix.fillmode && closefillmode!(coefmatrix) - for (varid, coef) in @view coefmatrix[constrid, :] + for (varid, coef) in @view getcoefmatrix(orig_form)[constrid, :] push!(terms, MOI.ScalarAffineTerm(coef, model.moi_varids[varid])) end return MOI.ScalarAffineFunction(terms, 0.0) From b4a66a697c6d3dc67c7dccdfc8f3faa530984b4e Mon Sep 17 00:00:00 2001 From: Lara di Cavalcanti Pontes <34425678+laradicp@users.noreply.github.com> Date: Tue, 10 Aug 2021 14:26:38 -0300 Subject: [PATCH 6/6] fix retrieval of sp specific disaggregated solution (#578) --- src/MOIwrapper.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/MOIwrapper.jl b/src/MOIwrapper.jl index 37a06fe51..e351b6bb7 100644 --- a/src/MOIwrapper.jl +++ b/src/MOIwrapper.jl @@ -635,11 +635,13 @@ mutable struct ColumnInfo <: BD.AbstractColumnInfo end function BD.getsolutions(model::Coluna.Optimizer, k) - ip_primal_sol = model.disagg_result.ip_primal_sols[k] + ip_primal_sol = get_best_ip_primal_sol(model.disagg_result) sp_columns_info = Vector{ColumnInfo}() for (varid, val) in ip_primal_sol if getduty(varid) <= MasterCol - push!(sp_columns_info, ColumnInfo(model, varid, val)) + if model.annotations.ann_per_form[getoriginformuid(varid)].axis_index_value == k + push!(sp_columns_info, ColumnInfo(model, varid, val)) + end end end return sp_columns_info