Optimality tolerance depending on the algorithm

src/Algorithm/Algorithm.jl

@@ -20,7 +20,7 @@ import ..MathProg: getterminationstatus,
     setterminationstatus!, get_ip_primal_bound, get_ip_dual_bound,
     get_lp_primal_bound, get_lp_dual_bound, update_ip_primal_bound!, update_ip_dual_bound!,
     update_lp_primal_bound!, update_lp_dual_bound!, set_ip_primal_bound!,
-    set_ip_dual_bound!, set_lp_primal_bound!, set_lp_dual_bound!, ip_gap, lp_gap
+    set_ip_dual_bound!, set_lp_primal_bound!, set_lp_dual_bound!, ip_gap, lp_gap, ip_gap_closed, lp_gap_closed
 
 # Utilities to build algorithms
 include("utilities/optimizationstate.jl")

src/Algorithm/benders.jl

@@ -2,7 +2,7 @@
     option_use_reduced_cost::Bool = false
     option_increase_cost_in_hybrid_phase::Bool = false
     feasibility_tol::Float64 = 1e-5
-    optimality_tol::Float64 = 1e-5
+    optimality_tol::Float64 = Coluna.DEF_OPTIMALITY_ATOL
     max_nb_iterations::Int = 100
 end
 

src/Algorithm/branching/branchingalgo.jl

@@ -233,7 +233,7 @@ function run!(algo::StrongBranching, data::ReformData, input::DivideInput)::Divi
             original_solution = get_best_lp_primal_sol(optstate)
         end
     else
-        @logmsg LogLevel(0) "Warning: no LP solution is passed to the branching algorithm. No children will be generated."
+        @warn "no LP solution is passed to the branching algorithm. No children will be generated."
         return DivideOutput(Vector{Node}(), optstate)
     end
 

src/Algorithm/colgen.jl

@@ -7,18 +7,18 @@
             log_level = 2
         ),
         max_nb_iterations::Int = 1000
-        optimality_tol::Float64 = 1e-5
         log_print_frequency::Int = 1
         store_all_ip_primal_sols::Bool = false
         redcost_tol::Float = 1e-5
         cleanup_threshold::Int = 10000
         cleanup_ratio::Float = 0.66
         smoothing_stabilization::Float64 = 0.0 # should be in [0, 1]
+        opt_atol::Float64 = DEF_OPTIMALITY_ATOL
+        opt_rtol::Float64 = DEF_OPTIMALITY_RTOL
     )
 
 Column generation algorithm. It applies `restr_master_solve_alg` to solve the linear
 restricted master and `pricing_prob_solve_alg` to solve the subproblems.
