diff --git a/cpp/src/branch_and_bound/branch_and_bound.cpp b/cpp/src/branch_and_bound/branch_and_bound.cpp index d204596c18..57388543a2 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.cpp +++ b/cpp/src/branch_and_bound/branch_and_bound.cpp @@ -1213,7 +1213,6 @@ struct deterministic_bfs_policy_t branch_direction_t dir) override { switch (status) { - case node_status_t::INFEASIBLE: this->worker.record_infeasible(node); break; case node_status_t::FATHOMED: this->worker.record_fathomed(node, node->lower_bound); break; case node_status_t::INTEGER_FEASIBLE: this->worker.record_integer_solution(node, node->lower_bound); @@ -1365,6 +1364,23 @@ void branch_and_bound_t::snap_to_lattice(const lp_problem_t& } } +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 lies on the lattice strictly below upper_bound. + f_t k = std::ceil((upper_bound - bias) / step - settings_.integer_tol); + return (k - 1) * step + bias; + } else if (lp.objective_is_integral) { + // Any improving feasible solution has an integer objective strictly below upper_bound. + return std::ceil(upper_bound - settings_.integer_tol) - 1; + } + return upper_bound; +} + template template std::pair branch_and_bound_t::update_tree_impl( @@ -1384,15 +1400,14 @@ std::pair branch_and_bound_t::updat worker->recompute_basis = true; if (lp_status == dual_status_t::DUAL_UNBOUNDED) { - node_ptr->lower_bound = inf; - policy.graphviz(search_tree, node_ptr, "infeasible", 0.0); - search_tree.update(node_ptr, node_status_t::INFEASIBLE); - status = node_status_t::INFEASIBLE; + node_ptr->lower_bound = upper_bound; + policy.graphviz(search_tree, node_ptr, "infeasible", upper_bound); + search_tree.update(node_ptr, node_status_t::FATHOMED); + status = node_status_t::FATHOMED; } else if (lp_status == dual_status_t::CUTOFF) { - f_t leaf_obj = compute_objective(leaf_problem, leaf_solution.x); node_ptr->lower_bound = upper_bound; - policy.graphviz(search_tree, node_ptr, "cut off", leaf_obj); + policy.graphviz(search_tree, node_ptr, "cut off", upper_bound); search_tree.update(node_ptr, node_status_t::FATHOMED); status = node_status_t::FATHOMED; @@ -1584,21 +1599,8 @@ 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; - } + const f_t cutoff = compute_node_cutoff(worker->leaf_problem, upper_bound_.load()); + lp_settings.cut_off = cutoff + settings_.dual_tol; 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; } @@ -1642,6 +1644,13 @@ dual_status_t branch_and_bound_t::solve_node_lp( 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; } + if (worker->search_strategy == search_strategy_t::BEST_FIRST) { + if (!worker->domain.propagate_objective( + worker->Arow, worker->var_types, settings_, worker->leaf_problem, cutoff)) { + return dual_status_t::CUTOFF; + } + } + i_t node_iter = 0; f_t lp_start_time = tic(); f_t node_work_estimate = 0.0; @@ -1729,7 +1738,6 @@ void branch_and_bound_t::plunge_with(bfs_worker_t* worker, worker->recompute_basis = true; worker->ensure_orbital_fixing(); - // The branching of the start node itself is applied when its LP is solved. if (start_node->parent != nullptr) { start_node->parent->rebuild_variable_bounds(worker->domain); @@ -1741,7 +1749,7 @@ void branch_and_bound_t::plunge_with(bfs_worker_t* worker, if (!is_start_feasible) { start_node->lower_bound = inf; search_tree_.graphviz_node(settings_.log, start_node, "infeasible", inf); - search_tree_.update(start_node, node_status_t::INFEASIBLE); + search_tree_.update(start_node, node_status_t::FATHOMED); ++exploration_stats_.nodes_explored; --exploration_stats_.nodes_unexplored; return; @@ -4260,7 +4268,6 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut return mip_status_t::OPTIMAL; } - is_running_ = true; lp_settings.concurrent_halt = settings_.concurrent_halt ? settings_.concurrent_halt : &node_concurrent_halt_; reduced_cost_bounds_t reduced_cost_bounds(original_lp_.num_cols); @@ -4355,6 +4362,7 @@ mip_status_t branch_and_bound_t::solve(mip_solution_t& solut } if (num_fractional != 0 && settings_.max_cut_passes > 0) { print_table_header(); } + is_running_ = true; cut_pool_t cut_pool(original_lp_.num_cols, settings_); cut_generation_t cut_generation(cut_pool, @@ -5269,7 +5277,8 @@ node_status_t branch_and_bound_t::solve_node_deterministic( simplex_solver_settings_t lp_settings = settings_; lp_settings.set_log(false); - lp_settings.cut_off = worker.local_upper_bound + settings_.dual_tol; + const f_t cutoff = compute_node_cutoff(worker.leaf_problem, worker.local_upper_bound); + lp_settings.cut_off = cutoff + settings_.dual_tol; lp_settings.inside_mip = 2; lp_settings.time_limit = remaining_time; lp_settings.scale_columns = false; @@ -5284,19 +5293,32 @@ node_status_t branch_and_bound_t::solve_node_deterministic( feasible = worker.rebuild_bounds_from_stack(original_lp_, settings_); } + const bool is_cutoff = + feasible && !worker.domain.propagate_objective( + worker.Arow, worker.var_types, settings_, worker.leaf_problem, cutoff); + if (settings_.deterministic) { // TEMP APPROXIMATION; worker.work_context.record_work_sync_on_horizon(worker.domain.last_nnz_processed / 1e8); } if (!feasible) { - node_ptr->lower_bound = std::numeric_limits::infinity(); - search_tree.update(node_ptr, node_status_t::INFEASIBLE); - worker.record_infeasible(node_ptr); + node_ptr->lower_bound = worker.local_upper_bound; + search_tree.update(node_ptr, node_status_t::FATHOMED); + worker.record_fathomed(node_ptr, worker.local_upper_bound); --exploration_stats_.nodes_unexplored; ++exploration_stats_.nodes_explored; worker.recompute_basis = true; - return node_status_t::INFEASIBLE; + return node_status_t::FATHOMED; + } + + if (is_cutoff) { + ++exploration_stats_.nodes_explored; + --exploration_stats_.nodes_unexplored; + deterministic_bfs_policy_t policy{*this, worker}; + auto [status, round_dir] = + update_tree_impl(node_ptr, search_tree, &worker, dual_status_t::CUTOFF, policy); + return status; } // Solve LP relaxation @@ -5894,12 +5916,15 @@ void branch_and_bound_t::deterministic_dive( // Setup LP settings simplex_solver_settings_t lp_settings = settings_; lp_settings.set_log(false); - lp_settings.cut_off = worker.local_upper_bound + settings_.dual_tol; + const f_t cutoff = compute_node_cutoff(worker.leaf_problem, worker.local_upper_bound); + lp_settings.cut_off = cutoff + settings_.dual_tol; lp_settings.inside_mip = 2; lp_settings.time_limit = remaining_time; lp_settings.scale_columns = false; - feasible = worker.update_variable_bounds(node_ptr, settings_); + feasible = worker.update_variable_bounds(node_ptr, settings_) && + worker.domain.propagate_objective( + worker.Arow, worker.var_types, settings_, worker.leaf_problem, cutoff); if (settings_.deterministic) { // TEMP APPROXIMATION; diff --git a/cpp/src/branch_and_bound/branch_and_bound.hpp b/cpp/src/branch_and_bound/branch_and_bound.hpp index e17e3a5043..83fa59471c 100644 --- a/cpp/src/branch_and_bound/branch_and_bound.hpp +++ b/cpp/src/branch_and_bound/branch_and_bound.hpp @@ -356,6 +356,10 @@ class branch_and_bound_t { mip_node_t* node_ptr, f_t leaf_obj); + // Largest objective value of a solution that improves on upper_bound. This is the value just + // below upper_bound on the objective lattice if the objective is integral or moves in steps. + f_t compute_node_cutoff(const simplex::lp_problem_t& lp, f_t upper_bound) const; + // Launch a new best-first worker from a given bfs worker. void launch_bfs_worker(bfs_worker_t* worker); diff --git a/cpp/src/branch_and_bound/deterministic_workers.hpp b/cpp/src/branch_and_bound/deterministic_workers.hpp index 6174713dd5..a0dde2b6fd 100644 --- a/cpp/src/branch_and_bound/deterministic_workers.hpp +++ b/cpp/src/branch_and_bound/deterministic_workers.hpp @@ -261,15 +261,6 @@ class deterministic_bfs_worker_t ++total_nodes_pruned; } - void record_infeasible(mip_node_t* node) - { - record_event( - bb_event_t::make_infeasible(this->clock, this->worker_id, node->creation_seq)); - ++nodes_processed_this_horizon; - ++this->total_nodes_processed; - ++total_nodes_infeasible; - } - void record_numerical(mip_node_t* node) { record_event( diff --git a/cpp/src/branch_and_bound/mip_node.hpp b/cpp/src/branch_and_bound/mip_node.hpp index bf06c83d9e..1ecc3e1287 100644 --- a/cpp/src/branch_and_bound/mip_node.hpp +++ b/cpp/src/branch_and_bound/mip_node.hpp @@ -27,16 +27,15 @@ namespace cuopt::mathematical_optimization::mip { enum class node_status_t : int { PENDING = 0, // Node is still in the tree, waiting to be solved INTEGER_FEASIBLE = 1, // Node has an integer feasible solution - INFEASIBLE = 2, // Node is infeasible - FATHOMED = 3, // Node objective is greater than the upper bound - HAS_CHILDREN = 4, // Node has children to explore - NUMERICAL = 5 // Encountered numerical issue when solving the LP relaxation + FATHOMED = 2, // Node is infeasible or has objective is greater than the upper bound + HAS_CHILDREN = 3, // Node has children to explore + NUMERICAL = 4 // Encountered numerical issue when solving the LP relaxation }; inline bool inactive_status(node_status_t status) { return (status == node_status_t::FATHOMED || status == node_status_t::INTEGER_FEASIBLE || - status == node_status_t::INFEASIBLE || status == node_status_t::NUMERICAL); + status == node_status_t::NUMERICAL); } template diff --git a/cpp/src/dual_simplex/domain.cpp b/cpp/src/dual_simplex/domain.cpp index 07e0ecb106..0f083dfe39 100644 --- a/cpp/src/dual_simplex/domain.cpp +++ b/cpp/src/dual_simplex/domain.cpp @@ -85,6 +85,12 @@ void domain_t::compute_activities(const csr_matrix_t& Arow, for (i_t i = 0; i < Arow.m; ++i) { compute_row_activity(i, Arow, lp.lower, lp.upper); } + + objective_vars.clear(); + objective_vars.reserve(lp.num_cols); + for (i_t j = 0; j < lp.num_cols; ++j) { + if (lp.objective[j] != 0) { objective_vars.push_back(j); } + } } template @@ -422,6 +428,95 @@ bool domain_t::apply_and_propagate(const csr_matrix_t& Arow, return run_bound_propagation(Arow, var_types, settings, lp); } +template +bool domain_t::propagate_objective(const csr_matrix_t& Arow, + const std::vector& var_types, + const simplex_solver_settings_t& settings, + lp_problem_t& lp, + f_t cutoff) +{ + if (!std::isfinite(cutoff) || objective_vars.empty()) { return true; } + + // The work is added to the one of the preceding propagation, so the callers record both. + const size_t prior_nnz = last_nnz_processed; + size_t work = 0; + const f_t rhs = cutoff + settings.primal_tol * std::max(1.0, std::abs(cutoff)); + bool feasible = true; + + for (i_t round = 0; round < params.max_objective_rounds; ++round) { + // Minimum objective over the bounds in lp, with the infinite contributions flagged as in + // compute_row_activity. + f_t min = 0; + f_t min_err = 0; + i_t min_inf = 0; + for (i_t j : objective_vars) { + const f_t c_j = lp.objective[j]; + const f_t bound_min = c_j < 0 ? lp.upper[j] : lp.lower[j]; + const f_t alpha_min = c_j * bound_min; + if (std::isfinite(alpha_min) && std::abs(alpha_min) < params.huge_value) { + dot2_add(c_j, bound_min, min, min_err); + } else { + ++min_inf; + } + } + work += objective_vars.size(); + + if (min_inf == 0 && min + min_err > rhs) { + feasible = false; + break; + } + if (min_inf > 1) { break; } + + const size_t num_changes = bound_changes.size(); + for (i_t j : objective_vars) { + const f_t lb = lp.lower[j]; + const f_t ub = lp.upper[j]; + if (lb == ub) { continue; } + + const f_t c_j = lp.objective[j]; + const f_t bound_min = c_j < 0 ? ub : lb; + const f_t alpha_min = c_j * bound_min; + const bool alpha_min_inf = + !std::isfinite(alpha_min) || std::abs(alpha_min) >= params.huge_value; + + // With a single infinite contribution, only its variable can be bounded. + if (min_inf == 1 && !alpha_min_inf) { continue; } + + // c_j x_j <= rhs - (minimum objective of the other variables): an upper bound if c_j > 0, a + // lower bound otherwise. The slack is formed with Dot2, as in the row propagation. + f_t slack = rhs; + f_t slack_err = -min_err; + dot2_add(-1.0, min, slack, slack_err); + if (min_inf == 0) { dot2_add(c_j, bound_min, slack, slack_err); } + const f_t gamma = (slack + slack_err) / c_j; + if (std::abs(gamma) >= params.max_derived_bound) { continue; } + + const f_t new_lower = c_j > 0 ? lb : gamma; + const f_t new_upper = c_j > 0 ? gamma : ub; + if (!tighten_bounds( + j, new_lower, new_upper, bound_change_origin_t::OBJECTIVE, var_types, settings, lp)) { + feasible = false; + break; + } + } + work += objective_vars.size(); + + if (!feasible || bound_changes.size() == num_changes) { break; } + + // Propagate the rows of the tightened variables, which may raise the minimum objective. + feasible = run_bound_propagation(Arow, var_types, settings, lp); + work += last_nnz_processed; + if (!feasible) { break; } + } + + if (!feasible) { clear_queues(); } + + work += nnz_processed; + nnz_processed = 0; + last_nnz_processed = prior_nnz + work; + return feasible; +} + template void domain_t::backtrack_to_parent(lp_problem_t& lp) { diff --git a/cpp/src/dual_simplex/domain.hpp b/cpp/src/dual_simplex/domain.hpp index ac823e1049..0c72530650 100644 --- a/cpp/src/dual_simplex/domain.hpp +++ b/cpp/src/dual_simplex/domain.hpp @@ -59,6 +59,8 @@ struct domain_params { double min_improvement_factor = 1e3; // Derived bounds larger than this in magnitude are discarded double max_derived_bound = 1e8; + // Maximum number of rounds of objective propagation followed by row propagation + int max_objective_rounds = 3; }; // The finite parts of the activities are accumulated with compensated (Dot2) summation: the @@ -164,7 +166,18 @@ class domain_t { bound_change_t bound_change, lp_problem_t& lp); + // Propagates c^T x <= cutoff, followed by the rows affected by its tightenings, for up to + // max_objective_rounds rounds. Returns false if no solution within the bounds in lp can reach + // the cutoff. + bool propagate_objective(const csr_matrix_t& Arow, + const std::vector& var_types, + const simplex_solver_settings_t& settings, + lp_problem_t& lp, + f_t cutoff); + size_t last_nnz_processed{0}; + size_t objective_fixings{0}; + size_t objective_cutoffs{0}; // Implications applied when an integer variable is fixed to 0 or 1. Either may be null, and the // clique table is skipped until it is ready. @@ -181,6 +194,9 @@ class domain_t { std::vector var_queued; circular_deque_t var_queue; + // Variables with a nonzero objective coefficient. + std::vector objective_vars; + std::vector> bound_changes; size_t nnz_processed{0};