From 4255e18c56220c4eb75911e310b79edd3da8462a Mon Sep 17 00:00:00 2001 From: Laurent Perron Date: Thu, 28 Sep 2023 14:28:27 +0200 Subject: [PATCH] [CP-SAT] improve lns for makespan objective; fix jeasibility jump check failed; fix non deterministic search with gap limits --- ortools/sat/cp_model_lns.cc | 233 ++++++++++++++++--------- ortools/sat/cp_model_lns.h | 29 +-- ortools/sat/cp_model_solver.cc | 16 +- ortools/sat/disjunctive.cc | 45 ++--- ortools/sat/disjunctive.h | 16 +- ortools/sat/docs/README.md | 13 +- ortools/sat/docs/boolean_logic.md | 23 --- ortools/sat/docs/channeling.md | 15 -- ortools/sat/docs/integer_arithmetic.md | 33 ---- ortools/sat/docs/model.md | 14 -- ortools/sat/docs/scheduling.md | 67 +++---- ortools/sat/docs/solver.md | 25 --- ortools/sat/feasibility_jump.cc | 162 ++++++++++------- ortools/sat/feasibility_jump.h | 47 ++++- ortools/sat/intervals.cc | 4 +- ortools/sat/sat_parameters.proto | 12 +- ortools/sat/synchronization.cc | 4 +- 17 files changed, 390 insertions(+), 368 deletions(-) diff --git a/ortools/sat/cp_model_lns.cc b/ortools/sat/cp_model_lns.cc index 5385f74968..6ba1b86011 100644 --- a/ortools/sat/cp_model_lns.cc +++ b/ortools/sat/cp_model_lns.cc @@ -35,8 +35,6 @@ #include "absl/strings/str_cat.h" #include "absl/strings/str_join.h" #include "absl/synchronization/mutex.h" -#include "absl/time/clock.h" -#include "absl/time/time.h" #include "absl/types/span.h" #include "ortools/base/logging.h" #include "ortools/base/stl_util.h" @@ -395,52 +393,48 @@ Neighborhood NeighborhoodGeneratorHelper::NoNeighborhood() const { return neighborhood; } +bool NeighborhoodGeneratorHelper::IntervalIsActive( + int index, const CpSolverResponse& initial_solution) const { + const ConstraintProto& interval_ct = ModelProto().constraints(index); + // We only look at intervals that are performed in the solution. The + // unperformed intervals should be automatically freed during the generation + // phase. + if (interval_ct.enforcement_literal().size() == 1) { + const int enforcement_ref = interval_ct.enforcement_literal(0); + const int enforcement_var = PositiveRef(enforcement_ref); + const int value = initial_solution.solution(enforcement_var); + if (RefIsPositive(enforcement_ref) == (value == 0)) return false; + } + + for (const int v : interval_ct.interval().start().vars()) { + if (!IsConstant(v)) return true; + } + for (const int v : interval_ct.interval().size().vars()) { + if (!IsConstant(v)) return true; + } + for (const int v : interval_ct.interval().end().vars()) { + if (!IsConstant(v)) return true; + } + return false; +} + +void NeighborhoodGeneratorHelper::FilterInactiveIntervals( + absl::Span unfiltered_intervals, + const CpSolverResponse& initial_solution, + std::vector* filtered_intervals) const { + filtered_intervals->clear(); + absl::ReaderMutexLock lock(&domain_mutex_); + for (const int i : unfiltered_intervals) { + if (IntervalIsActive(i, initial_solution)) filtered_intervals->push_back(i); + } +} + std::vector NeighborhoodGeneratorHelper::GetActiveIntervals( const CpSolverResponse& initial_solution) const { - std::vector active_intervals; - absl::ReaderMutexLock lock(&domain_mutex_); - for (const int i : TypeToConstraints(ConstraintProto::kInterval)) { - const ConstraintProto& interval_ct = ModelProto().constraints(i); - // We only look at intervals that are performed in the solution. The - // unperformed intervals should be automatically freed during the generation - // phase. - if (interval_ct.enforcement_literal().size() == 1) { - const int enforcement_ref = interval_ct.enforcement_literal(0); - const int enforcement_var = PositiveRef(enforcement_ref); - const int value = initial_solution.solution(enforcement_var); - if (RefIsPositive(enforcement_ref) == (value == 0)) { - continue; - } - } - - // We filter out fixed intervals. Because of presolve, if there is an - // enforcement literal, it cannot be fixed. - if (interval_ct.enforcement_literal().empty()) { - bool is_constant = true; - for (const int v : interval_ct.interval().start().vars()) { - if (!IsConstant(v)) { - is_constant = false; - break; - } - } - for (const int v : interval_ct.interval().size().vars()) { - if (!IsConstant(v)) { - is_constant = false; - break; - } - } - for (const int v : interval_ct.interval().end().vars()) { - if (!IsConstant(v)) { - is_constant = false; - break; - } - } - if (is_constant) continue; - } - - active_intervals.push_back(i); - } - return active_intervals; + std::vector result; + FilterInactiveIntervals(TypeToConstraints(ConstraintProto::kInterval), + initial_solution, &result); + return result; } std::vector> @@ -513,15 +507,22 @@ struct StartEndIndex { int64_t start; int64_t end; int index_in_input_vector; + double noise; bool operator<(const StartEndIndex& o) const { - return std::tie(start, end, index_in_input_vector) < - std::tie(o.start, o.end, o.index_in_input_vector); + return std::tie(start, end, noise, index_in_input_vector) < + std::tie(o.start, o.end, o.noise, o.index_in_input_vector); } }; +struct TimePartition { + std::vector indices_before_selected; + std::vector selected_indices; + std::vector indices_after_selected; +}; + // Selects all intervals in a random time window to meet the difficulty // requirement. -std::vector SelectIndicesInRandomTimeWindow( +TimePartition PartitionIndicesAroundRandomTimeWindow( const std::vector& intervals, const CpModelProto& model_proto, const CpSolverResponse& initial_solution, double difficulty, absl::BitGenRef random) { @@ -529,22 +530,12 @@ std::vector SelectIndicesInRandomTimeWindow( for (int index = 0; index < intervals.size(); ++index) { const int interval = intervals[index]; const ConstraintProto& interval_ct = model_proto.constraints(interval); - // We only look at intervals that are performed in the solution. The - // unperformed intervals should be automatically freed during the - // generation phase. - if (interval_ct.enforcement_literal().size() == 1) { - const int enforcement_ref = interval_ct.enforcement_literal(0); - const int enforcement_var = PositiveRef(enforcement_ref); - const int64_t value = initial_solution.solution(enforcement_var); - if (RefIsPositive(enforcement_ref) == (value == 0)) { - continue; - } - } const int64_t start_value = GetLinearExpressionValue( interval_ct.interval().start(), initial_solution); const int64_t end_value = GetLinearExpressionValue( interval_ct.interval().end(), initial_solution); - start_end_indices.push_back({start_value, end_value, index}); + start_end_indices.push_back( + {start_value, end_value, index, absl::Uniform(random, 0., 1.0)}); } if (start_end_indices.empty()) return {}; @@ -566,12 +557,22 @@ std::vector SelectIndicesInRandomTimeWindow( std::sort(start_end_indices.begin() + random_start_index, start_end_indices.end(), [](const StartEndIndex& a, const StartEndIndex& b) { - return std::tie(a.end, a.index_in_input_vector) < - std::tie(b.end, b.index_in_input_vector); + return std::tie(a.end, a.noise, a.index_in_input_vector) < + std::tie(b.end, b.noise, b.index_in_input_vector); }); - std::vector result; - for (int i = random_start_index; i < random_start_index + relaxed_size; ++i) { - result.push_back(start_end_indices[i].index_in_input_vector); + TimePartition result; + int i = 0; + for (; i < random_start_index; ++i) { + result.indices_before_selected.push_back( + start_end_indices[i].index_in_input_vector); + } + for (; i < random_start_index + relaxed_size; ++i) { + result.selected_indices.push_back( + start_end_indices[i].index_in_input_vector); + } + for (; i < start_end_indices.size(); ++i) { + result.indices_after_selected.push_back( + start_end_indices[i].index_in_input_vector); } return result; } @@ -999,6 +1000,13 @@ Neighborhood NeighborhoodGeneratorHelper::FixGivenVariables( const absl::flat_hash_set& variables_to_fix) const { Neighborhood neighborhood; + // TODO(user): Maybe relax all variables in the objective when the number + // is small or negligible compared to the number of variables. + const int unique_objective_variable = + model_proto_.has_objective() && model_proto_.objective().vars_size() == 1 + ? model_proto_.objective().vars(0) + : -1; + // Fill in neighborhood.delta all variable domains. { absl::ReaderMutexLock domain_lock(&domain_mutex_); @@ -1017,7 +1025,7 @@ Neighborhood NeighborhoodGeneratorHelper::FixGivenVariables( const Domain domain = ReadDomainFromProto(current_var); const int64_t base_value = base_solution.solution(i); - if (variables_to_fix.contains(i)) { + if (variables_to_fix.contains(i) && i != unique_objective_variable) { if (domain.Contains(base_value)) { new_var->add_domain(base_value); new_var->add_domain(base_value); @@ -1171,7 +1179,7 @@ double NeighborhoodGenerator::GetUCBScore(int64_t total_num_calls) const { return current_average_ + sqrt((2 * log(total_num_calls)) / num_calls_); } -void NeighborhoodGenerator::Synchronize() { +double NeighborhoodGenerator::Synchronize() { absl::MutexLock mutex_lock(&generator_mutex_); // To make the whole update process deterministic, we currently sort the @@ -1182,6 +1190,7 @@ void NeighborhoodGenerator::Synchronize() { int num_fully_solved_in_batch = 0; int num_not_fully_solved_in_batch = 0; + double total_dtime = 0.0; for (const SolveData& data : solve_data_) { ++num_calls_; @@ -1227,7 +1236,7 @@ void NeighborhoodGenerator::Synchronize() { current_average_ = 0.9 * current_average_ + 0.1 * gain_per_time_unit; } - deterministic_time_ += data.deterministic_time; + total_dtime += data.deterministic_time; } // Update the difficulty. @@ -1250,6 +1259,7 @@ void NeighborhoodGenerator::Synchronize() { } solve_data_.clear(); + return total_dtime; } namespace { @@ -1912,7 +1922,8 @@ Neighborhood GenerateSchedulingNeighborhoodFromIntervalPrecedences( } Neighborhood GenerateSchedulingNeighborhoodFromRelaxedIntervals( - const absl::Span intervals_to_relax, + absl::Span intervals_to_relax, + absl::Span variables_to_fix, const CpSolverResponse& initial_solution, absl::BitGenRef random, const NeighborhoodGeneratorHelper& helper) { Neighborhood neighborhood = helper.FullNeighborhood(); @@ -1977,6 +1988,14 @@ Neighborhood GenerateSchedulingNeighborhoodFromRelaxedIntervals( AddPrecedence(before_end, after_start, &neighborhood.delta); } + // fix the extra variables passed as parameters. + for (const int var : variables_to_fix) { + const int value = initial_solution.solution(var); + neighborhood.delta.mutable_variables(var)->clear_domain(); + neighborhood.delta.mutable_variables(var)->add_domain(value); + neighborhood.delta.mutable_variables(var)->add_domain(value); + } + // Set the current solution as a hint. helper.AddSolutionHinting(initial_solution, &neighborhood.delta); neighborhood.is_generated = true; @@ -1992,7 +2011,7 @@ Neighborhood RandomIntervalSchedulingNeighborhoodGenerator::Generate( GetRandomSubset(difficulty, &intervals_to_relax, random); return GenerateSchedulingNeighborhoodFromRelaxedIntervals( - intervals_to_relax, initial_solution, random, helper_); + intervals_to_relax, {}, initial_solution, random, helper_); } Neighborhood RandomPrecedenceSchedulingNeighborhoodGenerator::Generate( @@ -2013,34 +2032,72 @@ Neighborhood SchedulingTimeWindowNeighborhoodGenerator::Generate( if (active_intervals.empty()) return helper_.FullNeighborhood(); - const std::vector selected_indices = - SelectIndicesInRandomTimeWindow(active_intervals, helper_.ModelProto(), - initial_solution, difficulty, random); + const TimePartition partition = PartitionIndicesAroundRandomTimeWindow( + active_intervals, helper_.ModelProto(), initial_solution, difficulty, + random); std::vector intervals_to_relax; - intervals_to_relax.reserve(selected_indices.size()); - for (const int index : selected_indices) { + intervals_to_relax.reserve(partition.selected_indices.size()); + std::vector variables_to_fix; + for (const int index : partition.selected_indices) { intervals_to_relax.push_back(active_intervals[index]); } + + if (helper_.Parameters().push_all_tasks_toward_start()) { + for (const int index : partition.indices_before_selected) { + const int interval = active_intervals[index]; + // Remove the interval so that we do not use it for precedences. + intervals_to_relax.push_back(interval); + + // Fix all variables. + const std::vector vars = + UsedVariables(helper_.ModelProto().constraints(interval)); + variables_to_fix.insert(variables_to_fix.end(), vars.begin(), vars.end()); + } + } + + gtl::STLSortAndRemoveDuplicates(&intervals_to_relax); + gtl::STLSortAndRemoveDuplicates(&variables_to_fix); return GenerateSchedulingNeighborhoodFromRelaxedIntervals( - intervals_to_relax, initial_solution, random, helper_); + intervals_to_relax, variables_to_fix, initial_solution, random, helper_); } Neighborhood SchedulingResourceWindowsNeighborhoodGenerator::Generate( const CpSolverResponse& initial_solution, double difficulty, absl::BitGenRef random) { std::vector intervals_to_relax; + std::vector variables_to_fix; + std::vector active_intervals; for (const std::vector& intervals : intervals_in_constraints_) { - const std::vector selected_indices = SelectIndicesInRandomTimeWindow( - intervals, helper_.ModelProto(), initial_solution, difficulty, random); - for (const int index : selected_indices) { + helper_.FilterInactiveIntervals(intervals, initial_solution, + &active_intervals); + const TimePartition partition = PartitionIndicesAroundRandomTimeWindow( + active_intervals, helper_.ModelProto(), initial_solution, difficulty, + random); + for (const int index : partition.selected_indices) { intervals_to_relax.push_back(intervals[index]); } + + if (helper_.Parameters().push_all_tasks_toward_start()) { + for (const int index : partition.indices_before_selected) { + const int interval = intervals[index]; + // Remove the interval so that we do not use it for precedences. + intervals_to_relax.push_back(interval); + const std::vector vars = + UsedVariables(helper_.ModelProto().constraints(interval)); + variables_to_fix.insert(variables_to_fix.end(), vars.begin(), + vars.end()); + } + } + } + + if (intervals_to_relax.empty() && variables_to_fix.empty()) { + return helper_.FullNeighborhood(); } - if (intervals_to_relax.empty()) return helper_.FullNeighborhood(); gtl::STLSortAndRemoveDuplicates(&intervals_to_relax); + gtl::STLSortAndRemoveDuplicates(&variables_to_fix); return GenerateSchedulingNeighborhoodFromRelaxedIntervals( - intervals_to_relax, initial_solution, random, helper_); + intervals_to_relax, variables_to_fix, initial_solution, random, helper_); } Neighborhood RandomRectanglesPackingNeighborhoodGenerator::Generate( @@ -2073,7 +2130,7 @@ Neighborhood RandomPrecedencesPackingNeighborhoodGenerator::Generate( gtl::STLSortAndRemoveDuplicates(&intervals_to_relax); return GenerateSchedulingNeighborhoodFromRelaxedIntervals( - intervals_to_relax, initial_solution, random, helper_); + intervals_to_relax, {}, initial_solution, random, helper_); } Neighborhood SlicePackingNeighborhoodGenerator::Generate( @@ -2088,11 +2145,13 @@ Neighborhood SlicePackingNeighborhoodGenerator::Generate( projected_intervals.push_back(use_first_dimension ? x : y); } - const std::vector indices_to_relax = - SelectIndicesInRandomTimeWindow(projected_intervals, helper_.ModelProto(), - initial_solution, difficulty, random); + const TimePartition partition = PartitionIndicesAroundRandomTimeWindow( + projected_intervals, helper_.ModelProto(), initial_solution, difficulty, + random); std::vector indices_to_fix(active_rectangles.size(), true); - for (const int index : indices_to_relax) indices_to_fix[index] = false; + for (const int index : partition.selected_indices) { + indices_to_fix[index] = false; + } absl::flat_hash_set variables_to_freeze; for (int index = 0; index < active_rectangles.size(); ++index) { diff --git a/ortools/sat/cp_model_lns.h b/ortools/sat/cp_model_lns.h index 9b2f82f318..8fbd543865 100644 --- a/ortools/sat/cp_model_lns.h +++ b/ortools/sat/cp_model_lns.h @@ -28,13 +28,14 @@ #include "absl/random/bit_gen_ref.h" #include "absl/strings/string_view.h" #include "absl/synchronization/mutex.h" -#include "absl/time/time.h" #include "absl/types/span.h" #include "ortools/sat/cp_model.pb.h" #include "ortools/sat/integer.h" +#include "ortools/sat/model.h" #include "ortools/sat/sat_parameters.pb.h" #include "ortools/sat/subsolver.h" #include "ortools/sat/synchronization.h" +#include "ortools/sat/util.h" #include "ortools/util/adaptative_parameter_value.h" namespace operations_research { @@ -203,6 +204,19 @@ class NeighborhoodGeneratorHelper : public SubSolver { return absl::MakeSpan(type_to_constraints_[type]); } + // Checks if an interval is active w.r.t. the initial_solution. + // An interval is inactive if and only if it is either unperformed in the + // solution or constant in the model. + bool IntervalIsActive(int index, + const CpSolverResponse& initial_solution) const + ABSL_SHARED_LOCKS_REQUIRED(domain_mutex_); + + // Filters a vector of intervals against the initial_solution, and returns + // only the active intervals. + void FilterInactiveIntervals(absl::Span unfiltered_intervals, + const CpSolverResponse& initial_solution, + std::vector* filtered_intervals) const; + // Returns the list of indices of active interval constraints according // to the initial_solution and the parameter lns_focus_on_performed_intervals. // If true, this method returns the list of performed intervals in the @@ -411,8 +425,9 @@ class NeighborhoodGenerator { } // Process all the recently added solve data and update this generator - // score and difficulty. - void Synchronize(); + // score and difficulty. This returns the sum of the deterministic time of + // each SolveData. + double Synchronize(); // Returns a short description of the generator. std::string name() const { return name_; } @@ -455,12 +470,6 @@ class NeighborhoodGenerator { return deterministic_limit_; } - // The sum of the deterministic time spent in this generator. - double deterministic_time() const { - absl::MutexLock mutex_lock(&generator_mutex_); - return deterministic_time_; - } - protected: const std::string name_; const NeighborhoodGeneratorHelper& helper_; @@ -481,7 +490,6 @@ class NeighborhoodGenerator { int64_t num_improving_calls_ = 0; int64_t num_consecutive_non_improving_calls_ = 0; int64_t next_time_limit_bump_ = 50; - double deterministic_time_ = 0.0; double current_average_ = 0.0; }; @@ -605,6 +613,7 @@ class LocalBranchingLpBasedNeighborhoodGenerator // intervals. Neighborhood GenerateSchedulingNeighborhoodFromRelaxedIntervals( absl::Span intervals_to_relax, + absl::Span variables_to_fix, const CpSolverResponse& initial_solution, absl::BitGenRef random, const NeighborhoodGeneratorHelper& helper); diff --git a/ortools/sat/cp_model_solver.cc b/ortools/sat/cp_model_solver.cc index c9639e9584..ab36cc2a80 100644 --- a/ortools/sat/cp_model_solver.cc +++ b/ortools/sat/cp_model_solver.cc @@ -3333,10 +3333,9 @@ class LnsSolver : public SubSolver { } void Synchronize() override { - generator_->Synchronize(); - const double diff = generator_->deterministic_time() - deterministic_time(); - AddTaskDeterministicDuration(diff); - shared_->time_limit->AdvanceDeterministicTime(diff); + const double dtime = generator_->Synchronize(); + AddTaskDeterministicDuration(dtime); + shared_->time_limit->AdvanceDeterministicTime(dtime); } std::vector TableLineStats() const override { @@ -4153,8 +4152,13 @@ CpSolverResponse SolveCpModel(const CpModelProto& model_proto, Model* model) { } if (is_pure_sat) { for (const ConstraintProto& ct : model_proto.constraints()) { - if (ct.constraint_case() != ConstraintProto::ConstraintCase::kBoolOr && - ct.constraint_case() != ConstraintProto::ConstraintCase::kBoolAnd) { + if (ct.constraint_case() != ConstraintProto::kBoolOr && + ct.constraint_case() != ConstraintProto::kBoolAnd) { + is_pure_sat = false; + break; + } + if (ct.constraint_case() == ConstraintProto::kBoolAnd && + ct.enforcement_literal().size() > 1) { is_pure_sat = false; break; } diff --git a/ortools/sat/disjunctive.cc b/ortools/sat/disjunctive.cc index 5bc6e331ca..6c3b700de9 100644 --- a/ortools/sat/disjunctive.cc +++ b/ortools/sat/disjunctive.cc @@ -17,6 +17,7 @@ #include #include +#include "absl/algorithm/container.h" #include "absl/log/check.h" #include "ortools/sat/all_different.h" #include "ortools/sat/integer.h" @@ -508,8 +509,9 @@ bool CombinedDisjunctive::Propagate() { } bool DisjunctiveOverloadChecker::Propagate() { - if (!helper_->SynchronizeAndSetTimeDirection(/*is_forward=*/true)) + if (!helper_->SynchronizeAndSetTimeDirection(/*is_forward=*/true)) { return false; + } // Split problem into independent part. // @@ -522,20 +524,21 @@ bool DisjunctiveOverloadChecker::Propagate() { // all task in the window, [start_min, end_min] is inside the window, and the // end min of any set of task to the left is <= window_start, and the // start_min of any task to the right is >= end_min. - window_.clear(); IntegerValue window_end = kMinIntegerValue; IntegerValue relevant_end; + int window_size = 0; int relevant_size = 0; for (const TaskTime task_time : helper_->TaskByIncreasingShiftedStartMin()) { const int task = task_time.task_index; if (helper_->IsAbsent(task)) continue; + // Extend window. const IntegerValue start_min = task_time.time; if (start_min < window_end) { - window_.push_back(task_time); + window_[window_size++] = task_time; window_end += helper_->SizeMin(task); if (window_end > helper_->EndMax(task)) { - relevant_size = window_.size(); + relevant_size = window_size; relevant_end = window_end; } continue; @@ -544,21 +547,19 @@ bool DisjunctiveOverloadChecker::Propagate() { // Process current window. // We don't need to process the end of the window (after relevant_size) // because these interval can be greedily assembled in a feasible solution. - window_.resize(relevant_size); - if (relevant_size > 0 && !PropagateSubwindow(relevant_end)) { + if (relevant_size > 0 && !PropagateSubwindow(relevant_size, relevant_end)) { return false; } // Start of the next window. - window_.clear(); - window_.push_back(task_time); + window_size = 0; + window_[window_size++] = task_time; window_end = start_min + helper_->SizeMin(task); relevant_size = 0; } // Process last window. - window_.resize(relevant_size); - if (relevant_size > 0 && !PropagateSubwindow(relevant_end)) { + if (relevant_size > 0 && !PropagateSubwindow(relevant_size, relevant_end)) { return false; } @@ -573,24 +574,28 @@ bool DisjunctiveOverloadChecker::Propagate() { // overload after every insertion, but we could use an upper bound of the // theta envelope to save us from checking the actual value. bool DisjunctiveOverloadChecker::PropagateSubwindow( - IntegerValue global_window_end) { + int relevant_size, IntegerValue global_window_end) { // Set up theta tree and task_by_increasing_end_max_. - const int window_size = window_.size(); - theta_tree_.Reset(window_size); + int num_events = 0; task_by_increasing_end_max_.clear(); - for (int i = 0; i < window_size; ++i) { + for (int i = 0; i < relevant_size; ++i) { // No point adding a task if its end_max is too large. - const int task = window_[i].task_index; + const TaskTime& task_time = window_[i]; + const int task = task_time.task_index; const IntegerValue end_max = helper_->EndMax(task); if (end_max < global_window_end) { - task_to_event_[task] = i; + window_[num_events] = task_time; + task_to_event_[task] = num_events; task_by_increasing_end_max_.push_back({task, end_max}); + ++num_events; } } + theta_tree_.Reset(num_events); // Introduce events by increasing end_max, check for overloads. - std::sort(task_by_increasing_end_max_.begin(), - task_by_increasing_end_max_.end()); + // If end_max is the same, we want to add high start-min first. + absl::c_reverse(task_by_increasing_end_max_); + absl::c_stable_sort(task_by_increasing_end_max_); for (const auto task_time : task_by_increasing_end_max_) { const int current_task = task_time.task_index; @@ -624,7 +629,7 @@ bool DisjunctiveOverloadChecker::PropagateSubwindow( const IntegerValue window_start = window_[critical_event].time; const IntegerValue window_end = theta_tree_.GetEnvelopeOf(critical_event) - 1; - for (int event = critical_event; event < window_size; event++) { + for (int event = critical_event; event < num_events; event++) { const IntegerValue energy_min = theta_tree_.EnergyMin(event); if (energy_min > 0) { const int task = window_[event].task_index; @@ -657,7 +662,7 @@ bool DisjunctiveOverloadChecker::PropagateSubwindow( const IntegerValue window_start = window_[critical_event].time; const IntegerValue window_end = current_end + optional_size_min - available_energy - 1; - for (int event = critical_event; event < window_size; event++) { + for (int event = critical_event; event < num_events; event++) { const IntegerValue energy_min = theta_tree_.EnergyMin(event); if (energy_min > 0) { const int task = window_[event].task_index; diff --git a/ortools/sat/disjunctive.h b/ortools/sat/disjunctive.h index 471fa92543..42f39bcf0e 100644 --- a/ortools/sat/disjunctive.h +++ b/ortools/sat/disjunctive.h @@ -135,23 +135,25 @@ class TaskSet { class DisjunctiveOverloadChecker : public PropagatorInterface { public: explicit DisjunctiveOverloadChecker(SchedulingConstraintHelper* helper) - : helper_(helper) { - // Resize this once and for all. - task_to_event_.resize(helper_->NumTasks()); - } + : helper_(helper), + window_(new TaskTime[helper->NumTasks()]), + task_to_event_(new int[helper->NumTasks()]) {} + bool Propagate() final; int RegisterWith(GenericLiteralWatcher* watcher); private: - bool PropagateSubwindow(IntegerValue global_window_end); + bool PropagateSubwindow(int relevat_size, IntegerValue global_window_end); SchedulingConstraintHelper* helper_; - std::vector window_; + // Size assigned at construction, stay fixed afterwards. + std::unique_ptr window_; + std::unique_ptr task_to_event_; + std::vector task_by_increasing_end_max_; ThetaLambdaTree theta_tree_; - std::vector task_to_event_; }; class DisjunctiveDetectablePrecedences : public PropagatorInterface { diff --git a/ortools/sat/docs/README.md b/ortools/sat/docs/README.md index b19a30b139..aed756a545 100644 --- a/ortools/sat/docs/README.md +++ b/ortools/sat/docs/README.md @@ -1,4 +1,4 @@ -| [home](README.md) | [boolean logic](boolean_logic.md) | [integer arithmetic](integer_arithmetic.md) | [channeling constraints](channeling.md) | [scheduling](scheduling.md) | [Using the CP-SAT solver](solver.md) | [Model manipulation](model.md) | [Python API](https://or-tools.github.io/docs/python/ortools/sat/python/cp_model.html) | +| [home](README.md) | [boolean logic](boolean_logic.md) | [integer arithmetic](integer_arithmetic.md) | [channeling constraints](channeling.md) | [scheduling](scheduling.md) | [Using the CP-SAT solver](solver.md) | [Model manipulation](model.md) | [Python API](https://google.github.io/or-tools/python/ortools/sat/python/cp_model.html) | | ----------------- | --------------------------------- | ------------------------------------------- | --------------------------------------- | --------------------------- | ------------------------------------ | ------------------------------ | -------------------------------- | # Using the CP-SAT solver @@ -6,18 +6,7 @@ https://developers.google.com/optimization/cp/cp_solver - -* [Using the CP-SAT solver](#using-the-cp-sat-solver) - * [Documentation structure](#documentation-structure) - * [Searching for one (optimal) solution of a CP-SAT model](#searching-for-one-optimal-solution-of-a-cp-sat-model) - * [Python code samples](#python-code-samples) - * [C++ code samples](#c-code-samples) - * [Java code samples](#java-code-samples) - * [C# code samples](#c-code-samples-1) - - - ## Documentation structure diff --git a/ortools/sat/docs/boolean_logic.md b/ortools/sat/docs/boolean_logic.md index f11b9428cc..ab2b938f65 100644 --- a/ortools/sat/docs/boolean_logic.md +++ b/ortools/sat/docs/boolean_logic.md @@ -6,30 +6,7 @@ https://developers.google.com/optimization/ - -* [Boolean logic recipes for the CP-SAT solver.](#boolean-logic-recipes-for-the-cp-sat-solver) - * [Introduction](#introduction) - * [Boolean variables and literals](#boolean-variables-and-literals) - * [Python code](#python-code) - * [C++ code](#c-code) - * [Java code](#java-code) - * [C# code](#c-code-1) - * [Boolean constraints](#boolean-constraints) - * [Python code](#python-code-1) - * [C++ code](#c-code-2) - * [Java code](#java-code-1) - * [C# code](#c-code-3) - * [Reified constraints](#reified-constraints) - * [Python code](#python-code-2) - * [C++ code](#c-code-4) - * [Java code](#java-code-2) - * [C# code](#c-code-5) - * [Product of two Boolean Variables](#product-of-two-boolean-variables) - * [Python code](#python-code-3) - - - ## Introduction diff --git a/ortools/sat/docs/channeling.md b/ortools/sat/docs/channeling.md index 14232de3ce..a00049c049 100644 --- a/ortools/sat/docs/channeling.md +++ b/ortools/sat/docs/channeling.md @@ -5,22 +5,7 @@ https://developers.google.com/optimization/ - -* [Channeling constraints](#channeling-constraints) - * [If-Then-Else expressions](#if-then-else-expressions) - * [Python code](#python-code) - * [C++ code](#c-code) - * [Java code](#java-code) - * [C# code](#c-code-1) - * [A bin-packing problem](#a-bin-packing-problem) - * [Python code](#python-code-1) - * [C++ code](#c-code-2) - * [Java code](#java-code-1) - * [C# code](#c-code-3) - - - A *channeling constraint* links variables inside a model. They're used when you want to express a complicated relationship between variables, such as "if this diff --git a/ortools/sat/docs/integer_arithmetic.md b/ortools/sat/docs/integer_arithmetic.md index 4e20323c6c..63c9ecdc47 100644 --- a/ortools/sat/docs/integer_arithmetic.md +++ b/ortools/sat/docs/integer_arithmetic.md @@ -5,38 +5,7 @@ https://developers.google.com/optimization/ - -* [Integer arithmetic recipes for the CP-SAT solver.](#integer-arithmetic-recipes-for-the-cp-sat-solver) - * [Introduction](#introduction) - * [Integer variables](#integer-variables) - * [Interval domain](#interval-domain) - * [Non-contiguous domain](#non-contiguous-domain) - * [Boolean variables](#boolean-variables) - * [Other methods](#other-methods) - * [Linear constraints](#linear-constraints) - * [C++ and Java linear constraints and linear expressions](#c-and-java-linear-constraints-and-linear-expressions) - * [Python and C# linear constraints and linear expressions](#python-and-c-linear-constraints-and-linear-expressions) - * [Generic linear constraint](#generic-linear-constraint) - * [Limitations](#limitations) - * [Rabbits and Pheasants examples](#rabbits-and-pheasants-examples) - * [Python code](#python-code) - * [C++ code](#c-code) - * [Java code](#java-code) - * [C# code](#c-code-1) - * [Earliness-Tardiness cost function.](#earliness-tardiness-cost-function) - * [Python code](#python-code-1) - * [C++ code](#c-code-2) - * [Java code](#java-code-1) - * [C# code](#c-code-3) - * [Step function.](#step-function) - * [Python code](#python-code-2) - * [C++ code](#c-code-4) - * [Java code](#java-code-2) - * [C# code](#c-code-5) - - - ## Introduction @@ -764,8 +733,6 @@ step_function_sample_sat() ```cpp #include -#include - #include "absl/types/span.h" #include "ortools/base/logging.h" #include "ortools/sat/cp_model.h" diff --git a/ortools/sat/docs/model.md b/ortools/sat/docs/model.md index 5e28fe28cb..62fd2ff9ec 100644 --- a/ortools/sat/docs/model.md +++ b/ortools/sat/docs/model.md @@ -5,21 +5,7 @@ https://developers.google.com/optimization/ - -* [Model manipulation](#model-manipulation) - * [Introduction](#introduction) - * [Solution hinting](#solution-hinting) - * [Python code](#python-code) - * [C++ code](#c-code) - * [Java code](#java-code) - * [C# code](#c-code-1) - * [Model copy](#model-copy) - * [Python code](#python-code-1) - * [C++ code](#c-code-2) - - - ## Introduction diff --git a/ortools/sat/docs/scheduling.md b/ortools/sat/docs/scheduling.md index 4060d488d3..aacd2233ee 100644 --- a/ortools/sat/docs/scheduling.md +++ b/ortools/sat/docs/scheduling.md @@ -5,43 +5,7 @@ https://developers.google.com/optimization/ - -* [Scheduling recipes for the CP-SAT solver.](#scheduling-recipes-for-the-cp-sat-solver) - * [Introduction](#introduction) - * [Interval variables](#interval-variables) - * [Python code](#python-code) - * [C++ code](#c-code) - * [Java code](#java-code) - * [C# code](#c-code-1) - * [Optional intervals](#optional-intervals) - * [Python code](#python-code-1) - * [C++ code](#c-code-2) - * [Java code](#java-code-1) - * [C# code](#c-code-3) - * [NoOverlap constraint](#nooverlap-constraint) - * [Python code](#python-code-2) - * [C++ code](#c-code-4) - * [Java code](#java-code-2) - * [C# code](#c-code-5) - * [Cumulative constraint](#cumulative-constraint) - * [Alternative resources for one interval](#alternative-resources-for-one-interval) - * [Ranking tasks in a disjunctive resource](#ranking-tasks-in-a-disjunctive-resource) - * [Python code](#python-code-3) - * [C++ code](#c-code-6) - * [Java code](#java-code-3) - * [C# code](#c-code-7) - * [Intervals spanning over breaks in the calendar](#intervals-spanning-over-breaks-in-the-calendar) - * [Python code](#python-code-4) - * [Detecting if two intervals overlap.](#detecting-if-two-intervals-overlap) - * [Python code](#python-code-5) - * [Transitions in a disjunctive resource](#transitions-in-a-disjunctive-resource) - * [Precedences between intervals](#precedences-between-intervals) - * [Convex hull of a set of intervals](#convex-hull-of-a-set-of-intervals) - * [Reservoir constraint](#reservoir-constraint) - - - ## Introduction @@ -877,9 +841,9 @@ void RankingSampleSat() { const int kHorizon = 100; const int kNumTasks = 4; - auto add_task_ranking = [&cp_model](const std::vector& starts, - const std::vector& presences, - const std::vector& ranks) { + auto add_task_ranking = [&cp_model](absl::Span starts, + absl::Span presences, + absl::Span ranks) { const int num_tasks = starts.size(); // Creates precedence variables between pairs of intervals. @@ -1354,24 +1318,33 @@ need to take into account the case where no task is performed. #!/usr/bin/env python3 """Code sample to demonstrates how to rank intervals using a circuit.""" +from typing import List, Sequence + + from ortools.sat.python import cp_model -def rank_tasks_with_circuit(model, starts, durations, presences, ranks): - """This method uses a circuit constraints to rank tasks. +def rank_tasks_with_circuit( + model: cp_model.CpModel, + starts: Sequence[cp_model.IntVar], + durations: Sequence[int], + presences: Sequence[cp_model.IntVar], + ranks: Sequence[cp_model.IntVar], +): + """This method uses a circuit constraint to rank tasks. This method assumes that all starts are disjoint, meaning that all tasks have a strictly positive duration, and they appear in the same NoOverlap constraint. To implement this ranking, we will create a dense graph with num_tasks + 1 - node. - The extra node (with id 0) will be used to decide which tasks is first with + nodes. + The extra node (with id 0) will be used to decide which task is first with its only outgoing arc, and whhich task is last with its only incoming arc. Each task i will be associated with id i + 1, and an arc between i + 1 and j + 1 indicates that j is the immediate successor of i. - The circuit constraint will insure there is at most 1 hamiltonian path of + The circuit constraint ensures there is at most 1 hamiltonian path of length > 1. If no such path exists, then no tasks are active. The multiple enforced linear constraints are meant to ensure the compatibility @@ -1388,7 +1361,7 @@ def rank_tasks_with_circuit(model, starts, durations, presences, ranks): num_tasks = len(starts) all_tasks = range(num_tasks) - arcs = [] + arcs: List[cp_model.ArcT] = [] for i in all_tasks: # if node i is first. start_lit = model.NewBoolVar(f"start_{i}") @@ -1414,8 +1387,8 @@ def rank_tasks_with_circuit(model, starts, durations, presences, ranks): # To perform the transitive reduction from precedences to successors, # we need to tie the starts of the tasks with 'literal'. - # In a pure problem, the following inequation could be an equality. - # In is not true in general. + # In a pure problem, the following inequality could be an equality. + # It is not true in general. # # Note that we could use this literal to penalize the transition, add an # extra delay to the precedence. diff --git a/ortools/sat/docs/solver.md b/ortools/sat/docs/solver.md index 9d1d70509d..aaf6f991c7 100644 --- a/ortools/sat/docs/solver.md +++ b/ortools/sat/docs/solver.md @@ -5,32 +5,7 @@ https://developers.google.com/optimization/ - -* [Solving a CP-SAT model](#solving-a-cp-sat-model) - * [Changing the parameters of the solver](#changing-the-parameters-of-the-solver) - * [Specifying the time limit in Python](#specifying-the-time-limit-in-python) - * [Specifying the time limit in C++](#specifying-the-time-limit-in-c) - * [Specifying the time limit in Java](#specifying-the-time-limit-in-java) - * [Specifying the time limit in C#.](#specifying-the-time-limit-in-c-1) - * [Printing intermediate solutions](#printing-intermediate-solutions) - * [Python code](#python-code) - * [C++ code](#c-code) - * [Java code](#java-code) - * [C# code](#c-code-1) - * [Searching for all solutions in a satisfiability model](#searching-for-all-solutions-in-a-satisfiability-model) - * [Python code](#python-code-1) - * [C++ code](#c-code-2) - * [Java code](#java-code-1) - * [C# code](#c-code-3) - * [Stopping search early](#stopping-search-early) - * [Python code](#python-code-2) - * [C++ code](#c-code-4) - * [Java code](#java-code-2) - * [C# code](#c-code-5) - - - ## Changing the parameters of the solver diff --git a/ortools/sat/feasibility_jump.cc b/ortools/sat/feasibility_jump.cc index 9f227518db..08bad19de5 100644 --- a/ortools/sat/feasibility_jump.cc +++ b/ortools/sat/feasibility_jump.cc @@ -397,6 +397,7 @@ std::function FeasibilityJumpSolver::GenerateTask(int64_t /*task_id*/) { if (should_recompute_violations) { evaluator_->ComputeAllViolations(); + move_->Clear(); } if (reset_weights) { // Each time we reset the weight, we randomly choose if we do decay or @@ -404,9 +405,18 @@ std::function FeasibilityJumpSolver::GenerateTask(int64_t /*task_id*/) { bump_value_ = 1.0; weights_.assign(evaluator_->NumEvaluatorConstraints(), 1.0); use_decay_ = absl::Bernoulli(random_, 0.5); - } - if (params_.violation_ls_use_compound_moves()) { - use_compound_moves_ = absl::Bernoulli(random_, 0.25); + move_->Clear(); + if (params_.violation_ls_use_compound_moves()) { + use_compound_moves_ = absl::Bernoulli(random_, 0.25); + if (use_compound_moves_) { + compound_weight_changed_.clear(); + in_compound_weight_changed_.assign(weights_.size(), false); + compound_weights_.assign(weights_.size(), kCompoundDiscount); + for (int c = 0; c < evaluator_->NumEvaluatorConstraints(); ++c) { + if (evaluator_->IsViolated(c)) compound_weights_[c] = weights_[c]; + } + } + } } // Search for feasible solution. ++num_batches_; @@ -537,11 +547,6 @@ void FeasibilityJumpSolver::RecomputeJump(int var) { jump_deltas_[var] = best_jump.first - solution[var]; jump_scores_[var] = best_jump.second; } - - if (IsGood(var) && !in_good_jumps_[var]) { - in_good_jumps_[var] = true; - good_jumps_.push_back(var); - } } void FeasibilityJumpSolver::RecomputeAllJumps() { @@ -555,6 +560,10 @@ void FeasibilityJumpSolver::RecomputeAllJumps() { for (int var = 0; var < num_variables; ++var) { RecomputeJump(var); + if (IsGood(var)) { + in_good_jumps_[var] = true; + good_jumps_.push_back(var); + } } } @@ -639,19 +648,24 @@ void FeasibilityJumpSolver::UpdateViolatedConstraintWeights() { } // Important: This is for debugging, but unfortunately it currently change the -// deterministic time and hence the overall algo behavior. -// -// TODO(user): Because we keep updating the score incrementally and we might -// have large constraint weight, we might have a pretty bad precision on -// the score though, so it is possible this fail. +// deterministic time and so the overall algo behavior. bool FeasibilityJumpSolver::JumpIsUpToDate(int var) { + // With decay, we have a wide magnitude of weight, so the incremental update + // of the floating point score can be widely off. Here we just check that + // without decay (integer weight) the update should be reasonably precise. + if (use_decay_) return true; + const int64_t old_jump_delta = jump_deltas_[var]; const double old_score = jump_scores_[var]; RecomputeJump(var); CHECK_EQ(jump_deltas_[var], old_jump_delta); // No change const double relative = std::max({std::abs(jump_scores_[var]), std::abs(old_score), 1.0}); - return std::abs(jump_scores_[var] - old_score) / relative < 1e-2; + const bool result = std::abs(jump_scores_[var] - old_score) / relative < 1e-2; + + // To keep the same behavior as in -c opt, we restore the old value. + jump_scores_[var] = old_score; + return result; } bool FeasibilityJumpSolver::DoSomeLinearIterations() { @@ -726,12 +740,8 @@ bool FeasibilityJumpSolver::DoSomeLinearIterations() { } if (good_jumps_.empty()) break; - DCHECK_EQ(best_score, - evaluator_->LinearEvaluator().WeightedViolationDelta( - weights_, best_var, best_delta)); - - CHECK_NE(best_var, -1); - CHECK_NE(best_index, -1); + DCHECK_NE(best_var, -1); + DCHECK_NE(best_index, -1); // Perform the move. ++num_linear_moves_; @@ -790,7 +800,7 @@ void FeasibilityJumpSolver::MarkJumpsThatNeedsToBeRecomputed(int changed_var) { // We don't need to recompute score of binary variable, it should // already be correct. if (var_has_two_values_[var]) { - DCHECK(JumpIsUpToDate(var)); + DCHECK(JumpIsUpToDate(var)) << var << " " << jump_scores_[var]; if (IsGood(var) && !in_good_jumps_[var]) { in_good_jumps_[var] = true; good_jumps_.push_back(var); @@ -816,45 +826,25 @@ bool FeasibilityJumpSolver::DoSomeGeneralIterations() { evaluator_->UpdateAllNonLinearViolations(); evaluator_->RecomputeViolatedList(/*linear_only=*/false); RecomputeVarsToScan(); - - move_->Clear(); - absl::Span scan_weights = absl::MakeConstSpan(weights_); if (use_compound_moves_) { - compound_weight_changed_.clear(); - in_compound_weight_changed_.assign(weights_.size(), false); - compound_weights_.assign(weights_.begin(), weights_.end()); - for (int c = 0; c < evaluator_->NumEvaluatorConstraints(); ++c) { - if (!evaluator_->IsViolated(c)) compound_weights_[c] *= kCompoundDiscount; - } - scan_weights = absl::MakeConstSpan(compound_weights_); + compound_move_beam_search_width_ = absl::Bernoulli(random_, 0.75) ? 1 : 2; } auto effort = [&]() { return num_general_evals_ + num_weight_updates_ + num_general_moves_; }; const int64_t effort_limit = effort() + 100000; - // Check size to make sure we are at a local minimum when we terminate. - while (move_->Size() > 0 || effort() < effort_limit) { + while (effort() < effort_limit) { int var; int64_t value; double score; bool time_limit_crossed = false; DCHECK(!move_->IsImproving()); - // If we are past the effort limit stop looking for new moves. const bool found_move = - (effort() < effort_limit && - ScanRelevantVariables(scan_weights, &var, &value, &score, - &time_limit_crossed)); + ScanRelevantVariables(&var, &value, &score, &time_limit_crossed); const bool backtrack = !found_move && move_->Backtrack(&var, &value, &score); if (found_move || backtrack) { const int64_t prev_value = solution[var]; - // Score is wrong if we are using compound moves, recompute. - if (use_compound_moves_ && !backtrack) { - ++num_general_evals_; - score = evaluator_->WeightedViolationDelta(weights_, var, - value - prev_value); - } - // Perform the move. num_general_moves_++; @@ -877,7 +867,7 @@ bool FeasibilityJumpSolver::DoSomeGeneralIterations() { evaluator_->last_update_violation_changes()) { if (violation_delta == 0) continue; for (const int var : evaluator_->ConstraintToVars(c)) { - if (violation_delta > 0 && evaluator_->IsViolated(c)) { + if (violation_delta == evaluator_->Violation(c)) { num_violated_constraints_per_var_[var] += 1; } else if (violation_delta < 0 && !evaluator_->IsViolated(c)) { num_violated_constraints_per_var_[var] -= 1; @@ -901,6 +891,7 @@ bool FeasibilityJumpSolver::DoSomeGeneralIterations() { AddVarToScan(v); } } + DCHECK(SlowCheckNumViolatedConstraints()); if (use_compound_moves_ && !backtrack) { // Make sure we can undo the move. @@ -921,6 +912,7 @@ bool FeasibilityJumpSolver::DoSomeGeneralIterations() { if (evaluator_->ViolatedConstraints().empty()) return true; UpdateViolatedConstraintWeights(); + DCHECK(SlowCheckNumViolatedConstraints()); // Constraints with increased weight may lead to new negative score moves. for (const int c : evaluator_->ViolatedConstraints()) { for (const int v : evaluator_->ConstraintToVars(c)) { @@ -948,13 +940,36 @@ void FeasibilityJumpSolver::ResetChangedCompoundWeights() { compound_weight_changed_.clear(); } -bool FeasibilityJumpSolver::ScanRelevantVariables( - absl::Span weights, int* improving_var, - int64_t* improving_value, double* improving_score, - bool* time_limit_crossed) { +bool FeasibilityJumpSolver::ShouldAcceptUnitMove(int var, int64_t delta, + double score) { + if (score < 0) return true; + return score == 0 && evaluator_->ObjectiveDelta(var, delta) < 0; +} + +bool FeasibilityJumpSolver::ShouldExtendCompoundMove(int var, int64_t delta, + double score, + double novelty) { + if (move_->Score() + score - std::max(novelty, 0.0) < 0) { + return true; + } + if (move_->Score() + score == 0 && + move_->ObjectiveDelta() + evaluator_->ObjectiveDelta(var, delta) < 0) { + return true; + } + return score < move_->BestChildScore(); +} + +bool FeasibilityJumpSolver::ScanRelevantVariables(int* improving_var, + int64_t* improving_value, + double* improving_score, + bool* time_limit_crossed) { + if (move_->NumChildrenExplored() >= compound_move_beam_search_width_) { + return false; + } DCHECK_GE(move_->Score(), 0); absl::Span solution = evaluator_->current_solution(); - + absl::Span scan_weights = + use_compound_moves_ ? compound_weights_ : weights_; while (!vars_to_scan_.empty()) { const int idx = absl::Uniform(random_, 0, vars_to_scan_.size()); const int var = vars_to_scan_[idx]; @@ -965,21 +980,23 @@ bool FeasibilityJumpSolver::ScanRelevantVariables( if (!ShouldScan(var)) continue; int64_t new_value; - double score; + double scan_score; const int64_t current_value = solution[var]; if (!use_compound_moves_ && evaluator_->VariableOnlyInLinearConstraintWithConvexViolationChange( var)) { // We lazily update the jump value. if (jump_need_recomputation_[var]) { + // We don't use good_jumps_ and in_good_jumps_ here, so we don't update + // them. RecomputeJump(var); } else { DCHECK(JumpIsUpToDate(var)); } new_value = current_value + jump_deltas_[var]; - score = jump_scores_[var]; + scan_score = jump_scores_[var]; } else { - std::tie(new_value, score) = FindBestValue( + std::tie(new_value, scan_score) = FindBestValue( var_domains_[var], current_value, [&](int64_t value) -> double { // Check the time limit periodically. if (++num_general_evals_ % 1000 == 0 && @@ -988,16 +1005,24 @@ bool FeasibilityJumpSolver::ScanRelevantVariables( *time_limit_crossed = true; } if (*time_limit_crossed) return 0.0; - return evaluator_->WeightedViolationDelta(weights, var, + return evaluator_->WeightedViolationDelta(scan_weights, var, value - current_value); }); - score += move_->Score(); } if (*time_limit_crossed) return false; - if (score > 0) continue; - const double obj_delta = - evaluator_->ObjectiveDelta(var, new_value - current_value); - if (score == 0 && obj_delta >= 0) continue; + const int64_t delta = new_value - current_value; + if (delta == 0) continue; + const double score = use_compound_moves_ + ? evaluator_->WeightedViolationDelta( + weights_, var, new_value - current_value) + : scan_score; + if (use_compound_moves_) { + if (!ShouldExtendCompoundMove(var, delta, score, score - scan_score)) { + continue; + } + } else { + if (!ShouldAcceptUnitMove(var, delta, score)) continue; + } *improving_var = var; *improving_value = new_value; *improving_score = score; @@ -1038,6 +1063,20 @@ void FeasibilityJumpSolver::RecomputeVarsToScan() { } } +bool FeasibilityJumpSolver::SlowCheckNumViolatedConstraints() const { + std::vector result; + result.assign(var_domains_.size(), 0); + for (const int c : evaluator_->ViolatedConstraints()) { + for (const int v : evaluator_->ConstraintToVars(c)) { + ++result[v]; + } + } + for (int v = 0; v < result.size(); ++v) { + CHECK_EQ(result[v], num_violated_constraints_per_var_[v]); + } + return true; +} + bool CompoundMoveBuilder::IsImproving() const { return Score() < 0 || (Score() == 0 && ObjectiveDelta() < 0); } @@ -1070,6 +1109,11 @@ void CompoundMoveBuilder::Push(int var, int64_t prev_value, double score) { const double obj_delta = evaluator_->ObjectiveDelta( var, evaluator_->current_solution()[var] - prev_value); DCHECK(!var_on_stack_[var]); + if (!stack_.empty()) { + stack_.back().best_child_score = + std::min(stack_.back().best_child_score, score); + ++stack_.back().num_children; + } var_on_stack_[var] = true; stack_.push_back( {.var = var, diff --git a/ortools/sat/feasibility_jump.h b/ortools/sat/feasibility_jump.h index d177c2eb84..87e3111310 100644 --- a/ortools/sat/feasibility_jump.h +++ b/ortools/sat/feasibility_jump.h @@ -114,8 +114,7 @@ class FeasibilityJumpSolver : public SubSolver { bool DoSomeGeneralIterations(); // Returns true if an improving move was found. - bool ScanRelevantVariables(absl::Span weights, - int* improving_var, int64_t* improving_value, + bool ScanRelevantVariables(int* improving_var, int64_t* improving_value, double* improving_score, bool* time_limit_crossed); // Increases the weight of the currently infeasible constraints. @@ -140,6 +139,20 @@ class FeasibilityJumpSolver : public SubSolver { // weights_[c] otherwise. void ResetChangedCompoundWeights(); + // Returns true if we should accept a single variable move. + bool ShouldAcceptUnitMove(int var, int64_t delta, double score); + + // Returns true if we should push this change to move_. + // `novelty` is the total discount applied to the score due to using + // `cumulative_weights_`, should always be positive (modulo floating-point + // errors). + bool ShouldExtendCompoundMove(int var, int64_t delta, double score, + double novelty); + + // Validates each element in num_violated_constraints_per_var_ against + // evaluator_->ViolatedConstraints. + bool SlowCheckNumViolatedConstraints() const; + const LinearModel* linear_model_; SatParameters params_; ModelSharedTimeLimit* shared_time_limit_; @@ -199,12 +212,16 @@ class FeasibilityJumpSolver : public SubSolver { std::vector tmp_breakpoints_; - // Each time we reset the weight, we might randomly change this to update - // them with decay or not. + // Each time we reset the weights, randomly change this to update them with + // decay or not. bool use_decay_ = true; - // Each batch, randomly decide if we will use compound moves or not. + // Each time we reset the weights, randomly decide if we will use compound + // moves or not. bool use_compound_moves_ = false; + // We randomize the beam width each batch, choosing either 1 or 2, biased + // towards 1. + int compound_move_beam_search_width_ = 1; // Statistics int64_t num_batches_ = 0; @@ -269,6 +286,20 @@ class CompoundMoveBuilder { return stack_.empty() ? 0.0 : stack_.back().cumulative_objective_delta; } + // Returns the minimum score of any child of the current node in the stack, + // with a maximum of 0. + // NB: We assume all improving moves will immediately be committed, so this + // always returns 0 when the stack is empty. + double BestChildScore() const { + return stack_.empty() ? 0.0 : stack_.back().best_child_score; + } + + // Returns the number of times a child has been pushed to the current node. + // NB: Always returns 0 when the stack is empty. + int NumChildrenExplored() const { + return stack_.empty() ? 0 : stack_.back().num_children; + } + // Returns true if this move reduces weighted violation, or improves the // objective without increasing violation. bool IsImproving() const; @@ -286,6 +317,12 @@ class CompoundMoveBuilder { // to avoid numerical issues causing negative scores after backtracking. double cumulative_score; double cumulative_objective_delta; + + // Used to avoid infinite loops, this tracks the best score of any immediate + // children (and not deeper descendants) to avoid re-exploring the same + // prefix. + double best_child_score = 0.0; + int num_children = 0; }; LsEvaluator* evaluator_; std::vector var_on_stack_; diff --git a/ortools/sat/intervals.cc b/ortools/sat/intervals.cc index ecb8c7eeb2..57dcabcb92 100644 --- a/ortools/sat/intervals.cc +++ b/ortools/sat/intervals.cc @@ -259,13 +259,13 @@ bool SchedulingConstraintHelper::UpdateCachedValues(int t) { cached_size_min_[t] = dmin; // Note that we use the cached value here for EndMin()/StartMax(). - const IntegerValue new_shifted_start_min = EndMin(t) - dmin; + const IntegerValue new_shifted_start_min = emin - dmin; if (new_shifted_start_min != cached_shifted_start_min_[t]) { recompute_energy_profile_ = true; recompute_shifted_start_min_ = true; cached_shifted_start_min_[t] = new_shifted_start_min; } - const IntegerValue new_negated_shifted_end_max = -(StartMax(t) + dmin); + const IntegerValue new_negated_shifted_end_max = -(smax + dmin); if (new_negated_shifted_end_max != cached_negated_shifted_end_max_[t]) { recompute_negated_shifted_end_max_ = true; cached_negated_shifted_end_max_[t] = new_negated_shifted_end_max; diff --git a/ortools/sat/sat_parameters.proto b/ortools/sat/sat_parameters.proto index da4da241db..5109a7dd75 100644 --- a/ortools/sat/sat_parameters.proto +++ b/ortools/sat/sat_parameters.proto @@ -23,7 +23,7 @@ option csharp_namespace = "Google.OrTools.Sat"; // Contains the definitions for all the sat algorithm parameters and their // default values. // -// NEXT TAG: 261 +// NEXT TAG: 263 message SatParameters { // In some context, like in a portfolio of search, it makes sense to name a // given parameters set for logging purpose. @@ -1276,6 +1276,10 @@ message SatParameters { // The value of the variable is specific to each strategy. optional int64 search_random_variable_pool_size = 104 [default = 0]; + // Experimental code: specify if the objective pushes all tasks toward the + // start of the schedule. + optional bool push_all_tasks_toward_start = 262 [default = true]; + // If true, we automatically detect variables whose constraint are always // enforced by the same literal and we mark them as optional. This allows // to propagate them as if they were present in some situation. @@ -1432,4 +1436,10 @@ message SatParameters { // Any value in the input mip with a magnitude lower than this will be set to // zero. This is to avoid some issue in LP presolving. optional double mip_drop_tolerance = 232 [default = 1e-16]; + + // When solving a MIP, we do some basic floating point presolving before + // scaling the problem to integer to be handled by CP-SAT. This control how + // much of that presolve we do. It can help to better scale floating point + // model, but it is not always behaving nicely. + optional int32 mip_presolve_level = 261 [default = 2]; } diff --git a/ortools/sat/synchronization.cc b/ortools/sat/synchronization.cc index 50b3e8d7ad..7be1bcdc24 100644 --- a/ortools/sat/synchronization.cc +++ b/ortools/sat/synchronization.cc @@ -273,7 +273,7 @@ void SharedResponseManager::TestGapLimitsIfNeeded() { // Note(user): Some code path in single-thread assumes that the problem // can only be solved when they have proven infeasibility and do not check // the ProblemIsSolved() method. So we force a stop here. - shared_time_limit_->Stop(); + if (always_synchronize_) shared_time_limit_->Stop(); } if (gap / std::max(1.0, std::abs(user_best)) < relative_gap_limit_) { SOLVER_LOG(logger_, "Relative gap limit of ", relative_gap_limit_, @@ -281,7 +281,7 @@ void SharedResponseManager::TestGapLimitsIfNeeded() { UpdateBestStatus(CpSolverStatus::OPTIMAL); // Same as above. - shared_time_limit_->Stop(); + if (always_synchronize_) shared_time_limit_->Stop(); } }