Retrieve disaggregated solution (#547)

* first draft without BD * Merge master into sol_disaggregation * fix ColumnInfo and start link with BD * add disagg_result & sps_info to Optimizer and use BD.getsolutions * use BD.value methods * remove sps_info and MOI methods for SpsInfo * improve test Co-authored-by: Guillaume Marques <[email protected]>
atoptima · Jul 1, 2021 · 51db6d4 · 51db6d4
1 parent 68829cb
commit 51db6d4
Show file tree

Hide file tree

Showing 4 changed files with 80 additions and 7 deletions.
diff --git a/src/MOIwrapper.jl b/src/MOIwrapper.jl
@@ -29,6 +29,7 @@ mutable struct Optimizer <: MOI.AbstractOptimizer
     constrs_on_single_var_to_names::Dict{MOI.ConstraintIndex, String}
     names_to_constrs::Dict{String, MOI.ConstraintIndex}
     result::OptimizationState
+    disagg_result::Union{Nothing, OptimizationState}
     default_optimizer_builder::Union{Nothing, Function}
 
     feasibility_sense::Bool # Coluna supports only Max or Min.
@@ -47,6 +48,7 @@ mutable struct Optimizer <: MOI.AbstractOptimizer
         model.constrs_on_single_var_to_names = Dict{MOI.ConstraintIndex, String}()
         model.names_to_constrs = Dict{String, MOI.ConstraintIndex}()
         model.result = OptimizationState(get_optimization_target(model.inner))
+        model.disagg_result = nothing
         model.default_optimizer_builder = nothing
         model.feasibility_sense = false
         return model
@@ -95,7 +97,7 @@ end
 MOI.get(optimizer::Coluna.Optimizer, ::MOI.SolverName) = "Coluna"
 
 function MOI.optimize!(optimizer::Optimizer)
-    optimizer.result = optimize!(
+    optimizer.result, optimizer.disagg_result = optimize!(
         optimizer.env, optimizer.inner, optimizer.annotations
     )
     return
@@ -588,9 +590,36 @@ function MOI.empty!(model::Coluna.Optimizer)
         set_default_optimizer_builder!(model.inner, model.default_optimizer_builder)
     end
     model.result = OptimizationState(get_optimization_target(model.inner))
+    model.disagg_result = nothing
     return
 end
 
+mutable struct ColumnInfo <: BD.AbstractColumnInfo
+    optimizer::Coluna.Optimizer
+    column_var_id::VarId
+    column_val::Float64
+end
+
+function BD.getsolutions(model::Coluna.Optimizer, k)
+    ip_primal_sol = model.disagg_result.ip_primal_sols[k]
+    sp_columns_info = Vector{ColumnInfo}()
+    for (varid, val) in ip_primal_sol
+        if getduty(varid) <= MasterCol
+            push!(sp_columns_info, ColumnInfo(model, varid, val))
+        end
+    end
+    return sp_columns_info
+end
+
+BD.value(info::ColumnInfo) = info.column_val
+
+function BD.value(info::ColumnInfo, index::MOI.VariableIndex)
+    varid = info.optimizer.env.varids[index]
+    origin_form_uid = getoriginformuid(info.column_var_id)
+    spform = get_dw_pricing_sps(info.optimizer.inner.re_formulation)[origin_form_uid]
+    return info.column_val * getprimalsolmatrix(spform)[varid, info.column_var_id]
+end
+
 function MOI.get(model::Coluna.Optimizer, ::MOI.NumberOfVariables)
     orig_form = get_original_formulation(model.inner)
     return length(getvars(orig_form))

diff --git a/src/optimize.jl b/src/optimize.jl
@@ -60,7 +60,7 @@ function optimize!(env::Env, prob::MathProg.Problem, annotations::Annotations)
     @logmsg LogLevel(-1) env.params
 
     TO.@timeit _to "Coluna" begin
-        optstate = optimize!(get_optimization_target(prob), env, init_pb, init_db)
+        outstate, algstate = optimize!(get_optimization_target(prob), env, init_pb, init_db)
     end
 
     env.kpis.elapsed_optimization_time = elapsed_optim_time(env)
@@ -71,9 +71,9 @@ function optimize!(env::Env, prob::MathProg.Problem, annotations::Annotations)
     TO.reset_timer!(_to)
 
     @logmsg LogLevel(0) "Terminated"
-    @logmsg LogLevel(0) string("Primal bound: ", get_ip_primal_bound(optstate))
-    @logmsg LogLevel(0) string("Dual bound: ", get_ip_dual_bound(optstate))
-    return optstate
+    @logmsg LogLevel(0) string("Primal bound: ", get_ip_primal_bound(outstate))
+    @logmsg LogLevel(0) string("Dual bound: ", get_ip_dual_bound(outstate))
+    return outstate, algstate
 end
 
 function optimize!(
@@ -120,7 +120,7 @@ function optimize!(
         add_lp_primal_sol!(outstate, proj_cols_on_rep(lp_primal_sol, master))
     end
 
-    return outstate
+    return outstate, algstate
 end
 
 function optimize!(
@@ -134,7 +134,7 @@ function optimize!(
     )
     algorithm = env.params.solver
     output = Algorithm.run!(algorithm, env, form, Algorithm.OptimizationInput(initstate))
-    return Algorithm.getoptstate(output)
+    return Algorithm.getoptstate(output), nothing
 end
 
 """

diff --git a/test/runtests.jl b/test/runtests.jl
@@ -35,6 +35,7 @@ include("bound_callback_tests.jl")
 include("optimizer_with_attributes_test.jl")
 include("subproblem_solvers_tests.jl")
 include("custom_var_cuts_tests.jl")
+include("sol_disaggregation_tests.jl")
 
 rng = MersenneTwister(1234123)
 
@@ -80,3 +81,7 @@ end
 @testset "Custom Variables and Cuts" begin
     custom_var_cuts_test()
 end
+
+@testset "Solution Disaggregation" begin
+    sol_disaggregation_tests()
+end
diff --git a/test/sol_disaggregation_tests.jl b/test/sol_disaggregation_tests.jl
@@ -0,0 +1,39 @@
+function sol_disaggregation_tests()
+    I = 1:3
+    @axis(BinsType, [1])
+
+    w = [2, 5, 7]
+    Q = 8
+
+    coluna = JuMP.optimizer_with_attributes(
+        Coluna.Optimizer,
+        "params" => Coluna.Params(solver = ColumnGeneration()),
+        "default_optimizer" => GLPK.Optimizer
+    )
+
+    model = BlockModel(coluna)
+
+    @variable(model, x[k in BinsType, i in I], Bin)
+    @variable(model, y[k in BinsType], Bin)
+
+    @constraint(model, sp[i in I], sum(x[k, i] for k in BinsType) == 1)
+    @constraint(model, ks[k in BinsType], sum(w[i] * x[k, i] for i in I) - y[k] * Q <= 0)
+
+    @objective(model, Min, sum(y[k] for k in BinsType))
+
+    @dantzig_wolfe_decomposition(model, dec, BinsType)
+    subproblems = BlockDecomposition.getsubproblems(dec)
+    specify!.(subproblems, lower_multiplicity = 0, upper_multiplicity = BD.length(I)) # we use at most 3 bins
+
+    JuMP.optimize!(model)
+
+    for k in BinsType
+        bins = BD.getsolutions(model, k)
+        for bin in bins
+            @test BD.value(bin) == 1.0 # value of the master column variable
+            @test BD.value(bin, x[k, 1]) == BD.value(bin, x[k, 2]) # x[1,1] and x[1,2] in the same bin
+            @test BD.value(bin, x[k, 1]) != BD.value(bin, x[k, 3]) # only x[1,3] in its bin
+            @test BD.value(bin, x[k, 2]) != BD.value(bin, x[k, 3]) # only x[1,3] in its bin
+        end
+    end
+end