From a56f8fb57d68bd154652f8951e2a012c6501b99c Mon Sep 17 00:00:00 2001 From: "Nicolas L. Guidotti" Date: Wed, 14 Oct 2026 16:54:10 +0200 Subject: [PATCH 1/6] Uses the objective lattice to accept or reject new incumbents. Add a small delta to avoid accepting the incumbent over and over when exploring an empty or truncated sub-MIP. Signed-off-by: Nicolas L. Guidotti --- cpp/src/branch_and_bound/branch_and_bound.cpp | 122 +++++++++++------- cpp/src/branch_and_bound/branch_and_bound.hpp | 3 + cpp/src/mip_heuristics/solution/solution.cu | 9 +- cpp/src/mip_heuristics/solve.cu | 10 ++ 4 files changed, 93 insertions(+), 51 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 5880b3f924..098d09fd14 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -111,8 +111,13 @@ bool check_guess(const lp_problem_t& original_lp, matrix_vector_multiply(original_lp.A, 1.0, guess, -1.0, residual); primal_error = vector_norm_inf(residual); bound_error = 0.0; + f_t integrality_error = 0.0; constexpr bool verbose = false; for (i_t j = 0; j < original_lp.num_cols; j++) { + if (var_types[j] == variable_type_t::INTEGER) { + integrality_error = std::max(integrality_error, std::abs(guess[j] - std::round(guess[j]))); + } + // l_j <= x_j infeas means x_j < l_j or l_j - x_j > 0 const f_t low_bound_err = std::max(0.0, original_lp.lower[j] - guess[j]); // x_j <= u_j infeas means u_j < x_j or x_j - u_j > 0 @@ -139,6 +144,13 @@ bool check_guess(const lp_problem_t& original_lp, if (verbose) { settings.log.printf("Solution is feasible\n"); } feasible = true; } + + settings.log.debug_format( + "incumbent: primal_error={:.5g}, bound_error={:.5g}, integrality_error={:.5g}\n", + primal_error, + bound_error, + integrality_error); + return feasible; } @@ -460,6 +472,62 @@ void branch_and_bound_t::update_user_bound(const lp_problem_t +void branch_and_bound_t::snap_to_lattice(const lp_problem_t& lp, + mip_node_t* node_ptr, + f_t leaf_obj) +{ + if (lp.objective_step.has_step()) { + f_t step = lp.objective_step.step_size; + f_t bias = lp.objective_step.bias; + // Round up to next value on the lattice: k * step + bias >= leaf_obj + f_t k = std::ceil((leaf_obj - bias) / step - settings_.integer_tol); + node_ptr->lower_bound = k * step + bias; + } else if (lp.objective_is_integral) { + node_ptr->lower_bound = std::ceil(leaf_obj - settings_.integer_tol); + } +} + +template +f_t branch_and_bound_t::compute_node_cutoff(const lp_problem_t& lp, + f_t upper_bound) const +{ + if (lp.objective_step.has_step()) { + f_t step = lp.objective_step.step_size; + f_t bias = lp.objective_step.bias; + // Any improving feasible solution must have objective <= upper_bound - step. + f_t k = std::floor((upper_bound - bias) / step + settings_.integer_tol); + return (k - 1) * step + bias + settings_.dual_tol; + } else if (lp.objective_is_integral) { + // If the objective is integral, any feasible solution should produce an upper bound that is + // (approximately) integral. We add a small tolerance and floor this value to get an integer, + // we then subtract 1, to stop simplex on problems that cannot improve the primal objective. + return std::floor(upper_bound + settings_.integer_tol) - 1 + settings_.dual_tol; + } + return upper_bound + settings_.dual_tol; +} + +template +bool branch_and_bound_t::accept_new_incumbent(f_t obj) +{ + if (!incumbent_.has_incumbent) return true; + + // A new incumbent must have objective < threshold. + if (original_lp_.objective_step.has_step()) { + return obj < incumbent_.objective - 0.5 * original_lp_.objective_step.step_size; + } else if (original_lp_.objective_is_integral) { + return incumbent_.objective - 0.5; + } else { + return incumbent_.objective - settings_.primal_tol; + } +} + template bool branch_and_bound_t::set_solution_from_heuristics(const std::vector& solution, heuristics_origin_t origin) @@ -483,7 +551,7 @@ bool branch_and_bound_t::set_solution_from_heuristics(const std::vecto feasible_solution_symbol(origin), compute_user_objective(original_lp_, obj)); - if (!incumbent_.has_incumbent || obj < incumbent_.objective) { + if (accept_new_incumbent(obj)) { f_t primal_err; f_t bound_err; i_t num_fractional; @@ -496,9 +564,11 @@ bool branch_and_bound_t::set_solution_from_heuristics(const std::vecto } is_feasible = check_guess( original_lp_, settings_, var_types_, crushed_solution, primal_err, bound_err, num_fractional); + settings_.log.print_format("upper_bound={}, obj={}\n", incumbent_.objective, obj); mutex_original_lp_.unlock(); mutex_upper_.lock(); - if (is_feasible && (!incumbent_.has_incumbent || obj < incumbent_.objective)) { + + if (is_feasible && accept_new_incumbent(obj)) { f_t current_upper_bound = upper_bound_.load(); upper_bound_ = std::min(current_upper_bound, obj); incumbent_.set_incumbent_solution(obj, crushed_solution); @@ -745,8 +815,7 @@ void branch_and_bound_t::repair_heuristic_solutions() repair_solution(edge_norms_, crushed_solution, repaired_obj, repaired_solution); if (is_feasible) { mutex_upper_.lock(); - - if (!incumbent_.has_incumbent || repaired_obj < incumbent_.objective) { + if (accept_new_incumbent(repaired_obj)) { upper_bound_ = std::min(upper_bound_.load(), repaired_obj); incumbent_.set_incumbent_solution(repaired_obj, repaired_solution); report_heuristic(repaired_obj, heuristics_origin_t::HEURISTICS); @@ -920,7 +989,7 @@ void branch_and_bound_t::add_feasible_solution(const lp_problem_t -void branch_and_bound_t::snap_to_lattice(const lp_problem_t& lp, - mip_node_t* node_ptr, - f_t leaf_obj) -{ - if (lp.objective_step.has_step()) { - f_t step = lp.objective_step.step_size; - f_t bias = lp.objective_step.bias; - // Round up to next value on the lattice: k * step + bias >= leaf_obj - f_t k = std::ceil((leaf_obj - bias) / step - settings_.integer_tol); - node_ptr->lower_bound = k * step + bias; - } else if (lp.objective_is_integral) { - node_ptr->lower_bound = std::ceil(leaf_obj - settings_.integer_tol); - } -} - template template std::pair branch_and_bound_t::update_tree_impl( @@ -1581,21 +1628,7 @@ dual_status_t branch_and_bound_t::solve_node_lp( lp_settings.concurrent_halt = &node_concurrent_halt_; } lp_settings.set_log(false); - f_t cutoff = upper_bound_.load(); - if (worker->leaf_problem.objective_step.has_step()) { - f_t step = worker->leaf_problem.objective_step.step_size; - f_t bias = worker->leaf_problem.objective_step.bias; - // Any improving feasible solution must have objective <= cutoff - step. - f_t k = std::floor((cutoff - bias) / step + settings_.integer_tol); - lp_settings.cut_off = (k - 1) * step + bias + settings_.dual_tol; - } else if (worker->leaf_problem.objective_is_integral) { - // If the objective is integral, any feasible solution should produce an upper bound that is - // (approximately) integral. We add a small tolerance and floor this value to get an integer, - // we then subtract 1, to stop simplex on problems that cannot improve the primal objective. - lp_settings.cut_off = std::floor(cutoff + settings_.integer_tol) - 1 + settings_.dual_tol; - } else { - lp_settings.cut_off = cutoff + settings_.dual_tol; - } + lp_settings.cut_off = compute_node_cutoff(worker->leaf_problem, upper_bound_.load()); lp_settings.inside_mip = 2; lp_settings.time_limit = settings_.time_limit - toc(exploration_stats_.start_time); if (lp_settings.time_limit <= 0.0) { return dual_status_t::TIME_LIMIT; } @@ -5378,7 +5411,7 @@ void branch_and_bound_t::deterministic_process_worker_solutions( deterministic_current_horizon_); bool improved = false; - if (!incumbent_.has_incumbent || sol->objective < incumbent_.objective) { + if (accept_new_incumbent(sol->objective)) { upper_bound_ = std::min(upper_bound_.load(), sol->objective); incumbent_.set_incumbent_solution(sol->objective, sol->solution); current_upper = sol->objective; @@ -5541,8 +5574,7 @@ void branch_and_bound_t::deterministic_sort_replay_events( // Process heuristic solution at its correct work unit timestamp position f_t new_upper = std::numeric_limits::infinity(); - - if (!incumbent_.has_incumbent || hsol.objective < incumbent_.objective) { + if (accept_new_incumbent(hsol.objective)) { upper_bound_ = std::min(upper_bound_.load(), hsol.objective); incumbent_.set_incumbent_solution(hsol.objective, hsol.solution); new_upper = hsol.objective; diff --git a/cpp/src/branch_and_bound/branch_and_bound.hpp b/cpp/src/branch_and_bound/branch_and_bound.hpp index db62fa745c..c626c44356 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.hpp +++ b/cpp/src/branch_and_bound/branch_and_bound.hpp @@ -360,6 +360,9 @@ class branch_and_bound_t { mip_node_t* node_ptr, f_t leaf_obj); + f_t compute_node_cutoff(const simplex::lp_problem_t& lp, f_t upper_bound) const; + bool accept_new_incumbent(f_t obj); + // Launch a new best-first worker from a given bfs worker. void launch_bfs_worker(bfs_worker_t* worker); diff --git a/cpp/src/mip_heuristics/solution/solution.cu b/cpp/src/mip_heuristics/solution/solution.cu index 0d411a0ffc..6198b0e4c5 100644 --- a/cpp/src/mip_heuristics/solution/solution.cu +++ b/cpp/src/mip_heuristics/solution/solution.cu @@ -567,13 +567,10 @@ f_t solution_t::compute_max_int_violation() { return thrust::transform_reduce( handle_ptr->get_thrust_policy(), - thrust::make_counting_iterator(0), - thrust::make_counting_iterator(0) + problem_ptr->n_integer_vars, + problem_ptr->integer_indices.begin(), + problem_ptr->integer_indices.end(), cuda::proclaim_return_type([v = view()] __device__(i_t idx) -> f_t { - if (v.problem.variable_types[idx] == var_t::INTEGER) { - return abs(v.assignment[idx] - round(v.assignment[idx])); - } - return 0.; + return abs(v.assignment[idx] - round(v.assignment[idx])); }), 0., thrust::maximum()); diff --git a/cpp/src/mip_heuristics/solve.cu b/cpp/src/mip_heuristics/solve.cu index 1869976be3..9757df74ec 100644 --- a/cpp/src/mip_heuristics/solve.cu +++ b/cpp/src/mip_heuristics/solve.cu @@ -840,9 +840,19 @@ mip_solution_t solve_mip_helper( full_sol.copy_new_assignment( cuopt::host_copy(primal_solution, op_problem.get_handle_ptr()->get_stream())); full_sol.compute_feasibility(); + CUOPT_LOG_DEBUG( + "Post solve: max residual error=%g, max integrality error=%g, max variable bound " + "error=%g", + full_sol.compute_max_constraint_violation(), + full_sol.compute_max_int_violation(), + full_sol.compute_max_variable_violation()); if (!full_sol.get_feasible()) { CUOPT_LOG_WARN("The solution is not feasible after post solve"); } + if (mip::has_variable_bounds_violation( + full_problem.handle_ptr, full_sol.assignment, &full_problem)) { + CUOPT_LOG_WARN("The solution violates the variable bounds after post solve"); + } auto full_stats = sol.get_stats(); // add third party presolve time to cuopt presolve time From fca56c93371f8a7d94ba0145ad2c8ceda1925577 Mon Sep 17 00:00:00 2001 From: "Nicolas L. Guidotti" Date: Wed, 14 Oct 2026 17:12:59 +0200 Subject: [PATCH 2/6] Enabled debug messages Signed-off-by: Nicolas L. Guidotti --- cpp/src/branch_and_bound/branch_and_bound.cpp | 44 ++++++++++++++----- cpp/src/mip_heuristics/solve.cu | 2 +- 2 files changed, 35 insertions(+), 11 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 098d09fd14..a36ec97ffa 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -103,7 +103,8 @@ bool check_guess(const lp_problem_t& original_lp, const std::vector& guess, f_t& primal_error, f_t& bound_error, - i_t& num_fractional) + i_t& num_fractional, + f_t& integrality_error) { bool feasible = false; std::vector residual(original_lp.num_rows); @@ -111,7 +112,7 @@ bool check_guess(const lp_problem_t& original_lp, matrix_vector_multiply(original_lp.A, 1.0, guess, -1.0, residual); primal_error = vector_norm_inf(residual); bound_error = 0.0; - f_t integrality_error = 0.0; + integrality_error = 0.0; constexpr bool verbose = false; for (i_t j = 0; j < original_lp.num_cols; j++) { if (var_types[j] == variable_type_t::INTEGER) { @@ -145,7 +146,7 @@ bool check_guess(const lp_problem_t& original_lp, feasible = true; } - settings.log.debug_format( + settings.log.print_format( "incumbent: primal_error={:.5g}, bound_error={:.5g}, integrality_error={:.5g}\n", primal_error, bound_error, @@ -554,6 +555,7 @@ bool branch_and_bound_t::set_solution_from_heuristics(const std::vecto if (accept_new_incumbent(obj)) { f_t primal_err; f_t bound_err; + f_t integral_err; i_t num_fractional; mutex_original_lp_.lock(); if (crushed_solution.size() != original_lp_.num_cols) { @@ -562,8 +564,14 @@ bool branch_and_bound_t::set_solution_from_heuristics(const std::vecto crush_primal_solution( original_problem_, original_lp_, solution, new_slacks_, crushed_solution); } - is_feasible = check_guess( - original_lp_, settings_, var_types_, crushed_solution, primal_err, bound_err, num_fractional); + is_feasible = check_guess(original_lp_, + settings_, + var_types_, + crushed_solution, + primal_err, + bound_err, + num_fractional, + integral_err); settings_.log.print_format("upper_bound={}, obj={}\n", incumbent_.objective, obj); mutex_original_lp_.unlock(); mutex_upper_.lock(); @@ -639,8 +647,15 @@ void branch_and_bound_t::queue_external_solution_deterministic( f_t primal_err; f_t bound_err; i_t num_fractional; - bool is_feasible = check_guess( - original_lp_, settings_, var_types_, crushed_solution, primal_err, bound_err, num_fractional); + f_t integral_err; + bool is_feasible = check_guess(original_lp_, + settings_, + var_types_, + crushed_solution, + primal_err, + bound_err, + num_fractional, + integral_err); mutex_original_lp_.unlock(); if (!is_feasible) { @@ -767,13 +782,15 @@ bool branch_and_bound_t::repair_solution(const std::vector& edge_ f_t primal_error; f_t bound_error; i_t num_fractional; + f_t integral_err; feasible = check_guess(original_lp_, settings_, var_types_, lp_solution.x, primal_error, bound_error, - num_fractional); + num_fractional, + integral_err); repaired_obj = compute_objective(original_lp_, repaired_solution); constexpr bool verbose = false; if (verbose) { @@ -4066,8 +4083,15 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut f_t primal_err; f_t bound_err; i_t num_fractional; - const bool feasible = check_guess( - original_lp_, settings_, var_types_, crushed_guess, primal_err, bound_err, num_fractional); + f_t integral_err; + const bool feasible = check_guess(original_lp_, + settings_, + var_types_, + crushed_guess, + primal_err, + bound_err, + num_fractional, + integral_err); if (feasible) { const f_t computed_obj = compute_objective(original_lp_, crushed_guess); mutex_upper_.lock(); diff --git a/cpp/src/mip_heuristics/solve.cu b/cpp/src/mip_heuristics/solve.cu index 9757df74ec..0168ea8bdc 100644 --- a/cpp/src/mip_heuristics/solve.cu +++ b/cpp/src/mip_heuristics/solve.cu @@ -840,7 +840,7 @@ mip_solution_t solve_mip_helper( full_sol.copy_new_assignment( cuopt::host_copy(primal_solution, op_problem.get_handle_ptr()->get_stream())); full_sol.compute_feasibility(); - CUOPT_LOG_DEBUG( + CUOPT_LOG_INFO( "Post solve: max residual error=%g, max integrality error=%g, max variable bound " "error=%g", full_sol.compute_max_constraint_violation(), From 1e59d0129e109eaa41e027763a759133613cc3d3 Mon Sep 17 00:00:00 2001 From: "Nicolas L. Guidotti" Date: Wed, 14 Oct 2026 17:48:50 +0200 Subject: [PATCH 3/6] Fixed incorrect check Signed-off-by: Nicolas L. Guidotti --- cpp/src/branch_and_bound/branch_and_bound.cpp | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index a36ec97ffa..bfd979951f 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -523,9 +523,9 @@ bool branch_and_bound_t::accept_new_incumbent(f_t obj) if (original_lp_.objective_step.has_step()) { return obj < incumbent_.objective - 0.5 * original_lp_.objective_step.step_size; } else if (original_lp_.objective_is_integral) { - return incumbent_.objective - 0.5; + return obj < incumbent_.objective - 0.5; } else { - return incumbent_.objective - settings_.primal_tol; + return obj < incumbent_.objective - settings_.primal_tol; } } From ecb4a59c13103d7cb8e62555dde1261dfacf58f5 Mon Sep 17 00:00:00 2001 From: "Nicolas L. Guidotti" Date: Wed, 21 Oct 2026 17:13:12 +0200 Subject: [PATCH 4/6] removed debug logs Signed-off-by: Nicolas L. Guidotti --- cpp/src/branch_and_bound/branch_and_bound.cpp | 7 ------- cpp/src/mip_heuristics/solve.cu | 2 +- 2 files changed, 1 insertion(+), 8 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index bfd979951f..b602959784 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -145,13 +145,6 @@ bool check_guess(const lp_problem_t& original_lp, if (verbose) { settings.log.printf("Solution is feasible\n"); } feasible = true; } - - settings.log.print_format( - "incumbent: primal_error={:.5g}, bound_error={:.5g}, integrality_error={:.5g}\n", - primal_error, - bound_error, - integrality_error); - return feasible; } diff --git a/cpp/src/mip_heuristics/solve.cu b/cpp/src/mip_heuristics/solve.cu index 3df1c83d44..e5613f2484 100644 --- a/cpp/src/mip_heuristics/solve.cu +++ b/cpp/src/mip_heuristics/solve.cu @@ -847,7 +847,7 @@ mip_solution_t solve_mip_helper( full_sol.copy_new_assignment( cuopt::host_copy(primal_solution, op_problem.get_handle_ptr()->get_stream())); full_sol.compute_feasibility(); - CUOPT_LOG_INFO( + CUOPT_LOG_DEBUG( "Post solve: max residual error=%g, max integrality error=%g, max variable bound " "error=%g", full_sol.compute_max_constraint_violation(), From 47f0b6a058bac9eceb3b3fe65b84167c1a013b3f Mon Sep 17 00:00:00 2001 From: "Nicolas L. Guidotti" Date: Thu, 22 Oct 2026 11:43:56 +0200 Subject: [PATCH 5/6] Fixed cutoff calculation Signed-off-by: Nicolas L. Guidotti --- cpp/src/branch_and_bound/branch_and_bound.cpp | 8 ++------ 1 file changed, 2 insertions(+), 6 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index b602959784..35c201b68c 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -495,14 +495,10 @@ f_t branch_and_bound_t::compute_node_cutoff(const lp_problem_t Date: Thu, 22 Oct 2026 13:40:08 +0200 Subject: [PATCH 6/6] remove debug leftover. reject incumbent outside the lattice. Signed-off-by: Nicolas L. Guidotti --- cpp/src/branch_and_bound/branch_and_bound.cpp | 37 +++++++++++-------- 1 file changed, 21 insertions(+), 16 deletions(-) diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 35c201b68c..ed9728cc94 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -508,14 +508,20 @@ bool branch_and_bound_t::accept_new_incumbent(f_t obj) { if (!incumbent_.has_incumbent) return true; - // A new incumbent must have objective < threshold. + f_t threshold = incumbent_.objective - settings_.primal_tol; + f_t step = 0.0; if (original_lp_.objective_step.has_step()) { - return obj < incumbent_.objective - 0.5 * original_lp_.objective_step.step_size; + step = original_lp_.objective_step.step_size; } else if (original_lp_.objective_is_integral) { - return obj < incumbent_.objective - 0.5; - } else { - return obj < incumbent_.objective - settings_.primal_tol; + step = 1.0; + } + if (step > 0.0) { + // Half a step above the greatest lattice value strictly below the incumbent. + const f_t lattice_below = + compute_node_cutoff(original_lp_, incumbent_.objective) - settings_.dual_tol; + threshold = std::min(threshold, lattice_below + 0.5 * step); } + return obj < threshold; } template @@ -561,19 +567,19 @@ bool branch_and_bound_t::set_solution_from_heuristics(const std::vecto bound_err, num_fractional, integral_err); - settings_.log.print_format("upper_bound={}, obj={}\n", incumbent_.objective, obj); mutex_original_lp_.unlock(); - mutex_upper_.lock(); - if (is_feasible && accept_new_incumbent(obj)) { - f_t current_upper_bound = upper_bound_.load(); - upper_bound_ = std::min(current_upper_bound, obj); - incumbent_.set_incumbent_solution(obj, crushed_solution); - if (current_upper_bound > upper_bound_.load()) { - report_heuristic(obj, origin); - success = true; + if (is_feasible) { + std::lock_guard lock(mutex_upper_); + if (accept_new_incumbent(obj)) { + f_t current_upper_bound = upper_bound_.load(); + upper_bound_ = std::min(current_upper_bound, obj); + incumbent_.set_incumbent_solution(obj, crushed_solution); + if (current_upper_bound > upper_bound_.load()) { + report_heuristic(obj, origin); + success = true; + } } - } else { attempt_repair = true; constexpr bool verbose = false; @@ -586,7 +592,6 @@ bool branch_and_bound_t::set_solution_from_heuristics(const std::vecto num_fractional); } } - mutex_upper_.unlock(); } else { settings_.log.debug("Solution objective not better than current upper_bound_. Not accepted.\n"); }