Test Benders loop on toy example

src/Algorithm/benders/default.jl

@@ -276,6 +276,7 @@ struct BendersOutput <: Benders.AbstractBendersOutput
     infeasible_master::Bool
     infeasible_subproblem::Bool
     time_limit_reached::Bool
+    mlp::Float64
 end
 
 Benders.benders_output_type(::BendersContext) = BendersOutput
@@ -287,7 +288,8 @@ function Benders.new_output(
     return BendersOutput(
         benders_iter_output.infeasible_master,
         benders_iter_output.infeasible_subproblem,
-        benders_iter_output.time_limit_reached
+        benders_iter_output.time_limit_reached,
+        benders_iter_output.master
     )
 end
 

test/parser.jl

@@ -365,6 +365,7 @@ module Parser
                         duty = var.duty
                         if var.duty <= ClMP.MasterBendFirstStageVar
                             duty = ClMP.BendSpFirstStageRepVar
+                            explicit = false
                         end
                         if var.duty <= ClMP.MasterBendSecondStageCostVar
                             duty = ClMP.BendSpCostRepVar
@@ -395,33 +396,38 @@ module Parser
         return subproblems, all_spvars, constraints
     end
 
-    _sp_duty_data!(sp::ClMP.Formulation{ClMP.DwSp}, var, duty) = nothing
-    function _sp_duty_data!(sp::ClMP.Formulation{ClMP.BendersSp}, var, duty)
+    _sp_duty_data!(sp::ClMP.Formulation{ClMP.DwSp}, sp_var, master_var, duty) = nothing
+    function _sp_duty_data!(sp::ClMP.Formulation{ClMP.BendersSp}, sp_var, master_var, duty)
         if duty <= ClMP.MasterBendSecondStageCostVar
-            sp.duty_data.second_stage_cost_var = ClMP.getid(var)
+            sp.duty_data.second_stage_cost_var = ClMP.getid(master_var)
         end
     end
 
-    function add_master_vars!(master::ClMP.Formulation, all_spvars::Dict, cache::ReadCache)
+    function add_master_vars!(master::ClMP.Formulation, master_duty, all_spvars::Dict, cache::ReadCache)
         mastervars = Dict{String, ClMP.Variable}()
         for (varid, cost) in cache.master.objective.vars
             if haskey(cache.variables, varid)
                 var = cache.variables[varid]
-                if var.duty <= ClMP.AbstractOriginMasterVar || var.duty <= ClMP.AbstractAddedMasterVar
+                if (typeof(master_duty) <: ClMP.DwMaster && var.duty <= ClMP.AbstractOriginMasterVar) || var.duty <= ClMP.AbstractAddedMasterVar
                     is_explicit = !(var.duty <= ClMP.AbstractImplicitMasterVar)
                     v = ClMP.setvar!(master, varid, var.duty; lb = var.lb, ub = var.ub, kind = var.kind, is_explicit = is_explicit)
                     if haskey(all_spvars, varid)
-                        spvar, sp = all_spvars[varid]
-                        _sp_duty_data!(sp, spvar, var.duty)
+                        var, sp = all_spvars[varid]
+                        _sp_duty_data!(sp, var, v, ClMP.getduty(ClMP.getid(v)))
                     end
                 else
                     if haskey(all_spvars, varid)
                         var, sp = all_spvars[varid]
                         duty = ClMP.MasterRepPricingVar
+                        explicit = false
                         if ClMP.getduty(ClMP.getid(var)) <= ClMP.DwSpSetupVar
                             duty = ClMP.MasterRepPricingSetupVar
+                        elseif ClMP.getduty(ClMP.getid(var)) <= ClMP.BendSpFirstStageRepVar
+                            duty = ClMP.MasterPureVar
+                            explicit = true
                         end
-                        v = ClMP.clonevar!(sp, master, sp, var, duty; is_explicit = false)
+                        v = ClMP.clonevar!(sp, master, sp, var, duty; cost = ClMP.getcurcost(sp, var), is_explicit = explicit)
+                        _sp_duty_data!(sp, var, v, ClMP.getduty(ClMP.getid(var)))
                     else
                         throw(UndefVarParserError("Variable $varid not present in any subproblem"))
                     end
@@ -496,7 +502,7 @@ module Parser
             parent_formulation = reform
         )
         ClMP.setmaster!(reform, master)
