From 356e9de82caad82630d4730bb6554f3cfc71d971 Mon Sep 17 00:00:00 2001 From: Guillaume Marques Date: Wed, 19 Apr 2023 12:29:12 +0200 Subject: [PATCH] Benders reformulation --- src/decomposition.jl | 72 +------- test/unit/MathProg/benders_decomposition.jl | 181 ++++++++++++++++++++ 2 files changed, 187 insertions(+), 66 deletions(-) create mode 100644 test/unit/MathProg/benders_decomposition.jl diff --git a/src/decomposition.jl b/src/decomposition.jl index 32aefb85f..99f8806ef 100644 --- a/src/decomposition.jl +++ b/src/decomposition.jl @@ -291,26 +291,15 @@ function create_side_vars_constrs!( ::Env, ::Annotations ) - for (_, spform) in get_benders_sep_sps(masterform.parent_formulation) - nu_var = collect(values(filter( - v -> getduty(v.first) == BendSpSlackSecondStageCostVar, - getvars(spform) - )))[1] - - name = "η[$(split(getname(spform, nu_var), "[")[end])" - nu_id = VarId( - getid(nu_var); - duty = MasterBendSecondStageCostVar, - assigned_form_uid = getuid(masterform) - ) + for (spid, spform) in get_benders_sep_sps(masterform.parent_formulation) + name = "η[$(spid)]" setvar!( masterform, name, MasterBendSecondStageCostVar; cost = 1.0, - lb = getperenlb(spform, nu_var), - ub = getperenub(spform, nu_var), + lb = -Inf, + ub = Inf, kind = Continuous, - is_explicit = true, - id = nu_id + is_explicit = true ) end return @@ -327,7 +316,7 @@ function instantiate_orig_vars!( 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) + clonevar!(origform, spform, spform, var, BendSpSepVar, cost = getperencost(origform, var)) end end return @@ -416,55 +405,6 @@ function create_side_vars_constrs!( end end end - - sp_has_second_stage_cost = false - global_costprofit_ub = 0.0 - global_costprofit_lb = 0.0 - for (varid, _) in getvars(spform) - getduty(varid) == BendSpSepVar || continue - orig_var = getvar(origform, varid) - cost = getperencost(origform, orig_var) - if cost > 0 - global_costprofit_ub += cost * getcurub(origform, orig_var) - global_costprofit_lb += cost * getcurlb(origform, orig_var) - 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 || global_costprofit_lb < 0 - sp_has_second_stage_cost = true - end - - if sp_has_second_stage_cost - sp_id = getuid(spform) - # Cost constraint - 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_pos, 0.0) - setcurub!(spform, nu_pos, Inf) - - cost = setconstr!( - spform, "cost[$sp_id]", BendSpSecondStageCostConstr; - rhs = 0.0, - kind = Essential, - sense = getobjsense(spform) == MinSense ? Greater : Less, - is_explicit = true - ) - spcoef[getid(cost), getid(nu_pos)] = 1.0 - - for (varid, _) in getvars(spform) - getduty(varid) == BendSpSepVar || continue - spcoef[getid(cost), varid] = - getperencost(origform, varid) - end - end return end diff --git a/test/unit/MathProg/benders_decomposition.jl b/test/unit/MathProg/benders_decomposition.jl new file mode 100644 index 000000000..c14f12a29 --- /dev/null +++ b/test/unit/MathProg/benders_decomposition.jl @@ -0,0 +1,181 @@ +function benders_decomposition() + """ + original + min + x1 + 4x2 + 2y1 + 3y2 + s.t. + x1 + x2 >= 0 + - x1 + 3x2 - y1 + 2y2 >= 2 + x1 + 3x2 + y1 + y2 >= 3 + y1 + y2 >= 0 + + continuous + original + y1, y2 + + integer + original + x1, x2 + + bounds + x1 >= 0 + x2 >= 0 + y1 >= 0 + y2 >= 0 + """ + + env = Coluna.Env{Coluna.MathProg.VarId}(Coluna.Params()) + + origform = Coluna.MathProg.create_formulation!( + env, Coluna.MathProg.Original() + ) + + # Variables + vars = Dict{String, Coluna.MathProg.VarId}() + variables_infos = [ + ("x1", 1.0, Coluna.MathProg.Integ), + ("x2", 4.0, Coluna.MathProg.Integ), + ("y1", 2.0, Coluna.MathProg.Continuous), + ("y2", 3.0, Coluna.MathProg.Continuous) + ] + for (name, cost, kind) in variables_infos + vars[name] = Coluna.MathProg.getid(Coluna.MathProg.setvar!( + origform, name, Coluna.MathProg.OriginalVar; cost = cost, lb = 0.0, kind = kind + )) + end + + # Constraints + constrs = Dict{String, Coluna.MathProg.ConstrId}() + constraints_infos = [ + ("c1", 2.0, Coluna.MathProg.Greater, Dict(vars["x1"] => -1.0, vars["x2"] => 4.0, vars["y1"] => 2.0, vars["y2"] => 3.0)), + ("c2", 3.0, Coluna.MathProg.Greater, Dict(vars["x1"] => 1.0, vars["x2"] => 3.0, vars["y1"] => 1.0, vars["y2"] => 1.0)), + ("c3", 0.0, Coluna.MathProg.Greater, Dict(vars["x1"] => 1.0, vars["x2"] => 1.0)), + ("c4", 0.0, Coluna.MathProg.Greater, Dict(vars["y1"] => 1.0, vars["y2"] => 1.0)) + ] + for (name, rhs, sense, members) in constraints_infos + constrs[name] = Coluna.MathProg.getid(Coluna.MathProg.setconstr!( + origform, name, Coluna.MathProg.OriginalConstr; rhs = rhs, sense = sense, members = members + )) + end + + # Decomposition tree + m = JuMP.Model() + BlockDecomposition.@axis(axis, [1]) + tree = BlockDecomposition.Tree(m, BlockDecomposition.Benders, axis) + mast_ann = tree.root.master + sp_ann = BlockDecomposition.Annotation(tree, BlockDecomposition.BendersSepSp, BlockDecomposition.Benders, []) + BlockDecomposition.create_leaf!(BlockDecomposition.getroot(tree), axis[1], sp_ann) + + + # Benders annotations + ann = Coluna.Annotations() + ann.tree = tree + Coluna.store!(ann, mast_ann, Coluna.MathProg.getvar(origform, vars["x1"])) + Coluna.store!(ann, mast_ann, Coluna.MathProg.getvar(origform, vars["x2"])) + Coluna.store!(ann, sp_ann, Coluna.MathProg.getvar(origform, vars["y1"])) + Coluna.store!(ann, sp_ann, Coluna.MathProg.getvar(origform, vars["y2"])) + Coluna.store!(ann, sp_ann, Coluna.MathProg.getconstr(origform, constrs["c1"])) + Coluna.store!(ann, sp_ann, Coluna.MathProg.getconstr(origform, constrs["c2"])) + Coluna.store!(ann, mast_ann, Coluna.MathProg.getconstr(origform, constrs["c3"])) + Coluna.store!(ann, sp_ann, Coluna.MathProg.getconstr(origform, constrs["c4"])) + + problem = Coluna.MathProg.Problem(env) + Coluna.MathProg.set_original_formulation!(problem, origform) + + Coluna.reformulate!(problem, ann, env) + reform = Coluna.MathProg.get_reformulation(problem) + + # Test first stage variables & constraints + # Coluna.MathProg.MinSense + 1.0 x1 + 4.0 x2 + 1.0 η[5] + # c3 : + 1.0 x1 + 1.0 x2 >= 0.0 (MasterPureConstrConstraintu3 | true) + # 0.0 <= x1 <= Inf (Continuous | MasterPureVar | true) + # 0.0 <= x2 <= Inf (Continuous | MasterPureVar | true) + # 0.0 <= η[5] <= Inf (Continuous | MasterBendSecondStageCostVar | true) + + master = Coluna.MathProg.getmaster(reform) + fs_vars = Dict(getname(master, varid) => var for (varid, var) in Coluna.MathProg.getvars(master)) + fs_constrs = Dict(getname(master, constrid) => constr for (constrid, constr) in Coluna.MathProg.getconstrs(master)) + + @test length(fs_vars) == 3 + + @test Coluna.MathProg.getduty(Coluna.MathProg.getid(fs_vars["x1"])) <= Coluna.MathProg.MasterPureVar + @test Coluna.MathProg.getduty(Coluna.MathProg.getid(fs_vars["x2"])) <= Coluna.MathProg.MasterPureVar + @test Coluna.MathProg.getduty(Coluna.MathProg.getid(fs_vars["η[5]"])) <= Coluna.MathProg.MasterBendSecondStageCostVar + + @test Coluna.MathProg.getcurlb(master, fs_vars["x1"]) == 0.0 + @test Coluna.MathProg.getcurlb(master, fs_vars["x2"]) == 0.0 + @test Coluna.MathProg.getcurlb(master, fs_vars["η[5]"]) == -Inf + + @test Coluna.MathProg.getcurub(master, fs_vars["x1"]) == Inf + @test Coluna.MathProg.getcurub(master, fs_vars["x2"]) == Inf + @test Coluna.MathProg.getcurub(master, fs_vars["η[5]"]) == Inf + + @test Coluna.MathProg.getcurcost(master, fs_vars["x1"]) == 1.0 + @test Coluna.MathProg.getcurcost(master, fs_vars["x2"]) == 4.0 + @test Coluna.MathProg.getcurcost(master, fs_vars["η[5]"]) == 1.0 + + @test length(fs_constrs) == 1 + + @test Coluna.MathProg.getduty(Coluna.MathProg.getid(fs_constrs["c3"])) <= Coluna.MathProg.MasterPureConstr + + @test Coluna.MathProg.getcurrhs(master, fs_constrs["c3"]) == 0.0 + + # Test second stage variables & Constraints + # Coluna.MathProg.MinSense + 2.0 y1 + 3.0 y2 + 1.0 μ⁺[x1] + 1.0 μ⁻[x1] + 4.0 μ⁺[x2] + 4.0 μ⁻[x2] + # c1 : - 1.0 μ⁺[x1] + 4.0 μ⁺[x2] + 2.0 y1 + 3.0 y2 - 4.0 μ⁻[x2] + 1.0 μ⁻[x1] >= 2.0 (BendSpTechnologicalConstrConstraintu1 | true) + # c2 : + 1.0 μ⁺[x1] + 3.0 μ⁺[x2] + 1.0 y1 + 1.0 y2 - 3.0 μ⁻[x2] - 1.0 μ⁻[x1] >= 3.0 (BendSpTechnologicalConstrConstraintu2 | true) + # c4 : + 1.0 y1 + 1.0 y2 >= 0.0 (BendSpPureConstrConstraintu4 | true) + # 0.0 <= y1 <= Inf (Continuous | BendSpSepVar | true) + # 0.0 <= y2 <= Inf (Continuous | BendSpSepVar | true) + # 0.0 <= μ⁺[x1] <= Inf (Continuous | BendSpPosSlackFirstStageVar | true) + # 0.0 <= μ⁺[x2] <= Inf (Continuous | BendSpPosSlackFirstStageVar | true) + # 0.0 <= μ⁻[x1] <= Inf (Continuous | BendSpNegSlackFirstStageVar | true) + # 0.0 <= μ⁻[x2] <= Inf (Continuous | BendSpNegSlackFirstStageVar | true) + + subprob = first(values(Coluna.MathProg.get_benders_sep_sps(reform))) + ss_vars = Dict(getname(subprob, varid) => var for (varid, var) in Coluna.MathProg.getvars(subprob)) + ss_constrs = Dict(getname(subprob, constrid) => constr for (constrid, constr) in Coluna.MathProg.getconstrs(subprob)) + + @test length(ss_vars) == 6 + + @test Coluna.MathProg.getduty(Coluna.MathProg.getid(ss_vars["y1"])) <= Coluna.MathProg.BendSpSepVar + @test Coluna.MathProg.getduty(Coluna.MathProg.getid(ss_vars["y2"])) <= Coluna.MathProg.BendSpSepVar + @test Coluna.MathProg.getduty(Coluna.MathProg.getid(ss_vars["μ⁺[x1]"])) <= Coluna.MathProg.BendSpPosSlackFirstStageVar + @test Coluna.MathProg.getduty(Coluna.MathProg.getid(ss_vars["μ⁺[x2]"])) <= Coluna.MathProg.BendSpPosSlackFirstStageVar + @test Coluna.MathProg.getduty(Coluna.MathProg.getid(ss_vars["μ⁻[x1]"])) <= Coluna.MathProg.BendSpNegSlackFirstStageVar + @test Coluna.MathProg.getduty(Coluna.MathProg.getid(ss_vars["μ⁻[x2]"])) <= Coluna.MathProg.BendSpNegSlackFirstStageVar + + @test Coluna.MathProg.getcurlb(subprob, ss_vars["y1"]) == 0.0 + @test Coluna.MathProg.getcurlb(subprob, ss_vars["y2"]) == 0.0 + @test Coluna.MathProg.getcurlb(subprob, ss_vars["μ⁺[x1]"]) == 0.0 + @test Coluna.MathProg.getcurlb(subprob, ss_vars["μ⁺[x2]"]) == 0.0 + @test Coluna.MathProg.getcurlb(subprob, ss_vars["μ⁻[x1]"]) == 0.0 + @test Coluna.MathProg.getcurlb(subprob, ss_vars["μ⁻[x2]"]) == 0.0 + + @test Coluna.MathProg.getcurub(subprob, ss_vars["y1"]) == Inf + @test Coluna.MathProg.getcurub(subprob, ss_vars["y2"]) == Inf + @test Coluna.MathProg.getcurub(subprob, ss_vars["μ⁺[x1]"]) == Inf + @test Coluna.MathProg.getcurub(subprob, ss_vars["μ⁺[x2]"]) == Inf + @test Coluna.MathProg.getcurub(subprob, ss_vars["μ⁻[x1]"]) == Inf + @test Coluna.MathProg.getcurub(subprob, ss_vars["μ⁻[x2]"]) == Inf + + @test Coluna.MathProg.getcurcost(subprob, ss_vars["y1"]) == 2.0 + @test Coluna.MathProg.getcurcost(subprob, ss_vars["y2"]) == 3.0 + @test Coluna.MathProg.getcurcost(subprob, ss_vars["μ⁺[x1]"]) == 1.0 + @test Coluna.MathProg.getcurcost(subprob, ss_vars["μ⁺[x2]"]) == 4.0 + @test Coluna.MathProg.getcurcost(subprob, ss_vars["μ⁻[x1]"]) == 1.0 + @test Coluna.MathProg.getcurcost(subprob, ss_vars["μ⁻[x2]"]) == 4.0 + + @test Coluna.MathProg.getduty(Coluna.MathProg.getid(ss_constrs["c1"])) <= Coluna.MathProg.BendSpTechnologicalConstr + @test Coluna.MathProg.getduty(Coluna.MathProg.getid(ss_constrs["c2"])) <= Coluna.MathProg.BendSpTechnologicalConstr + @test Coluna.MathProg.getduty(Coluna.MathProg.getid(ss_constrs["c4"])) <= Coluna.MathProg.BendSpPureConstr + + @test Coluna.MathProg.getcurrhs(subprob, ss_constrs["c1"]) == 2.0 + @test Coluna.MathProg.getcurrhs(subprob, ss_constrs["c2"]) == 3.0 + @test Coluna.MathProg.getcurrhs(subprob, ss_constrs["c4"]) == 0.0 + + @test length(ss_constrs) == 3 + return +end +register!(unit_tests, "benders_decomposition", benders_decomposition) +