diff --git a/ortools/sat/disjunctive.h b/ortools/sat/disjunctive.h index 7c5e3cf340..6760b3680d 100644 --- a/ortools/sat/disjunctive.h +++ b/ortools/sat/disjunctive.h @@ -16,6 +16,7 @@ #include #include +#include #include #include "ortools/base/macros.h" diff --git a/ortools/sat/docs/troubleshooting.md b/ortools/sat/docs/troubleshooting.md index 04dc6e36eb..5720abe66f 100644 --- a/ortools/sat/docs/troubleshooting.md +++ b/ortools/sat/docs/troubleshooting.md @@ -11,17 +11,16 @@ A good explanation of the log can be found in the ## Improving performance with multiple workers -CP-SAT is built with parallelism in mind. While single worker is a good -solution, the solver is tuned to reached its best performance with multiple -workers. +CP-SAT is built with parallelism in mind. While you can achieve a good solution +with a single worker, you'll get the best results when using multiple workers. -We distinguish multiple thresholds: +There are several tiers of behavior: -- **[8 workers]** This is the minimal setting for parallel search. It blends - workers with different linear relaxations (none, default, maximal), core - based search if applicable, a quick_restart subsolver, a feasibility_jump - first solution subsolver, and dedicated Large Neighborhood Search subsolvers - to find improving solutions. +- **[8 workers]** This is the minimum number of workers needed to trigger + parallel search. It blends workers with different linear relaxations (none, + default, maximal), core-based search if applicable, a quick_restart + subsolver, a feasibility_jump first solution subsolver, and dedicated Large + Neighborhood Search subsolvers to find improved solutions. - **[16 workers]** Bumping to 16 workers adds a continuous probing subsolver, more first solution subsolvers (random search and feasibility_jump), two dual subsolvers dedicated to improving the lower bound of the objective @@ -46,11 +45,11 @@ Solving a problem yields the following possible status (CpSolverStatus): solutions of a feasibility problem (if asked). - **[INFEASIBLE]** The problem has been proven infeasible. - **[OPTIMAL]** An optimal feasible solution has been found. More generally, - this status represent a success. So we also return OPTIMAL if we find a + this status represents a success. So we also return OPTIMAL if we find a solution for a pure feasibility problem or if a gap limit has been specified and we return a solution within this limit. In the case where we need to return all the feasible solution, this status will only be returned if - solve() enumerated all of them; If we stopped before, we will return + solve() enumerated all of them; If we stopped earlier, we will return FEASIBLE. ## Debugging infeasible models @@ -58,21 +57,21 @@ Solving a problem yields the following possible status (CpSolverStatus): Oftentimes, solving a model yields an infeasibility status. While this can be a valid answer, it can happen because of either a data issue, or a model issue. -Here a simple common sense processes to help diagnose the source of the -infeasibility. +Here are some ways to diagnose the source of the infeasibility. -- reduce the model size if this is a parameter -- remove all constraints while keeping the problem infeasible -- play with the domain of the variables to enlarge them +- reduce the model size, if possible +- remove all constraints, which making sure that the resulting model is still + infeasible +- enlarge the domains of variables - inject a known feasible solution and try to find where it breaks -- add assumptions to check the soundness of your data (capacity is >= 0, no - crazy ranges of values, ...) +- check the soundness of your data (e.g., do you have negative capacities, or + unreasonable ranges?) ## Using assumptions to explain infeasibility -CP-SAT implements a system of assumptions to find the root cause of an -infeasible problem. To use it, one must instrument constraints with literals and -add these literals to the set of assumptions. +CP-SAT implements a mechanism of assumptions to find the root cause of an +infeasible problem. To use it, one must add enforcement literals to constraints +and add these literals to the set of assumptions. Note that this set is minimized but not guaranteed to be minimal. To find a minimal unsatisfiable set (MUS), you must minimize the (weighted) sum of these diff --git a/ortools/sat/feasibility_jump.cc b/ortools/sat/feasibility_jump.cc index b32ec5074e..8ece56ff64 100644 --- a/ortools/sat/feasibility_jump.cc +++ b/ortools/sat/feasibility_jump.cc @@ -408,6 +408,8 @@ std::function FeasibilityJumpSolver::GenerateTask(int64_t /*task_id*/) { if (ub < lb) return; // Search is finished. if (evaluator_->ReduceObjectiveBounds(lb.value(), ub.value())) { should_recompute_violations = true; + // The scores in the current move may now be wrong. + move_->Clear(); } } @@ -437,6 +439,10 @@ std::function FeasibilityJumpSolver::GenerateTask(int64_t /*task_id*/) { should_recompute_violations = true; } } + // Check if compound move search might backtrack out of the new domains. + if (!move_->StackValuesInDomains(var_domains_)) { + move_->Clear(); + } } if (should_recompute_violations) { @@ -450,17 +456,15 @@ std::function FeasibilityJumpSolver::GenerateTask(int64_t /*task_id*/) { weights_.assign(evaluator_->NumEvaluatorConstraints(), 1.0); use_decay_ = absl::Bernoulli(random_, 0.5); move_->Clear(); - if (params_.violation_ls_use_compound_moves()) { - use_compound_moves_ = absl::Bernoulli(random_, 0.25); - if (use_compound_moves_) { - compound_move_max_discrepancy_ = 0; - 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]; - } - } + use_compound_moves_ = absl::Bernoulli( + random_, params_.violation_ls_compound_move_probability()); + if (use_compound_moves_) { + compound_move_max_discrepancy_ = 0; + compound_weight_changed_.clear(); + in_compound_weight_changed_.assign(weights_.size(), false); + // Compound weights for violated constraints will be set to the + // correct values in InitializeCompoundWeights. + compound_weights_.assign(weights_.size(), kCompoundDiscount); } } // Search for feasible solution. @@ -495,15 +499,14 @@ std::function FeasibilityJumpSolver::GenerateTask(int64_t /*task_id*/) { }; } -double FeasibilityJumpSolver::ComputeScore( - absl::Span scan_weights, int var, int64_t delta, - bool linear_only) const { +double FeasibilityJumpSolver::ComputeScore(absl::Span weights, + int var, int64_t delta, + bool linear_only) const { constexpr double kEpsilon = 1.0 / std::numeric_limits::max(); - double score = evaluator_->LinearEvaluator().WeightedViolationDelta( - scan_weights, var, delta); + double score = + evaluator_->LinearEvaluator().WeightedViolationDelta(weights, var, delta); if (!linear_only) { - score += - evaluator_->WeightedNonLinearViolationDelta(scan_weights, var, delta); + score += evaluator_->WeightedNonLinearViolationDelta(weights, var, delta); } score += kEpsilon * evaluator_->ObjectiveDelta(var, delta); return score; @@ -524,8 +527,8 @@ std::pair FeasibilityJumpSolver::ComputeLinearJump(int var) { const int64_t max_value = var_domains_[var].Max(); const int64_t delta = solution[var] == min_value ? max_value - min_value : min_value - max_value; - return std::make_pair(delta, - ComputeScore(var, delta, /*linear_only=*/true)); + return std::make_pair( + delta, ComputeScore(ScanWeights(), var, delta, /*linear_only=*/true)); } // In practice, after a few iterations, the chance of finding an improving @@ -537,10 +540,10 @@ std::pair FeasibilityJumpSolver::ComputeLinearJump(int var) { const int64_t p2 = var_domains_[var].ValueAtOrAfter(solution[var] + 1); std::pair best_jump; - const double v1 = - var_domains_[var].Contains(p1) - ? ComputeScore(var, p1 - solution[var], /*linear_only=*/true) - : std::numeric_limits::infinity(); + const double v1 = var_domains_[var].Contains(p1) + ? ComputeScore(ScanWeights(), var, p1 - solution[var], + /*linear_only=*/true) + : std::numeric_limits::infinity(); if (v1 < 0.0) { // Point p1 is improving. Look for best before it. // Note that we can exclude all point after solution[var] since it is @@ -555,15 +558,15 @@ std::pair FeasibilityJumpSolver::ComputeLinearJump(int var) { best_jump = ConvexMinimum( /*is_to_the_right=*/true, {p1, v1}, tmp_breakpoints_, [this, var, &solution](int64_t jump_value) { - return ComputeScore(var, jump_value - solution[var], + return ComputeScore(ScanWeights(), var, jump_value - solution[var], /*linear_only=*/true); }); } } else { - const double v2 = - var_domains_[var].Contains(p2) - ? ComputeScore(var, p2 - solution[var], /*linear_only=*/true) - : std::numeric_limits::infinity(); + const double v2 = var_domains_[var].Contains(p2) + ? ComputeScore(ScanWeights(), var, p2 - solution[var], + /*linear_only=*/true) + : std::numeric_limits::infinity(); if (v2 < 0.0) { // Point p2 is improving. Look for best after it. // Similarly, we exclude the other points by convexity. @@ -577,7 +580,8 @@ std::pair FeasibilityJumpSolver::ComputeLinearJump(int var) { best_jump = ConvexMinimum( /*is_to_the_right=*/false, {p2, v2}, tmp_breakpoints_, [this, var, &solution](int64_t jump_value) { - return ComputeScore(var, jump_value - solution[var], + return ComputeScore(ScanWeights(), var, + jump_value - solution[var], /*linear_only=*/true); }); } @@ -611,14 +615,14 @@ std::pair FeasibilityJumpSolver::ComputeGeneralJump(int var) { std::pair result = RangeConvexMinimum( domain[0].start - current_value, domain[0].end - current_value + 1, [&](int64_t delta) -> double { - return ComputeScore(var, delta, /*linear_only=*/false); + return ComputeScore(ScanWeights(), var, delta, /*linear_only=*/false); }); for (int i = 1; i < domain.NumIntervals(); ++i) { const int64_t min_delta = domain[i].start - current_value; const int64_t max_delta = domain[i].end - current_value; const auto& [delta, score] = RangeConvexMinimum( result, min_delta, max_delta + 1, [&](int64_t delta) -> double { - return ComputeScore(var, delta, /*linear_only=*/false); + return ComputeScore(ScanWeights(), var, delta, /*linear_only=*/false); }); if (score < result.second) result = std::make_pair(delta, score); } @@ -632,7 +636,7 @@ void FeasibilityJumpSolver::UpdateViolatedConstraintWeights(JumpTable& jumps) { ++num_weight_updates_; // Because we update the weight incrementally, it is better to not have a - // super hight magnitude, otherwise doing +max_weight and then -max_weight + // super high magnitude, otherwise doing +max_weight and then -max_weight // will just ignore any constraint with a small weight and our // DCHECK(JumpIsUpToDate(var)) will fail more often. const double kMaxWeight = 1e10; @@ -647,6 +651,7 @@ void FeasibilityJumpSolver::UpdateViolatedConstraintWeights(JumpTable& jumps) { bool rescale = false; for (const int c : evaluator_->ViolatedConstraints()) { DCHECK(evaluator_->IsViolated(c)); + if (use_compound_moves_) DCHECK_EQ(compound_weights_[c], weights_[c]); weights_[c] += bump_value_; if (use_compound_moves_) compound_weights_[c] = weights_[c]; if (weights_[c] > kMaxWeight) rescale = true; @@ -674,8 +679,7 @@ void FeasibilityJumpSolver::UpdateViolatedConstraintWeights(JumpTable& jumps) { linear_evaluator->ClearAffectedVariables(); for_weight_update_.resize(num_variables); for (const int c : evaluator_->ViolatedConstraints()) { - // TODO(user): We can probably handle compound moves too. - if (c < evaluator_->NumLinearConstraints() && !use_compound_moves_) { + if (c < evaluator_->NumLinearConstraints()) { linear_evaluator->UpdateScoreOnWeightUpdate( c, jumps.Deltas(), absl::MakeSpan(for_weight_update_)); } else { @@ -694,8 +698,7 @@ void FeasibilityJumpSolver::UpdateViolatedConstraintWeights(JumpTable& jumps) { // TODO(user): We could compute the minimal bump that would lead to a // good move. That might change depending on the jump value though, so // we can only do that easily for Booleans. - // TODO(user): We can probably handle compound moves too. - if (!var_has_two_values_[var] || use_compound_moves_) { + if (!var_has_two_values_[var]) { jumps.Recompute(var); } else { // We may know the correct score for binary vars. @@ -780,9 +783,7 @@ bool FeasibilityJumpSolver::DoSomeLinearIterations() { void FeasibilityJumpSolver::MarkJumpsThatNeedToBeRecomputed(int changed_var, JumpTable& jumps) { for (const int var : evaluator_->VariablesAffectedByLastLinearUpdate()) { - // TODO(user): We can probably handle compound moves too. - if (var != changed_var && - (!var_has_two_values_[var] || use_compound_moves_)) { + if (var != changed_var && !var_has_two_values_[var]) { jumps.Recompute(var); } AddVarToScan(jumps, var); @@ -791,7 +792,9 @@ void FeasibilityJumpSolver::MarkJumpsThatNeedToBeRecomputed(int changed_var, evaluator_->last_update_violation_changes()) { if (c < evaluator_->NumLinearConstraints()) continue; for (const int var : evaluator_->ConstraintToVars(c)) { - if (var != changed_var) jumps.Recompute(var); + if (var != changed_var) { + jumps.Recompute(var); + } AddVarToScan(jumps, var); } } @@ -806,6 +809,7 @@ bool FeasibilityJumpSolver::DoSomeGeneralIterations() { evaluator_->UpdateAllNonLinearViolations(); evaluator_->RecomputeViolatedList(/*linear_only=*/false); RecomputeVarsToScan(general_jumps_); + InitializeCompoundWeights(); auto effort = [&]() { return num_general_evals_ + num_linear_evals_ + num_weight_updates_ + num_general_moves_; @@ -826,13 +830,12 @@ bool FeasibilityJumpSolver::DoSomeGeneralIterations() { const int64_t prev_value = solution[var]; DCHECK_NE(prev_value, value); // Update the linear part. - evaluator_->UpdateLinearScores(var, value, weights_, + evaluator_->UpdateLinearScores(var, value, ScanWeights(), general_jumps_.Deltas(), general_jumps_.MutableScores()); // Update the non-linear part. Note it also commits the move. evaluator_->UpdateNonLinearViolations(var, value); evaluator_->UpdateVariableValue(var, value); - // TODO(user): Handle compound moves? if (use_compound_moves_ && !backtrack) { // `!backtrack` is just an optimisation - we can never break any new // constraints on backtrack, so we can never change any @@ -840,7 +843,8 @@ bool FeasibilityJumpSolver::DoSomeGeneralIterations() { for (const auto& [c, violation_delta] : evaluator_->last_update_violation_changes()) { if (violation_delta == 0) continue; - if (violation_delta > 0 && compound_weights_[c] != weights_[c]) { + if (evaluator_->IsViolated(c) && + compound_weights_[c] != weights_[c]) { compound_weights_[c] = weights_[c]; if (!in_compound_weight_changed_[c]) { in_compound_weight_changed_[c] = true; @@ -850,16 +854,21 @@ bool FeasibilityJumpSolver::DoSomeGeneralIterations() { general_jumps_.Recompute(v); // Vars will be added in MarkJumpsThatNeedToBeRecomputed. } - } else if (violation_delta < 0 && !in_compound_weight_changed_[c] && - !evaluator_->IsViolated(c)) { - DCHECK_EQ(compound_weights_[c], weights_[c]); + } else if (!evaluator_->IsViolated(c) && + !in_compound_weight_changed_[c]) { in_compound_weight_changed_[c] = true; compound_weight_changed_.push_back(c); } } } if (!use_decay_) { - DCHECK_EQ(-score, ComputeScore(var, prev_value - value, false)); + // Check that the score for undoing the move is -score with both the + // default weights (which may be `weights_` or `compound_weights_`), and + // with `weights_` explicitly. + DCHECK_EQ(-score, + ComputeScore(ScanWeights(), var, prev_value - value, false)); + DCHECK_EQ(-score, + ComputeScore(weights_, var, prev_value - value, false)); } if (var_has_two_values_[var]) { // We already know the score of the only possible move (undoing what we @@ -923,6 +932,7 @@ bool FeasibilityJumpSolver::ScanRelevantVariables(int num_to_scan, int* best_var, int64_t* best_value, double* best_score) { + if (time_limit_crossed_) return false; if (move_->Discrepancy() > compound_move_max_discrepancy_) { return false; } @@ -956,15 +966,17 @@ bool FeasibilityJumpSolver::ScanRelevantVariables(int num_to_scan, time_limit_crossed_ = true; return false; } - if (time_limit_crossed_) return false; const int64_t current_value = evaluator_->current_solution()[var]; - DCHECK(var_domains_[var].Contains(current_value + delta)); + DCHECK(var_domains_[var].Contains(current_value + delta)) + << var << " " << current_value << "+" << delta << " not in " + << var_domains_[var].ToString(); DCHECK(!var_domains_[var].IsFixed()); // Note that this will likely fail if you use decaying weights as they // will have large magnitudes and the incremental update will be // imprecise. DCHECK(use_decay_ || jumps.JumpIsUpToDate(var)) - << var << " " << var_domains_[var].ToString(); + << var << " " << var_domains_[var].ToString() << " " + << ComputeScore(ScanWeights(), var, delta, (&jumps == &linear_jumps_)); if (scan_score >= 0) { remove_var_to_scan_at_index(index); continue; @@ -1031,6 +1043,13 @@ void FeasibilityJumpSolver::RecomputeVarsToScan(JumpTable& jumps) { } } +void FeasibilityJumpSolver::InitializeCompoundWeights() { + if (!use_compound_moves_) return; + for (const int c : evaluator_->ViolatedConstraints()) { + compound_weights_[c] = weights_[c]; + } +} + bool FeasibilityJumpSolver::SlowCheckNumViolatedConstraints() const { std::vector result; result.assign(var_domains_.size(), 0); @@ -1087,4 +1106,13 @@ void CompoundMoveBuilder::Push(int var, int64_t prev_value, double score) { .discrepancy = Discrepancy(), }); } + +bool CompoundMoveBuilder::StackValuesInDomains( + absl::Span var_domains) const { + for (const UnitMove& move : stack_) { + if (!var_domains[move.var].Contains(move.prev_value)) return false; + } + return true; +} + } // namespace operations_research::sat diff --git a/ortools/sat/feasibility_jump.h b/ortools/sat/feasibility_jump.h index d618c37b98..a754c5e8fa 100644 --- a/ortools/sat/feasibility_jump.h +++ b/ortools/sat/feasibility_jump.h @@ -158,16 +158,18 @@ class FeasibilityJumpSolver : public SubSolver { void PerturbateCurrentSolution(); std::string OneLineStats() const; - // Returns the weighted violation delta plus epsilon * the objective delta. - double ComputeScore(absl::Span scan_weights, int var, - int64_t delta, bool linear_only) const; - - // As above, but uses appropriate scan weights based on the current value of - // `use_compound_moves_`. - double ComputeScore(int var, int64_t delta, bool linear_only) const { - return ComputeScore(use_compound_moves_ ? compound_weights_ : weights_, var, - delta, linear_only); + absl::Span ScanWeights() { + return absl::MakeSpan(use_compound_moves_ ? compound_weights_ : weights_); } + absl::Span ScanWeights() const { + return absl::MakeConstSpan(use_compound_moves_ ? compound_weights_ + : weights_); + } + + // Returns the weighted violation delta plus epsilon * the objective delta. + double ComputeScore(absl::Span weights, int var, int64_t delta, + bool linear_only) const; + // Computes the optimal value for variable v, considering only the violation // of linear constraints. std::pair ComputeLinearJump(int var); @@ -195,6 +197,12 @@ class FeasibilityJumpSolver : public SubSolver { void UpdateNumViolatedConstraintsPerVar(); void RecomputeVarsToScan(JumpTable&); + + // Ensures that all currently violated constraints have compound_weight_[c] == + // weight_[c]. Mostly only necessary for the first batch with new weights or a + // new imported solution or if the objective bounds get tightened. + void InitializeCompoundWeights(); + // Returns true if it is possible that `var` may have value that reduces // weighted violation or improve the objective. // Note that this is independent of the actual weights used. @@ -351,6 +359,9 @@ class CompoundMoveBuilder { // Returns the number of backtracking moves that have been applied. int NumBacktracks() const { return num_backtracks_; } + // Returns true if all prev_values on the stack are in the appropriate domain. + bool StackValuesInDomains(absl::Span var_domains) const; + private: struct UnitMove { int var; diff --git a/ortools/sat/parameters_validation.cc b/ortools/sat/parameters_validation.cc index 9fe93e14cb..ab114956ef 100644 --- a/ortools/sat/parameters_validation.cc +++ b/ortools/sat/parameters_validation.cc @@ -116,8 +116,10 @@ std::string ValidateParameters(const SatParameters& params) { TEST_IN_RANGE(feasibility_jump_var_perburbation_range_ratio, 0.0, 1.0); // Violation ls. + TEST_NOT_NAN(violation_ls_compound_move_probability); TEST_IN_RANGE(num_violation_ls, 0, kMaxReasonableParallelism); TEST_IN_RANGE(violation_ls_perturbation_period, 1, 1'000'000'000); + TEST_IN_RANGE(violation_ls_compound_move_probability, 0.0, 1.0); TEST_POSITIVE(glucose_decay_increment_period); TEST_POSITIVE(shared_tree_max_nodes_per_worker); diff --git a/ortools/sat/sat_parameters.proto b/ortools/sat/sat_parameters.proto index 8e5ec22416..19cdb7bbeb 100644 --- a/ortools/sat/sat_parameters.proto +++ b/ortools/sat/sat_parameters.proto @@ -1121,8 +1121,9 @@ message SatParameters { // How long violation_ls should wait before perturbating a solution. optional int32 violation_ls_perturbation_period = 249 [default = 100]; - // Search for ejection-chain-like moves in jump and violation_ls workers. - optional bool violation_ls_use_compound_moves = 259 [default = true]; + // Probability of using compound move search each restart. + // TODO(user): Add reference to paper when published. + optional double violation_ls_compound_move_probability = 259 [default = 0.5]; // Enables experimental workstealing-like shared tree search. // If non-zero, start this many complete worker threads to explore a shared