diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index 6ce9a4f4d0..e1ff6c15f4 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -70,6 +70,29 @@ i_t fractional_variables(const simplex_solver_settings_t& settings, return fractional.size(); } +template +i_t prune_fixed_fractional_variables(const std::vector& lower_bounds, + const std::vector& upper_bounds, + const simplex_solver_settings_t& settings, + std::vector& fractional) +{ + std::vector new_fractional; + new_fractional.reserve(fractional.size()); + + i_t num_fixed = 0; + for (i_t k = 0; k < (i_t)fractional.size(); k++) { + const i_t j = fractional[k]; + if (std::abs(upper_bounds[j] - lower_bounds[j]) < settings.fixed_tol) { + num_fixed++; + } else { + new_fractional.push_back(j); + } + } + + fractional = std::move(new_fractional); + return num_fixed; +} + template void full_variable_types(const user_problem_t& original_problem, const lp_problem_t& original_lp, @@ -2363,6 +2386,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut root_objective_, root_vstatus_, edge_norms_, + upper_bound_.load(), pc_); } @@ -2372,6 +2396,168 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut return solver_status_; } + // Exploit infeasible/fathomed branches from strong branching for bounds tightening. + // If branching down on x_j is infeasible, we can tighten lb[j] = ceil(x_j*). + // If branching up on x_j is infeasible, we can tighten ub[j] = floor(x_j*). + // With an incumbent, branches whose objective exceeds the cutoff yield the same deductions. + { + const f_t current_upper = upper_bound_.load(); + i_t num_tightened = 0; + i_t num_infeasible = 0; + i_t num_cutoff = 0; + for (i_t k = 0; k < (i_t)fractional.size(); k++) { + const i_t j = fractional[k]; + const f_t sb_down = pc_.strong_branch_down[k]; + const f_t sb_up = pc_.strong_branch_up[k]; + bool down_infeasible = std::isinf(sb_down); + bool up_infeasible = std::isinf(sb_up); + bool down_cutoff = false; + bool up_cutoff = false; + + if (!down_infeasible && std::isfinite(sb_down) && std::isfinite(current_upper)) { + down_cutoff = (sb_down + root_objective_ > current_upper + settings_.dual_tol); + down_infeasible = down_cutoff; + } + if (!up_infeasible && std::isfinite(sb_up) && std::isfinite(current_upper)) { + up_cutoff = (sb_up + root_objective_ > current_upper + settings_.dual_tol); + up_infeasible = up_cutoff; + } + + if (down_infeasible && up_infeasible) { + bool truly_infeasible = std::isinf(sb_down) && std::isinf(sb_up); + if (truly_infeasible) { + settings_.log.printf("Strong branching: both branches infeasible for variable %d\n", j); + return mip_status_t::INFEASIBLE; + } + // Might happen if the incumbent is already the optimal + settings_.log.printf("Strong branching: both branches fathomed for variable %d\n", j); + bool has_incumbent = false; + mutex_upper_.lock(); + has_incumbent = incumbent_.has_incumbent; + mutex_upper_.unlock(); + assert(has_incumbent); + solver_status_ = mip_status_t::OPTIMAL; + set_final_solution(solution, upper_bound_.load()); + return solver_status_; + } + if (down_infeasible) { + mutex_original_lp_.lock(); + f_t new_lb = std::ceil(root_relax_soln_.x[j]); + if (new_lb > original_lp_.lower[j]) { + settings_.log.debug("SB tighten var %d: lb %e -> %e (%s)", + j, + original_lp_.lower[j], + new_lb, + down_cutoff ? "cutoff" : "infeasible"); + original_lp_.lower[j] = new_lb; + num_tightened++; + if (down_cutoff) { + num_cutoff++; + } else { + num_infeasible++; + } + } + mutex_original_lp_.unlock(); + } + if (up_infeasible) { + mutex_original_lp_.lock(); + f_t new_ub = std::floor(root_relax_soln_.x[j]); + if (new_ub < original_lp_.upper[j]) { + settings_.log.debug("SB tighten var %d: ub %e -> %e (%s)", + j, + original_lp_.upper[j], + new_ub, + up_cutoff ? "cutoff" : "infeasible"); + original_lp_.upper[j] = new_ub; + num_tightened++; + if (up_cutoff) { + num_cutoff++; + } else { + num_infeasible++; + } + } + mutex_original_lp_.unlock(); + } + } + if (num_tightened > 0) { + settings_.log.printf( + "Strong branching bounds tightening: %d tightened (%d infeasible, %d cutoff)\n", + num_tightened, + num_infeasible, + num_cutoff); + + std::vector bounds_changed(original_lp_.num_cols, true); + std::vector row_sense; + std::vector new_lower; + std::vector new_upper; + mutex_original_lp_.lock(); + new_lower = original_lp_.lower; + new_upper = original_lp_.upper; + mutex_original_lp_.unlock(); + bounds_strengthening_t sb_presolve(original_lp_, Arow_, row_sense, var_types_); + bool feasible = + sb_presolve.bounds_strengthening(settings_, bounds_changed, new_lower, new_upper); + i_t num_fixed = 0; + if (feasible) { + num_fixed = prune_fixed_fractional_variables(new_lower, new_upper, settings_, fractional); + } + mutex_original_lp_.lock(); + original_lp_.lower = new_lower; + original_lp_.upper = new_upper; + mutex_original_lp_.unlock(); + if (!feasible) { + settings_.log.printf("Strong branching bounds propagation detected infeasibility\n"); + return mip_status_t::INFEASIBLE; + } + if (num_fixed > 0) { + settings_.log.printf( + "Strong branching bounds tightening: %d variables fixed (%d from propagation)\n", + num_fixed, + num_fixed - num_tightened); + num_fractional = fractional.size(); + } + + // If no fractionals remain after the fixings - perform a resolve + // to get fractionals to branch on, or return optimality if the root relaxation is integer + if (num_fractional == 0) { + lp_settings.concurrent_halt = NULL; + i_t iter = 0; + bool initialize_basis = false; + dual::status_t lp_status = dual_phase2_with_advanced_basis(2, + 0, + initialize_basis, + exploration_stats_.start_time, + original_lp_, + lp_settings, + root_vstatus_, + basis_update, + basic_list, + nonbasic_list, + root_relax_soln_, + iter, + edge_norms_); + exploration_stats_.total_lp_iters += iter; + root_objective_ = compute_objective(original_lp_, root_relax_soln_.x); + if (lp_status == dual::status_t::OPTIMAL) { + fractional.clear(); + num_fractional = + fractional_variables(settings_, root_relax_soln_.x, var_types_, fractional); + if (num_fractional == 0) { + set_solution_at_root(solution, cut_info); + return mip_status_t::OPTIMAL; + } + } else if (lp_status == dual::status_t::TIME_LIMIT) { + solver_status_ = mip_status_t::TIME_LIMIT; + set_final_solution(solution, root_objective_); + return solver_status_; + } else { + settings_.log.printf("LP re-solve after SB tightening returned status %d\n", lp_status); + return mip_status_t::NUMERICAL; + } + } + } + } + if (settings_.reduced_cost_strengthening >= 2 && upper_bound_.load() < last_upper_bound) { std::vector lower_bounds; std::vector upper_bounds; @@ -2379,39 +2565,59 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut if (num_fixed > 0) { std::vector bounds_changed(original_lp_.num_cols, true); std::vector row_sense; - + std::vector new_lower = lower_bounds; + std::vector new_upper = upper_bounds; bounds_strengthening_t node_presolve(original_lp_, Arow_, row_sense, var_types_); - - mutex_original_lp_.lock(); - original_lp_.lower = lower_bounds; - original_lp_.upper = upper_bounds; - bool feasible = node_presolve.bounds_strengthening( - settings_, bounds_changed, original_lp_.lower, original_lp_.upper); - mutex_original_lp_.unlock(); + bool feasible = + node_presolve.bounds_strengthening(settings_, bounds_changed, new_lower, new_upper); if (!feasible) { settings_.log.printf("Bound strengthening failed\n"); return mip_status_t::NUMERICAL; // We had a feasible integer solution, but bound // strengthening thinks we are infeasible. } - // Go through and check the fractional variables and remove any that are now fixed to their - // bounds - std::vector to_remove(fractional.size(), 0); - i_t num_to_remove = 0; - for (i_t k = 0; k < fractional.size(); k++) { - const i_t j = fractional[k]; - if (std::abs(original_lp_.upper[j] - original_lp_.lower[j]) < settings_.fixed_tol) { - to_remove[k] = 1; - num_to_remove++; - } - } - if (num_to_remove > 0) { - std::vector new_fractional; - new_fractional.reserve(fractional.size() - num_to_remove); - for (i_t k = 0; k < fractional.size(); k++) { - if (!to_remove[k]) { new_fractional.push_back(fractional[k]); } - } - fractional = new_fractional; + i_t num_fixed = prune_fixed_fractional_variables(new_lower, new_upper, settings_, fractional); + mutex_original_lp_.lock(); + original_lp_.lower = new_lower; + original_lp_.upper = new_upper; + mutex_original_lp_.unlock(); + if (num_fixed > 0) { num_fractional = fractional.size(); + if (num_fractional == 0) { + lp_settings.concurrent_halt = NULL; + i_t iter = 0; + bool initialize_basis = false; + dual::status_t lp_status = dual_phase2_with_advanced_basis(2, + 0, + initialize_basis, + exploration_stats_.start_time, + original_lp_, + lp_settings, + root_vstatus_, + basis_update, + basic_list, + nonbasic_list, + root_relax_soln_, + iter, + edge_norms_); + exploration_stats_.total_lp_iters += iter; + root_objective_ = compute_objective(original_lp_, root_relax_soln_.x); + if (lp_status == dual::status_t::OPTIMAL) { + fractional.clear(); + num_fractional = + fractional_variables(settings_, root_relax_soln_.x, var_types_, fractional); + if (num_fractional == 0) { + set_solution_at_root(solution, cut_info); + return mip_status_t::OPTIMAL; + } + } else if (lp_status == dual::status_t::TIME_LIMIT) { + solver_status_ = mip_status_t::TIME_LIMIT; + set_final_solution(solution, root_objective_); + return solver_status_; + } else { + settings_.log.printf("LP re-solve after RC tightening returned status %d\n", lp_status); + return mip_status_t::NUMERICAL; + } + } } } } diff --git a/cpp/src/branch_and_bound/pseudo_costs.cpp b/cpp/src/branch_and_bound/pseudo_costs.cpp index ee7e2f7803..7633a33990 100644 --- a/cpp/src/branch_and_bound/pseudo_costs.cpp +++ b/cpp/src/branch_and_bound/pseudo_costs.cpp @@ -34,6 +34,7 @@ void strong_branch_helper(i_t start, const std::vector& root_soln, const std::vector& root_vstatus, const std::vector& edge_norms, + f_t upper_bound, pseudo_costs_t& pc) { raft::common::nvtx::range scope("BB::strong_branch_helper"); @@ -62,6 +63,7 @@ void strong_branch_helper(i_t start, if (elapsed_time > settings.time_limit) { break; } child_settings.time_limit = std::max(0.0, settings.time_limit - elapsed_time); child_settings.iteration_limit = 200; + child_settings.cut_off = upper_bound + settings.dual_tol; lp_solution_t solution(original_lp.num_rows, original_lp.num_cols); i_t iter = 0; std::vector vstatus = root_vstatus; @@ -80,7 +82,8 @@ void strong_branch_helper(i_t start, if (status == dual::status_t::DUAL_UNBOUNDED) { // LP was infeasible obj = std::numeric_limits::infinity(); - } else if (status == dual::status_t::OPTIMAL || status == dual::status_t::ITERATION_LIMIT) { + } else if (status == dual::status_t::OPTIMAL || status == dual::status_t::ITERATION_LIMIT || + status == dual::status_t::CUTOFF) { obj = compute_objective(child_problem, solution.x); } else { settings.log.debug("Thread id %2d remaining %d variable %d branch %d status %d\n", @@ -307,6 +310,7 @@ void strong_branching(const user_problem_t& original_problem, f_t root_obj, const std::vector& root_vstatus, const std::vector& edge_norms, + f_t upper_bound, pseudo_costs_t& pc) { pc.resize(original_lp.num_cols); @@ -432,6 +436,7 @@ void strong_branching(const user_problem_t& original_problem, root_soln, root_vstatus, edge_norms, + upper_bound, pc); } } @@ -786,6 +791,7 @@ template void strong_branching(const user_problem_t& o double root_obj, const std::vector& root_vstatus, const std::vector& edge_norms, + double upper_bound, pseudo_costs_t& pc); #endif diff --git a/cpp/src/branch_and_bound/pseudo_costs.hpp b/cpp/src/branch_and_bound/pseudo_costs.hpp index 6b6c6917b6..4b7a852d31 100644 --- a/cpp/src/branch_and_bound/pseudo_costs.hpp +++ b/cpp/src/branch_and_bound/pseudo_costs.hpp @@ -527,6 +527,7 @@ void strong_branching(const user_problem_t& original_problem, f_t root_obj, const std::vector& root_vstatus, const std::vector& edge_norms, + f_t upper_bound, pseudo_costs_t& pc); } // namespace cuopt::linear_programming::dual_simplex diff --git a/cpp/src/mip_heuristics/solver.cu b/cpp/src/mip_heuristics/solver.cu index 235d4500d2..16262ff634 100644 --- a/cpp/src/mip_heuristics/solver.cu +++ b/cpp/src/mip_heuristics/solver.cu @@ -218,8 +218,10 @@ solution_t mip_solver_t::run_solver() branch_and_bound_settings.knapsack_cuts = context.settings.knapsack_cuts; branch_and_bound_settings.strong_chvatal_gomory_cuts = context.settings.strong_chvatal_gomory_cuts; + // default is reduced cost strengthening ON branch_and_bound_settings.reduced_cost_strengthening = - context.settings.reduced_cost_strengthening; + context.settings.reduced_cost_strengthening < 0 ? 2 + : context.settings.reduced_cost_strengthening; branch_and_bound_settings.cut_change_threshold = context.settings.cut_change_threshold; branch_and_bound_settings.cut_min_orthogonality = context.settings.cut_min_orthogonality; branch_and_bound_settings.mip_batch_pdlp_strong_branching =