Give access to the dual solution of the master to the user (#594)

* Give dual solution of the master to the user

* wip

* ok

* forgot to push a file

src/Algorithm/basic/cutcallback.jl

@@ -4,8 +4,6 @@ todo
 @with_kw struct CutCallbacks <: AbstractAlgorithm
     call_robust_facultative = true
     call_robust_essential = true
-    #call_nonrobust_facultative = false
-    #call_nonrobust_core = false
     tol = 1e-6
 end
 

src/Algorithm/colgen.jl

@@ -57,7 +57,6 @@ Here are their meanings :
 - `DB` is the dual bound of the master LP at the current iteration
 - `mlp` is the objective value of the master LP at the current iteration
 - `PB` is the objective value of the best primal solution found by Coluna at the current iteration
-
 """
 @with_kw struct ColumnGeneration <: AbstractOptimizationAlgorithm
     restr_master_solve_alg = SolveLpForm(get_dual_solution=true)
@@ -691,6 +690,11 @@ function cg_main_loop!(
             lp_bound = get_lp_primal_bound(rm_optstate) + getvalue(partial_solution)
             set_lp_primal_bound!(cg_optstate, lp_bound)
 
+            dual_rm_sol = get_best_lp_dual_sol(rm_optstate)
+            if dual_rm_sol !== nothing
+                set_lp_dual_sol!(cg_optstate, dual_rm_sol)
+            end
+
             if phase != 1 && !contains(rm_sol, varid -> isanArtificialDuty(getduty(varid)))
                 proj_sol = proj_cols_on_rep(rm_sol, masterform)
                 if isinteger(proj_sol) && isbetter(lp_bound, get_ip_primal_bound(cg_optstate))
@@ -774,8 +778,6 @@ function cg_main_loop!(
         update_stab_after_colgen_iteration!(stabunit)
 
         dual_bound = get_ip_dual_bound(cg_optstate)
-        primal_bound = get_lp_primal_bound(cg_optstate)
-        ip_primal_bound = get_ip_primal_bound(cg_optstate)
 
         if ip_gap_closed(cg_optstate, atol=algo.opt_atol, rtol=algo.opt_rtol)
             setterminationstatus!(cg_optstate, OPTIMAL)
@@ -794,6 +796,7 @@ function cg_main_loop!(
         if lp_gap_closed(cg_optstate, atol=algo.opt_atol, rtol=algo.opt_rtol) && !essential_cuts_separated
             @logmsg LogLevel(0) "Column generation algorithm has converged."
             setterminationstatus!(cg_optstate, OPTIMAL)
+            set_lp_dual_sol!(cg_optstate, get_best_lp_dual_sol(cg_optstate))
             return false, false
         end
         if nb_new_columns == 0 && !essential_cuts_separated

src/Algorithm/treesearch.jl

@@ -133,6 +133,7 @@ end
 function nb_open_nodes(data::TreeSearchRuntimeData)
     return nb_open_nodes(data.primary_tree) + nb_open_nodes(data.secondary_tree)
 end
+
 get_tree_order(data::TreeSearchRuntimeData) = data.tree_order
 getoptstate(data::TreeSearchRuntimeData) = data.optstate
 
@@ -259,11 +260,17 @@ function run_conquer_algorithm!(
     )        
 
     add_ip_primal_sols!(treestate, get_ip_primal_sols(nodestate)...)
-
+
+    # TreeSearchAlgorithm returns the primal LP & the dual solution found at the root node
     best_lp_primal_sol = get_best_lp_primal_sol(nodestate)
     if algo.storelpsolution && isrootnode(node) && best_lp_primal_sol !== nothing
         set_lp_primal_sol!(treestate, best_lp_primal_sol) 
-    end 
+    end
+
+    best_lp_dual_sol = get_best_lp_dual_sol(nodestate)
+    if isrootnode(node) && best_lp_dual_sol !== nothing
+        set_lp_dual_sol!(treestate, best_lp_dual_sol)
+    end
     return
 end
 
@@ -283,7 +290,7 @@ function run_divide_algorithm!(
 
     first_child_with_runconquer = true
     for child in children
-        if (child.conquerwasrun)
+        if child.conquerwasrun
             set_tree_order!(child, tsdata.tree_order)
             tsdata.tree_order += 1
             if first_child_with_runconquer
@@ -301,22 +308,25 @@ end
 function updatedualbound!(data::TreeSearchRuntimeData)
     treestate = getoptstate(data)
     bound_value = getvalue(get_ip_primal_bound(treestate))
-    worst_bound = DualBound{data.Sense}(bound_value)  
-    for (node, priority) in getnodes(data.primary_tree)
+    worst_bound = DualBound{data.Sense}(bound_value)
+    for (node, _) in getnodes(data.primary_tree)
         db = get_ip_dual_bound(getoptstate(node))
         if isbetter(worst_bound, db)
             worst_bound = db
         end
     end
-    for (node, priority) in getnodes(data.secondary_tree)
+
+    for (node, _) in getnodes(data.secondary_tree)
         db = get_ip_dual_bound(getoptstate(node))
         if isbetter(worst_bound, db)
             worst_bound = db
         end
     end
+
     if isbetter(worst_bound, data.worst_db_of_pruned_node)
         worst_bound = data.worst_db_of_pruned_node
     end
+
     set_ip_dual_bound!(treestate, worst_bound)
     return
 end
@@ -334,7 +344,7 @@ function run!(
         # dual bound of the optstate only at the root node.
         run_conquer_algorithm!(algo, env, tsdata, reform, node)
         print_node_in_branching_tree_file(algo, env, tsdata, node)
-               
+
         nodestatus = getterminationstatus(node.optstate)
         if nodestatus == OPTIMAL || nodestatus == INFEASIBLE ||
            ip_gap_closed(node.optstate, rtol = algo.opt_rtol, atol = algo.opt_atol)             

src/Algorithm/utilities/optimizationstate.jl

@@ -243,10 +243,17 @@ function update!(dest_state::OptimizationState, orig_state::OptimizationState)
     update_ip_dual_bound!(dest_state, get_ip_dual_bound(orig_state))
     update_lp_dual_bound!(dest_state, get_lp_dual_bound(orig_state))
     set_lp_primal_bound!(dest_state, get_lp_primal_bound(orig_state))
+
     best_lp_primal_sol = get_best_lp_primal_sol(orig_state) 
     if best_lp_primal_sol !== nothing
         set_lp_primal_sol!(dest_state, best_lp_primal_sol)
-    end        
+    end     
+
+    best_lp_dual_sol = get_best_lp_dual_sol(orig_state)
+    if best_lp_dual_sol !== nothing
+        set_lp_dual_sol!(dest_state, best_lp_dual_sol)
+    end
+    return
 end
 
 """

src/MOIcallbacks.jl

@@ -73,12 +73,12 @@ end
 ############################################################################################
 #  Robust Constraints Callback                                                             #
 ############################################################################################
-function _register_callback!(form::Formulation, attr::MOI.UserCutCallback, sep::Function)
+function _register_callback!(form::Formulation, ::MOI.UserCutCallback, sep::Function)
     set_robust_constr_generator!(form, Facultative, sep)
     return
 end
 
-function _register_callback!(form::Formulation, attr::MOI.LazyConstraintCallback, sep::Function)
+function _register_callback!(form::Formulation, ::MOI.LazyConstraintCallback, sep::Function)
     set_robust_constr_generator!(form, Essential, sep)
     return
 end

src/MOIwrapper.jl

@@ -879,6 +879,18 @@ function MOI.get(
     return error("Invalid result index.")
 end
 
+# Useful method to retrieve dual values of generated cuts because they don't 
+# have MOI.ConstraintIndex
+function MOI.get(
+    optimizer::Optimizer, attr::MOI.ConstraintDual, constrid::ConstrId
+)
+    dualsols = get_lp_dual_sols(optimizer.result)
+    if 1 <= attr.N <= length(dualsols)
+        return get(dualsols[attr.N], constrid, 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}

src/optimize.jl

@@ -43,7 +43,7 @@ function optimize!(env::Env, prob::MathProg.Problem, annotations::Annotations)
 
     # Apply decomposition
     reformulate!(prob, annotations, env)
-    
+
     # Coluna ready to start
     set_optim_start_time!(env)
     @logmsg LogLevel(-1) "Coluna ready to start."
@@ -79,7 +79,7 @@ function optimize!(
 
     algorithm = env.params.solver
 
-    #this will initialize all the units used by the algorithm and its child algorithms
+    # initialize all the units used by the algorithm and its child algorithms
     Algorithm.initialize_storage_units!(reform, algorithm)
 
     output = Algorithm.run!(algorithm, env, reform, Algorithm.OptimizationInput(initstate))
@@ -110,6 +110,16 @@ function optimize!(
         add_lp_primal_sol!(outstate, proj_cols_on_rep(lp_primal_sol, master))
     end
 
+    # lp_dual_sol to retrieve, for instance, the dual value of generated cuts
+    lp_dual_sol = get_best_lp_dual_sol(algstate)
+    if lp_dual_sol !== nothing
+        add_lp_dual_sol!(outstate, lp_dual_sol)
+    end
+
+    # It returns two optimisation states.
+    # The first one contains the solutions projected on the original formulation.
+    # The second one contains the solutions to the master formulation so the user can
+    # retrieve the disagreggated solution.
     return outstate, algstate
 end
 

test/issues_tests.jl

@@ -404,6 +404,53 @@ function simple_benders()
     @test objective_value(model) == 1.0
 end
 
+# See https://github.com/atoptima/Coluna.jl/issues/591
+function get_dual_of_generated_cuts()
+    coluna = JuMP.optimizer_with_attributes(
+        Coluna.Optimizer,
+        "params" => Coluna.Params(
+            solver=Coluna.Algorithm.TreeSearchAlgorithm(),
+        ),
+        "default_optimizer" => GLPK.Optimizer 
+    );
+
+    model = BlockModel(coluna, direct_model=true)
+
+    @axis(I, 1:7)
+
+    @variable(model, 0<= x[i in I] <= 1) # subproblem variables & constraints
+    @variable(model, y[1:2] >= 0) # master
+    @variable(model, u >=0) # master
+
+    @constraint(model, xCon, sum(x[i] for i = I) <= 1)
+    @constraint(model, yCon, sum(y[i] for i = 1:2) == 1)
+    @constraint(model, initCon1, u >= 0.9*y[1] + y[2] - x[1] - x[2] - x[3])
+    @constraint(model, initCon2, u >= y[1] + y[2] - x[7])
+
+    @objective(model, Min, u)
+
+    callback_called = false
+    constrid = nothing
+    function my_callback_function(cbdata)
+        if !callback_called
+            con = @build_constraint(u >= y[1] + 0.9*y[2] - x[5] - x[6])
+            constrid = MOI.submit(model, MOI.LazyConstraint(cbdata), con)
+            callback_called = true
+        end
+        return
+    end
+
+    MOI.set(model, MOI.LazyConstraintCallback(), my_callback_function)
+
+    @dantzig_wolfe_decomposition(model, dec, I)
+
+    optimize!(model)
+
+    @test objective_value(model) ≈ 0.63333333
+    @test MOI.get(JuMP.unsafe_backend(model), MOI.ConstraintDual(), constrid) ≈ 0.33333333
+    return
+end
+
 function test_issues_fixed()
     @testset "no_decomposition" begin
         solve_with_no_decomposition()
@@ -444,4 +491,8 @@ function test_issues_fixed()
     @testset "simple benders decomposition" begin
         simple_benders()
     end
+
+    @testset "dual of generated cuts" begin
+        get_dual_of_generated_cuts()
+    end
 end