-
 """
 @with_kw struct ColumnGeneration <: AbstractOptimizationAlgorithm
     restr_master_solve_alg = SolveLpForm(get_dual_solution = true)
@@ -28,14 +28,15 @@ restricted master and `pricing_prob_solve_alg` to solve the subproblems.
                                          enforce_integrality = false,
                                          log_level = 2)
     max_nb_iterations::Int64 = 1000
-    optimality_tol::Float64 = 1e-5
     log_print_frequency::Int64 = 1
     store_all_ip_primal_sols::Bool = false
     redcost_tol::Float64 = 1e-5
     solve_subproblems_parallel::Bool = false
     cleanup_threshold::Int64 = 10000
     cleanup_ratio::Float64 = 0.66
     smoothing_stabilization::Float64 = 0.0 # should be in [0, 1]
+    opt_atol::Float64 = Coluna.DEF_OPTIMALITY_ATOL
+    opt_rtol::Float64 = Coluna.DEF_OPTIMALITY_RTOL
 end
 
 stabilization_is_used(algo::ColumnGeneration) = !iszero(algo.smoothing_stabilization)
@@ -462,8 +463,8 @@ function cleanup_columns(algo::ColumnGeneration, iteration::Int64, data::ReformD
     return
 end
 
-ph_one_infeasible_db(algo, db::DualBound{MinSense}) = getvalue(db) > algo.optimality_tol
-ph_one_infeasible_db(algo, db::DualBound{MaxSense}) = getvalue(db) < - algo.optimality_tol
+ph_one_infeasible_db(algo, db::DualBound{MinSense}) = getvalue(db) > algo.opt_atol
+ph_one_infeasible_db(algo, db::DualBound{MaxSense}) = getvalue(db) < - algo.opt_atol
 
 function update_lagrangian_dual_bound!(
     stabstorage::ColGenStabilizationStorage, optstate::OptimizationState{F, S}, algo::ColumnGeneration,
@@ -632,7 +633,7 @@ function cg_main_loop!(
 
         if phase != 1 && getterminationstatus(rm_optstate) == INFEASIBLE
             @warn string("Solver returned that LP restricted master is infeasible or unbounded ",
-            "(feasibility status = " , status, ") during phase != 1.")
+            "(termination status = INFEASIBLE) during phase != 1.")
             setterminationstatus!(cg_optstate, INFEASIBLE)
             return true
         end
@@ -722,7 +723,7 @@ function cg_main_loop!(
         primal_bound = get_lp_primal_bound(cg_optstate)
         ip_primal_bound = get_ip_primal_bound(cg_optstate)
 
-        if ip_gap(cg_optstate) < algo.optimality_tol
+        if ip_gap_closed(cg_optstate, atol = algo.opt_atol, rtol = algo.opt_rtol)
             setterminationstatus!(cg_optstate, OPTIMAL)
             @logmsg LogLevel(0) "Dual bound reached primal bound."
             return true
@@ -736,11 +737,16 @@ function cg_main_loop!(
             @logmsg LogLevel(0) "Phase one determines infeasibility."
             return true
         end
-        if nb_new_columns == 0 || lp_gap(cg_optstate) < algo.optimality_tol
-            @logmsg LogLevel(0) "Column Generation Algorithm has converged."
+        if lp_gap_closed(cg_optstate, atol = algo.opt_atol, rtol = algo.opt_rtol)
+            @logmsg LogLevel(0) "Column generation algorithm has converged."
             setterminationstatus!(cg_optstate, OPTIMAL)
             return false
         end
+        if nb_new_columns == 0
+            @logmsg LogLevel(0) "No new column generated by the pricing problems."
+            setterminationstatus!(cg_optstate, OTHER_LIMIT)
+            return false
+        end
         if iteration > algo.max_nb_iterations
             setterminationstatus!(cg_optstate, OTHER_LIMIT)
             @warn "Maximum number of column generation iteration is reached."

src/Algorithm/conquer.jl

@@ -52,11 +52,9 @@ function apply_conquer_alg_to_node!(
         @logmsg LogLevel(-1) string("Node IP DB: ", get_ip_dual_bound(nodestate))
         @logmsg LogLevel(-1) string("Tree IP PB: ", get_ip_primal_bound(nodestate))
     end
-    if (ip_gap(nodestate) <= 0.0 + 0.00000001)
-        isverbose(algo) && @logmsg LogLevel(-1) string(
-            "IP Gap is non-positive: ", ip_gap(getincumbents(node)), ". Abort treatment."
-        )
-    else    
+    if ip_gap_closed(nodestate)
+        @info "IP Gap is closed: $(ip_gap(getincumbents(node))). Abort treatment."
+    else
         isverbose(algo) && @logmsg LogLevel(-1) string("IP Gap is positive. Need to treat node.")
 
         run!(algo, data, ConquerInput(node, storages_to_restore))
@@ -120,6 +118,8 @@ to optimize the integer restricted master.
     max_nb_cut_rounds::Int = 3 # TODO : tailing-off ?
     run_mastipheur::Bool = true
     run_preprocessing::Bool = false
+    opt_atol::Float64 = colgen.opt_atol # TODO : force this value in an init() method
+    opt_rtol::Float64 = colgen.opt_rtol # TODO : force this value in an init() method
 end
 
 isverbose(algo::ColCutGenConquer) = algo.colgen.log_print_frequency > 0
@@ -150,7 +150,10 @@ function run!(algo::ColCutGenConquer, data::ReformData, input::ConquerInput)
     colgen_output = run!(algo.colgen, data, OptimizationInput(nodestate))
     update!(nodestate, getoptstate(colgen_output))
 
-    while !to_be_pruned(node) && nb_tightening_rounds < algo.max_nb_cut_rounds
+    node_pruned = getterminationstatus(nodestate) == INFEASIBLE ||
+        ip_gap_closed(nodestate, atol = algo.opt_atol, rtol = algo.opt_rtol)
+
+    while !node_pruned && nb_tightening_rounds < algo.max_nb_cut_rounds
         sol = get_best_lp_primal_sol(getoptstate(colgen_output))
         if sol !== nothing
             cutcb_input = CutCallbacksInput(sol)
@@ -166,18 +169,30 @@ function run!(algo::ColCutGenConquer, data::ReformData, input::ConquerInput)
         colgen_output = run!(algo.colgen, data, OptimizationInput(nodestate))
         update!(nodestate, getoptstate(colgen_output))
 
+        node_pruned = getterminationstatus(nodestate) == INFEASIBLE ||
+            ip_gap_closed(nodestate, atol = algo.opt_atol, rtol = algo.opt_rtol)
+
         nb_tightening_rounds += 1
     end
 
-    if !to_be_pruned(node) && algo.run_mastipheur 
+    if !node_pruned && algo.run_mastipheur 
         @logmsg LogLevel(0) "Run IP restricted master heuristic."
         TO.@timeit Coluna._to "RestMasterHeur" begin
             heur_output = run!(
                 algo.mastipheur, getmasterdata(data), OptimizationInput(nodestate)
             )
             update_all_ip_primal_solutions!(nodestate, getoptstate(heur_output))
         end
-    end 
+        node_pruned = ip_gap_closed(
+            nodestate, atol = algo.opt_atol, rtol = algo.opt_rtol
+        )
+    end
+
+    if node_pruned
+        setterminationstatus!(nodestate, OPTIMAL)
+    else
+        setterminationstatus!(nodestate, OTHER_LIMIT)
+    end
     return
 end
 
@@ -212,6 +227,13 @@ function run!(algo::RestrMasterLPConquer, data::ReformData, input::ConquerInput)
     node = getnode(input)
     nodestate = getoptstate(node)
     output = run!(algo.masterlpalgo, getmasterdata(data), OptimizationInput(nodestate))
-    update!(nodestate, getoptstate(output))
+    masterlp_state =  getoptstate(output)
+    update!(nodestate, masterlp_state)
+    if ip_gap_closed(masterlp_state)
+        setterminationstatus!(nodestate, OPTIMAL)
+    else
+        setterminationstatus!(nodestate, OTHER_LIMIT)
+    end
+    return
 end
 

src/Algorithm/node.jl

@@ -57,9 +57,9 @@ isrootnode(n::Node) = n.tree_order == 1
 getinfeasible(n::Node) = n.infesible
 setinfeasible(n::Node, status::Bool) = n.infeasible = status
 
+# TODO remove
 function to_be_pruned(node::Node)
     nodestate = getoptstate(node)
     getterminationstatus(nodestate) == INFEASIBLE && return true
-    bounds_ratio = get_ip_primal_bound(nodestate) / get_ip_dual_bound(nodestate)
-    return isapprox(bounds_ratio, 1) || ip_gap(nodestate) < 0
+    return ip_gap_closed(nodestate)
 end