-        mastervars = add_master_vars!(master, all_spvars, cache)
+        mastervars = add_master_vars!(master, master_duty, all_spvars, cache)
         ClMP.setobjconst!(master, cache.master.objective.constant)
         add_master_constraints!(reform, master, mastervars, constraints, cache)
 

test/unit/Benders/benders_default.jl

@@ -1,72 +1,115 @@
-function benders_simple_example()
-    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
+# function benders_simple_example()
+#     env = Coluna.Env{Coluna.MathProg.VarId}(Coluna.Params())
+#     origform = Coluna.MathProg.create_formulation!(env, Coluna.MathProg.Original())
 
-    # 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)),
-    ]
-    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
+#     # 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
 
-    @show origform
-
-    # 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"]))
-
-    problem = Coluna.MathProg.Problem(env)
-    Coluna.MathProg.set_original_formulation!(problem, origform)
-
-    Coluna.reformulate!(problem, ann, env)
-    reform = Coluna.MathProg.get_reformulation(problem)
-    return env, reform
-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)),
+#     ]
+#     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
+
+#     @show origform
+
+#     # 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"]))
+
+#     problem = Coluna.MathProg.Problem(env)
+#     Coluna.MathProg.set_original_formulation!(problem, origform)
+
+#     Coluna.reformulate!(problem, ann, env)
+#     reform = Coluna.MathProg.get_reformulation(problem)
+#     return env, reform
+# end
 
-function benders_iteration_default()
+
+function benders_form_A()
     # 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)
+    # optimize!(m)
+    # objective_value(m)
+    # value.(x)
+    # value.(y)
+
+    form = """
+    master
+        min
+        x1 + 4x2 + z
+        s.t.
+        x1 + x2 >= 0
+
+    benders_sp
+        min
+        0x1 + 0x2 + 2y1 + 3y2 + z
+        s.t.
+        -x1 + 4x2 + 2y1 + 3y2 >= 2 {BendTechConstr}
+        x1 + 3x2 + y1 + y2 >= 3 {BendTechConstr}
+        y1 + y2 >= 0
+
+    integer
+        first_stage
+            x1, x2
 
-    env, reform = benders_simple_example()
+    continuous
+        second_stage_cost
+            z
+        second_stage
+            y1, y2
+
+    bounds
+        -Inf <= z <= Inf
+        x1 >= 0
+        x2 >= 0
+        y1 >= 0
+        y2 >= 0
+    """
+    env, _, _, _, reform = reformfromstring(form)
+    return env, reform
+end
+
+# A with continuous first stage finds optimal solution
+function benders_iteration_default_1()
+    #env, reform = benders_simple_example()
+    env, reform = benders_form_A()
 
     master = Coluna.MathProg.getmaster(reform)
     master.optimizers = Coluna.MathProg.AbstractOptimizer[] # dirty
@@ -83,25 +126,39 @@ function benders_iteration_default()
     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)
+    result = Coluna.Benders.run_benders_loop!(ctx, env)
+    @test result.mlp ≈ 3.7142857142857144
+end
+register!(unit_tests, "benders_default", benders_iteration_default_1)
+
+# A with integer first stage finds optimal solution
+function benders_iteration_default_2()
+    #env, reform = benders_simple_example()
+    env, reform = benders_form_A()
 
+    master = Coluna.MathProg.getmaster(reform)
+    master.optimizers = Coluna.MathProg.AbstractOptimizer[] # dirty
+    ClMP.push_optimizer!(master, () -> ClA.MoiOptimizer(GLPK.Optimizer()))
+    for (sp_id, sp) in Coluna.MathProg.get_benders_sep_sps(reform)
+        sp.optimizers = Coluna.MathProg.AbstractOptimizer[] # dirty
+        ClMP.push_optimizer!(sp, () -> ClA.MoiOptimizer(GLPK.Optimizer()))
+    end
 
-    Coluna.Benders.run_benders_loop!(ctx, env)
+    alg = Coluna.Algorithm.BendersCutGeneration(
+        max_nb_iterations = 10,
+        separation_solve_alg = Coluna.Algorithm.SolveIpForm()
+    )
+    ctx = Coluna.Algorithm.BendersPrinterContext(
+        reform, alg;
+        print = true
+    )
+    Coluna.set_optim_start_time!(env)
+
+    result = Coluna.Benders.run_benders_loop!(ctx, env)
+    @test result.mlp ≈ 3.7142857142857144
 end
-register!(unit_tests, "benders_default", benders_iteration_default)
+register!(unit_tests, "benders_default", benders_iteration_default_2; x = true)
+