diff --git a/ortools/sat/cp_model_lns.cc b/ortools/sat/cp_model_lns.cc new file mode 100644 index 0000000000..5f1ebe0561 --- /dev/null +++ b/ortools/sat/cp_model_lns.cc @@ -0,0 +1,201 @@ +// Copyright 2010-2017 Google +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#include "ortools/sat/cp_model_lns.h" + +#include +#include "ortools/sat/cp_model_utils.h" +#include "ortools/util/random_engine.h" + +namespace operations_research { +namespace sat { + +// Base class constructor. +NeighborhoodGenerator::NeighborhoodGenerator(CpModelProto const* model) + : model_proto_(*model) { + var_to_constraint_.resize(model_proto_.variables_size()); + constraint_to_var_.resize(model_proto_.constraints_size()); + for (int ct_index = 0; ct_index < model_proto_.constraints_size(); + ++ct_index) { + for (const int var : UsedVariables(model_proto_.constraints(ct_index))) { + var_to_constraint_[var].push_back(ct_index); + constraint_to_var_[ct_index].push_back(var); + CHECK_GE(var, 0); + CHECK_LT(var, model_proto_.variables_size()); + } + } +} + +namespace { + +// Returns a CpModelProto where the variables at given position where fixed to +// the value they take in the given response. +CpModelProto FixGivenPosition(const CpSolverResponse& response, + const std::vector variables_to_fix, + CpModelProto model_proto) { + CHECK_EQ(response.solution_size(), model_proto.variables_size()); + for (const int var : variables_to_fix) { + model_proto.mutable_variables(var)->clear_domain(); + model_proto.mutable_variables(var)->add_domain(response.solution(var)); + model_proto.mutable_variables(var)->add_domain(response.solution(var)); + } + + // Set the current solution as a hint. + model_proto.clear_solution_hint(); + for (int var = 0; var < model_proto.variables_size(); ++var) { + model_proto.mutable_solution_hint()->add_vars(var); + model_proto.mutable_solution_hint()->add_values(response.solution(var)); + } + + // TODO(user): force better objective? Note that this is already done when the + // hint above is sucessfully loaded (i.e. if it passes the presolve correctly) + // since the solver will try to find better solution than the current one. + return model_proto; +} + +} // namespace + +CpModelProto SimpleNeighborhoodGenerator::Generate( + const CpSolverResponse& initial_solution, int64 seed, + double difficulty) const { + random_engine_t random; + random.seed(seed); + + // TODO(user): we could generate this more efficiently that using random + // shuffle. + const int num_vars = model_proto_.variables_size(); + std::vector fixed_variables(num_vars); + for (int i = 0; i < num_vars; ++i) fixed_variables[i] = i; + std::shuffle(fixed_variables.begin(), fixed_variables.end(), random); + + const int target_size = std::ceil((1.0 - difficulty) * num_vars); + fixed_variables.resize(target_size); + return FixGivenPosition(initial_solution, fixed_variables, model_proto_); +} + +CpModelProto VariableGraphNeighborhoodGenerator::Generate( + const CpSolverResponse& initial_solution, int64 seed, + double difficulty) const { + const int num_variables = model_proto_.variables_size(); + const int target_size = std::ceil(difficulty * num_variables); + if (target_size == num_variables) return model_proto_; + CHECK_GT(target_size, 0); + + random_engine_t random; + random.seed(seed); + + std::uniform_int_distribution random_var(0, num_variables - 1); + std::vector used_variables(num_variables, false); + std::vector relaxed_variables; + + relaxed_variables.push_back(random_var(random)); + used_variables[relaxed_variables.back()] = true; + + std::vector random_variables; + for (int i = 0; i < relaxed_variables.size(); ++i) { + random_variables.clear(); + // Collect all the variables that appears in the same constraints as + // relaxed_variables[i]. + for (const int ct : var_to_constraint_[relaxed_variables[i]]) { + for (const int var : constraint_to_var_[ct]) { + if (used_variables[var]) continue; + used_variables[var] = true; + random_variables.push_back(var); + } + } + if (relaxed_variables.size() + random_variables.size() > target_size) { + random_variables.resize(target_size - relaxed_variables.size()); + } + std::shuffle(random_variables.begin(), random_variables.end(), random); + relaxed_variables.insert(relaxed_variables.end(), random_variables.begin(), + random_variables.end()); + if (relaxed_variables.size() >= target_size) break; + } + + // Compute the complement of relaxed_variables in fixed_variables. + std::vector fixed_variables; + { + used_variables.assign(num_variables, false); + for (const int var : relaxed_variables) used_variables[var] = true; + for (int var = 0; var < num_variables; ++var) { + if (!used_variables[var]) fixed_variables.push_back(var); + } + } + return FixGivenPosition(initial_solution, fixed_variables, model_proto_); +} + +CpModelProto ConstraintGraphNeighborhoodGenerator::Generate( + const CpSolverResponse& initial_solution, int64 seed, + double difficulty) const { + const int num_vars = model_proto_.variables_size(); + const int target_size = std::ceil(difficulty * num_vars); + if (target_size == num_vars) return model_proto_; + CHECK_GT(target_size, 0); + + random_engine_t random; + random.seed(seed); + + int num_used_variables = 0; + std::vector relaxed_variables(var_to_constraint_.size(), false); + std::vector added_constraints(constraint_to_var_.size(), false); + std::vector next_constraints; + + // Start by a random constraint. + { + std::uniform_int_distribution random_start( + 0, constraint_to_var_.size() - 1); + next_constraints.push_back(random_start(random)); + added_constraints[next_constraints.back()] = true; + } + + std::vector random_variables; + while (num_used_variables < target_size) { + // Stop if we have a full connected component. + if (next_constraints.empty()) break; + + // Pick a random unprocessed constraint. + std::uniform_int_distribution random_constraint( + 0, next_constraints.size() - 1); + const int i = random_constraint(random); + const int contraint_index = next_constraints[i]; + std::swap(next_constraints[i], next_constraints.back()); + next_constraints.pop_back(); + + // Add all the variable of this constraint and increase the set of next + // possible constraints. + CHECK_LT(contraint_index, constraint_to_var_.size()); + random_variables = constraint_to_var_[contraint_index]; + std::shuffle(random_variables.begin(), random_variables.end(), random); + for (const int var : random_variables) { + if (!relaxed_variables[var]) { + ++num_used_variables; + } + relaxed_variables[var] = true; + if (num_used_variables == target_size) break; + for (const int ct : var_to_constraint_[var]) { + if (added_constraints[ct]) continue; + added_constraints[ct] = true; + next_constraints.push_back(ct); + } + } + } + + std::vector fixed_variables; + for (int i = 0; i < num_vars; ++i) { + if (!relaxed_variables[i]) fixed_variables.push_back(i); + } + return FixGivenPosition(initial_solution, fixed_variables, model_proto_); +} + +} // namespace sat +} // namespace operations_research diff --git a/ortools/sat/cp_model_lns.h b/ortools/sat/cp_model_lns.h new file mode 100644 index 0000000000..11eb7c998e --- /dev/null +++ b/ortools/sat/cp_model_lns.h @@ -0,0 +1,92 @@ +// Copyright 2010-2017 Google +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#ifndef OR_TOOLS_SAT_CP_MODEL_LNS_H_ +#define OR_TOOLS_SAT_CP_MODEL_LNS_H_ + +#include + +#include "ortools/base/integral_types.h" +#include "ortools/sat/cp_model.pb.h" + +namespace operations_research { +namespace sat { + +// Base class for a CpModelProto neighborhood generator. +class NeighborhoodGenerator { + public: + explicit NeighborhoodGenerator(CpModelProto const* model); + virtual ~NeighborhoodGenerator() {} + + // Generates a "local" subproblem for the given seed. + // + // The difficulty will be in [0, 1] and is related to the asked neighborhood + // size (and thus local problem difficulty). A difficulty of 0.0 means empty + // neighborhood and a difficulty of 1.0 means the full problem. The algorithm + // should try to generate a neighborhood according to this difficulty which + // will be dynamically adjusted depending on whether or not we can solve the + // subproblem in a given time limit. + // + // The given initial_solution should contains a feasible solution to the + // initial CpModelProto given to this class. Any solution to the returned + // CPModelProto should also be valid solution to the same initial model. + // + // This function should be thread-safe. + virtual CpModelProto Generate(const CpSolverResponse& initial_solution, + int64 seed, double difficulty) const = 0; + + protected: + // TODO(user): for now and for convenience, we generate the + // variable-constraint graph even if not all subclass will need it. + const CpModelProto& model_proto_; + std::vector> constraint_to_var_; + std::vector> var_to_constraint_; +}; + +// Pick a random subset of variables. +class SimpleNeighborhoodGenerator : public NeighborhoodGenerator { + public: + explicit SimpleNeighborhoodGenerator(CpModelProto const* model) + : NeighborhoodGenerator(model) {} + CpModelProto Generate(const CpSolverResponse& initial_solution, int64 seed, + double difficulty) const final; +}; + +// Pick a random subset of variables that are constructed by a BFS in the +// variable <-> constraint graph. That is, pick a random variable, then all the +// variable connected by some constraint to the first one, and so on. The +// variable of the last "level" are selected randomly. +class VariableGraphNeighborhoodGenerator : public NeighborhoodGenerator { + public: + explicit VariableGraphNeighborhoodGenerator(CpModelProto const* model) + : NeighborhoodGenerator(model) {} + CpModelProto Generate(const CpSolverResponse& initial_solution, int64 seed, + double difficulty) const final; +}; + +// Pick a random subset of constraint and relax all of their variables. We are a +// bit smarter than this because after the first contraint is selected, we only +// select constraints that share at least one variable with the already selected +// constraints. The variable from the "last" constraint are selected randomly. +class ConstraintGraphNeighborhoodGenerator : public NeighborhoodGenerator { + public: + explicit ConstraintGraphNeighborhoodGenerator(CpModelProto const* model) + : NeighborhoodGenerator(model) {} + CpModelProto Generate(const CpSolverResponse& initial_solution, int64 seed, + double difficulty) const final; +}; + +} // namespace sat +} // namespace operations_research + +#endif // OR_TOOLS_SAT_CP_MODEL_LNS_H_ diff --git a/ortools/sat/cp_model_search.cc b/ortools/sat/cp_model_search.cc index 8403a99da6..6648fea50c 100644 --- a/ortools/sat/cp_model_search.cc +++ b/ortools/sat/cp_model_search.cc @@ -13,9 +13,11 @@ #include "ortools/sat/cp_model_search.h" +#include #include #include "ortools/sat/cp_model_utils.h" +#include "ortools/util/random_engine.h" namespace operations_research { namespace sat { @@ -29,6 +31,19 @@ struct Strategy { DecisionStrategyProto::VariableSelectionStrategy var_strategy; DecisionStrategyProto::DomainReductionStrategy domain_strategy; }; + +// Stores one variable and its strategy value. +struct VarValue { + IntegerVariable var; + IntegerValue value; +}; + +// Thin wrapper around a random_engine. This helps store a random engine +// inside a Model as template programming messes up compilation. +struct SearchStrategyRandomGenerator { + random_engine_t rand; +}; + const std::function ConstructSearchStrategyInternal( const std::unordered_map>& var_to_coeff_offset_pair, @@ -38,87 +53,126 @@ const std::function ConstructSearchStrategyInternal( // Note that we copy strategies to keep the return function validity // independently of the life of the passed vector. - return - [integer_encoder, integer_trail, strategies, var_to_coeff_offset_pair]() { - for (const Strategy& strategy : strategies) { - IntegerVariable candidate = kNoIntegerVariable; - IntegerValue candidate_value = kMaxIntegerValue; - IntegerValue candidate_lb; - IntegerValue candidate_ub; + return [integer_encoder, integer_trail, strategies, var_to_coeff_offset_pair, + model]() { + const SatParameters* const parameters = model->GetOrCreate(); - // TODO(user): Improve the complexity if this becomes an issue which - // may be the case if we do a fixed_search. - for (const IntegerVariable var : strategy.variables) { - if (integer_trail->IsCurrentlyIgnored(var)) continue; - const IntegerValue lb = integer_trail->LowerBound(var); - const IntegerValue ub = integer_trail->UpperBound(var); - if (lb == ub) continue; - IntegerValue value(0); - IntegerValue coeff(1); - IntegerValue offset(0); - if (gtl::ContainsKey(var_to_coeff_offset_pair, var.value())) { - const auto coeff_offset = - gtl::FindOrDie(var_to_coeff_offset_pair, var.value()); - coeff = coeff_offset.first; - offset = coeff_offset.second; - } - DCHECK_GT(coeff, 0); - switch (strategy.var_strategy) { - case DecisionStrategyProto::CHOOSE_FIRST: - break; - case DecisionStrategyProto::CHOOSE_LOWEST_MIN: - value = coeff * lb + offset; - break; - case DecisionStrategyProto::CHOOSE_HIGHEST_MAX: - value = -(coeff * ub + offset); - break; - case DecisionStrategyProto::CHOOSE_MIN_DOMAIN_SIZE: - // TODO(user): Evaluate an exact domain computation. - value = coeff * (ub - lb + 1); - break; - case DecisionStrategyProto::CHOOSE_MAX_DOMAIN_SIZE: - // TODO(user): Evaluate an exact domain computation. - value = -coeff * (ub - lb + 1); - break; - default: - LOG(FATAL) << "Unknown VariableSelectionStrategy " - << strategy.var_strategy; - } - if (value < candidate_value) { - candidate = var; - candidate_lb = lb; - candidate_ub = ub; - candidate_value = value; - } - if (strategy.var_strategy == DecisionStrategyProto::CHOOSE_FIRST) - break; - } - if (candidate == kNoIntegerVariable) continue; + for (const Strategy& strategy : strategies) { + IntegerVariable candidate = kNoIntegerVariable; + IntegerValue candidate_value = kMaxIntegerValue; + IntegerValue candidate_lb; + IntegerValue candidate_ub; - IntegerLiteral literal; - switch (strategy.domain_strategy) { - case DecisionStrategyProto::SELECT_MIN_VALUE: - literal = IntegerLiteral::LowerOrEqual(candidate, candidate_lb); - break; - case DecisionStrategyProto::SELECT_MAX_VALUE: - literal = IntegerLiteral::GreaterOrEqual(candidate, candidate_ub); - break; - case DecisionStrategyProto::SELECT_LOWER_HALF: - literal = IntegerLiteral::LowerOrEqual( - candidate, candidate_lb + (candidate_ub - candidate_lb) / 2); - break; - case DecisionStrategyProto::SELECT_UPPER_HALF: - literal = IntegerLiteral::GreaterOrEqual( - candidate, candidate_ub - (candidate_ub - candidate_lb) / 2); - break; - default: - LOG(FATAL) << "Unknown DomainReductionStrategy " - << strategy.domain_strategy; - } - return integer_encoder->GetOrCreateAssociatedLiteral(literal).Index(); + // TODO(user): Improve the complexity if this becomes an issue which + // may be the case if we do a fixed_search. + + // To store equivalent variables in randomized search. + std::vector active_vars; + + for (const IntegerVariable var : strategy.variables) { + if (integer_trail->IsCurrentlyIgnored(var)) continue; + const IntegerValue lb = integer_trail->LowerBound(var); + const IntegerValue ub = integer_trail->UpperBound(var); + if (lb == ub) continue; + IntegerValue value(0); + IntegerValue coeff(1); + IntegerValue offset(0); + if (gtl::ContainsKey(var_to_coeff_offset_pair, var.value())) { + const auto coeff_offset = + gtl::FindOrDie(var_to_coeff_offset_pair, var.value()); + coeff = coeff_offset.first; + offset = coeff_offset.second; } - return kNoLiteralIndex; - }; + DCHECK_GT(coeff, 0); + switch (strategy.var_strategy) { + case DecisionStrategyProto::CHOOSE_FIRST: + break; + case DecisionStrategyProto::CHOOSE_LOWEST_MIN: + value = coeff * lb + offset; + break; + case DecisionStrategyProto::CHOOSE_HIGHEST_MAX: + value = -(coeff * ub + offset); + break; + case DecisionStrategyProto::CHOOSE_MIN_DOMAIN_SIZE: + // TODO(user): Evaluate an exact domain computation. + value = coeff * (ub - lb + 1); + break; + case DecisionStrategyProto::CHOOSE_MAX_DOMAIN_SIZE: + // TODO(user): Evaluate an exact domain computation. + value = -coeff * (ub - lb + 1); + break; + default: + LOG(FATAL) << "Unknown VariableSelectionStrategy " + << strategy.var_strategy; + } + if (value < candidate_value) { + candidate = var; + candidate_lb = lb; + candidate_ub = ub; + candidate_value = value; + } + if (strategy.var_strategy == DecisionStrategyProto::CHOOSE_FIRST && + !parameters->randomize_search()) { + break; + } else if (parameters->randomize_search()) { + if (active_vars.empty() || + value <= candidate_value + + parameters->search_randomization_tolerance()) { + active_vars.push_back({var, value}); + } + } + } + if (candidate == kNoIntegerVariable) continue; + if (parameters->randomize_search()) { + CHECK(!active_vars.empty()); + const IntegerValue threshold( + candidate_value + parameters->search_randomization_tolerance()); + auto is_above_tolerance = [threshold](const VarValue& entry) { + return entry.value > threshold; + }; + // Remove all values above tolerance. + active_vars.erase(std::remove_if(active_vars.begin(), active_vars.end(), + is_above_tolerance), + active_vars.end()); + // TODO(user): Move this to its own method. + const bool seeded = + model->Get() != nullptr; + SearchStrategyRandomGenerator* const rw = + model->GetOrCreate(); + if (!seeded) { + rw->rand.seed(parameters->random_seed()); + } + const int winner = std::uniform_int_distribution( + 0, active_vars.size() - 1)(rw->rand); + candidate = active_vars[winner].var; + candidate_lb = integer_trail->LowerBound(candidate); + candidate_ub = integer_trail->UpperBound(candidate); + } + + IntegerLiteral literal; + switch (strategy.domain_strategy) { + case DecisionStrategyProto::SELECT_MIN_VALUE: + literal = IntegerLiteral::LowerOrEqual(candidate, candidate_lb); + break; + case DecisionStrategyProto::SELECT_MAX_VALUE: + literal = IntegerLiteral::GreaterOrEqual(candidate, candidate_ub); + break; + case DecisionStrategyProto::SELECT_LOWER_HALF: + literal = IntegerLiteral::LowerOrEqual( + candidate, candidate_lb + (candidate_ub - candidate_lb) / 2); + break; + case DecisionStrategyProto::SELECT_UPPER_HALF: + literal = IntegerLiteral::GreaterOrEqual( + candidate, candidate_ub - (candidate_ub - candidate_lb) / 2); + break; + default: + LOG(FATAL) << "Unknown DomainReductionStrategy " + << strategy.domain_strategy; + } + return integer_encoder->GetOrCreateAssociatedLiteral(literal).Index(); + } + return kNoLiteralIndex; + }; } std::function ConstructSearchStrategy( diff --git a/ortools/sat/cp_model_solver.cc b/ortools/sat/cp_model_solver.cc index 9ae4b1e8ee..d6ab4e5cbf 100644 --- a/ortools/sat/cp_model_solver.cc +++ b/ortools/sat/cp_model_solver.cc @@ -35,6 +35,7 @@ #include "ortools/base/iterator_adaptors.h" #include "ortools/base/join.h" #include "ortools/base/map_util.h" +#include "ortools/base/memory.h" #include "ortools/base/stl_util.h" #include "ortools/graph/connectivity.h" #include "ortools/port/proto_utils.h" @@ -43,6 +44,7 @@ #include "ortools/sat/cp_constraints.h" #include "ortools/sat/cp_model_checker.h" #include "ortools/sat/cp_model_expand.h" +#include "ortools/sat/cp_model_lns.h" #include "ortools/sat/cp_model_presolve.h" #include "ortools/sat/cp_model_search.h" #include "ortools/sat/cp_model_utils.h" @@ -53,6 +55,7 @@ #include "ortools/sat/intervals.h" #include "ortools/sat/linear_programming_constraint.h" #include "ortools/sat/linear_relaxation.h" +#include "ortools/sat/lns.h" #include "ortools/sat/optimization.h" #include "ortools/sat/pb_constraint.h" #include "ortools/sat/precedences.h" @@ -2433,7 +2436,6 @@ CpSolverResponse SolveCpModelInternal( // reset the solver to its initial state, but then with phase saving it // should still follow the same path again. if (model_proto.has_solution_hint()) { - LOG(INFO) << "Loading solution hint ... "; const int64 old_conflict_limit = parameters.max_number_of_conflicts(); model->GetOrCreate()->set_max_number_of_conflicts(10); std::vector vars; @@ -2457,9 +2459,9 @@ CpSolverResponse SolveCpModelInternal( status = SolveProblemWithPortfolioSearch(decision_policies, {no_restart}, model); if (status == SatSolver::Status::MODEL_SAT) { - LOG(INFO) << "Solution hint: success, feasible solution found."; + VLOG(1) << "Solution hint: success, feasible solution found."; } else { - LOG(INFO) << "Solution: failure, no feasible solution found."; + VLOG(1) << "Solution: failure, no feasible solution found."; } model->GetOrCreate()->set_max_number_of_conflicts( old_conflict_limit); @@ -2651,10 +2653,10 @@ CpSolverResponse SolvePureSatModel(const CpModelProto& model_proto, if (!FLAGS_drat_output.empty()) { File* output; CHECK_OK(file::Open(FLAGS_drat_output, "w", &output, file::Defaults())); - drat_proof_handler.reset(new DratProofHandler( - /*in_binary_format=*/false, output, FLAGS_drat_check)); + drat_proof_handler = absl::make_unique( + /*in_binary_format=*/false, output, FLAGS_drat_check); } else { - drat_proof_handler.reset(new DratProofHandler()); + drat_proof_handler = absl::make_unique(); } solver->SetDratProofHandler(drat_proof_handler.get()); } @@ -2804,6 +2806,198 @@ CpSolverResponse SolvePureSatModel(const CpModelProto& model_proto, return response; } +// The model_proto is just used for solution checking and reconstruction of the +// solution to the original model. +CpSolverResponse PresolveAndSolve(const CpModelProto& model_proto, + CpModelProto* current_model, Model* model) { + // Solve without presolving ? + const SatParameters& params = *model->GetOrCreate(); + const auto& observers = model->GetOrCreate()->observers; + const int num_original_variables = model_proto.variables_size(); + if (!params.cp_model_presolve() || params.enumerate_all_solutions()) { + CpSolverResponse response = SolveCpModelInternal( + *current_model, true, + [&](const CpSolverResponse& intermediate_response) { + if (observers.empty()) return; + // Truncate the solution in case model expansion added more variables. + CpSolverResponse truncated_response = intermediate_response; + truncated_response.mutable_solution()->Truncate( + num_original_variables); + DCHECK(SolutionIsFeasible( + model_proto, + std::vector(truncated_response.solution().begin(), + truncated_response.solution().end()))); + for (const auto& observer : observers) { + observer(truncated_response); + } + }, + model); + if (response.status() == CpSolverStatus::MODEL_SAT || + response.status() == CpSolverStatus::OPTIMAL) { + // Truncate the solution in case model expansion added more variables. + response.mutable_solution()->Truncate(num_original_variables); + CHECK(SolutionIsFeasible(model_proto, + std::vector(response.solution().begin(), + response.solution().end()))); + } + return response; + } + + // Do the actual presolve. + CpModelProto mapping_proto; + std::vector postsolve_mapping; + PresolveCpModel(current_model, &mapping_proto, &postsolve_mapping); + VLOG(1) << CpModelStats(*current_model); + + // Note that it is okay to use the initial model_proto in the postsolve even + // though we called PresolveCpModel() on the expanded proto. This is because + // PostsolveResponse() only use the proto to known the number of variables to + // fill in the response and to check the solution feasibility of these + // variables. + CpSolverResponse response = SolveCpModelInternal( + *current_model, true, + [&](const CpSolverResponse& response) { + if (observers.empty()) return; + CpSolverResponse copy = response; + PostsolveResponse(model_proto, mapping_proto, postsolve_mapping, ©); + for (const auto& observer : observers) { + observer(copy); + } + }, + model); + PostsolveResponse(model_proto, mapping_proto, postsolve_mapping, &response); + return response; +} + +CpSolverResponse SolveCpModelWithLNS(const CpModelProto& model_proto, + const CpModelProto& expanded_model, + Model* model) { + SatParameters* parameters = model->GetOrCreate(); + parameters->set_stop_after_first_solution(true); + CpSolverResponse response; + { + CpModelProto copy = expanded_model; + response = PresolveAndSolve(model_proto, ©, model); + } + if (response.status() != CpSolverStatus::MODEL_SAT) { + return response; + } + LOG(INFO) << "LNS First solution: " << response.objective_value(); + + // For now we will just alternate between our possible neighborhoods. + // + // TODO(user): work on the presolved global problem rather than just the + // expanded problem? + std::vector> generators; + generators.push_back( + absl::make_unique(&expanded_model)); + generators.push_back( + absl::make_unique(&expanded_model)); + generators.push_back( + absl::make_unique(&expanded_model)); + + // The "optimal" difficulties do not have to be the same for different + // generators. TODO(user): move this inside the generator API? + std::vector difficulties(generators.size(), + AdaptiveParameterValue(0.5)); + + TimeLimit* limit = model->GetOrCreate(); + double deterministic_time = 0.1; + int num_no_progress = 0; + + const int num_threads = std::max(1, parameters->lns_num_threads()); + OptimizeWithLNS( + num_threads, + [&]() { + // If we didn't see any progress recently, bump the time limit. + // TODO(user): Tune the logic and expose the parameters. + if (num_no_progress > 100) { + deterministic_time *= 1.1; + num_no_progress = 0; + } + return limit->LimitReached() || + response.objective_value() == response.best_objective_bound(); + }, + [&](int64 seed) { + AdaptiveParameterValue& difficulty = + difficulties[seed % generators.size()]; + const double saved_difficulty = difficulty.value(); + CpModelProto local_problem = + generators[seed % generators.size()]->Generate(response, seed, + saved_difficulty); + + Model local_model; + { + SatParameters local_parameters; + local_parameters = *parameters; + local_parameters.set_max_deterministic_time(deterministic_time); + local_parameters.set_stop_after_first_solution(false); + local_model.Add(NewSatParameters(local_parameters)); + } + const CpSolverResponse local_response = + PresolveAndSolve(model_proto, &local_problem, &local_model); + + return [&num_no_progress, &model_proto, &response, &difficulty, + &deterministic_time, saved_difficulty, local_response]() { + // TODO(user): This is not ideal in multithread because even though + // the saved_difficulty will be the same for all thread, we will + // Increase()/Decrease() the difficuty sequentially more than once. + if (local_response.status() == CpSolverStatus::OPTIMAL || + local_response.status() == CpSolverStatus::MODEL_UNSAT) { + difficulty.Increase(); + } else { + difficulty.Decrease(); + } + if (local_response.status() == CpSolverStatus::MODEL_SAT || + local_response.status() == CpSolverStatus::OPTIMAL) { + LOG(INFO) << "dtime: " << deterministic_time + << " difficulty: " << saved_difficulty + << " local_obj: " << local_response.objective_value() + << " global_obj: " << response.objective_value(); + + // If the objective are the same, we override the solution, + // otherwise we just ignore this local solution and increment + // num_no_progress. + double coeff = model_proto.objective().scaling_factor(); + if (coeff == 0.0) coeff = 1.0; + if (local_response.objective_value() * coeff >= + response.objective_value() * coeff) { + if (local_response.objective_value() * coeff > + response.objective_value() * coeff) { + return; + } + ++num_no_progress; + } else { + num_no_progress = 0; + } + + // Update the global response. + *(response.mutable_solution()) = local_response.solution(); + response.set_objective_value(local_response.objective_value()); + response.set_wall_time(response.wall_time() + + local_response.wall_time()); + response.set_user_time(response.user_time() + + local_response.user_time()); + response.set_deterministic_time( + response.deterministic_time() + + local_response.deterministic_time()); + DCHECK(SolutionIsFeasible( + model_proto, + std::vector(local_response.solution().begin(), + local_response.solution().end()))); + } + }; + }); + + if (response.status() == CpSolverStatus::MODEL_SAT) { + if (response.objective_value() == response.best_objective_bound()) { + response.set_status(CpSolverStatus::OPTIMAL); + } + } + + return response; +} + } // namespace CpSolverResponse SolveCpModel(const CpModelProto& model_proto, Model* model) { @@ -2846,7 +3040,7 @@ CpSolverResponse SolveCpModel(const CpModelProto& model_proto, Model* model) { // automatic. const SatParameters& params = *model->GetOrCreate(); if (!model_proto.has_objective() && !model_proto.has_solution_hint() && - !params.enumerate_all_solutions()) { + !params.enumerate_all_solutions() && !params.use_lns()) { bool is_pure_sat = true; for (const IntegerVariableProto& var : model_proto.variables()) { if (var.domain_size() != 2 || var.domain(0) < 0 || var.domain(1) > 1) { @@ -2867,64 +3061,11 @@ CpSolverResponse SolveCpModel(const CpModelProto& model_proto, Model* model) { } // Starts by expanding some constraints if needed. - CpModelProto presolved_proto = ExpandCpModel(model_proto); - - // Solve without presolving ? - const auto& observers = model->GetOrCreate()->observers; - const int num_original_variables = model_proto.variables_size(); - if (!params.cp_model_presolve() || params.enumerate_all_solutions()) { - CpSolverResponse response = SolveCpModelInternal( - presolved_proto, true, - [&](const CpSolverResponse& intermediate_response) { - if (observers.empty()) return; - // Truncate the solution in case model expansion added more variables. - CpSolverResponse truncated_response = intermediate_response; - truncated_response.mutable_solution()->Truncate( - num_original_variables); - DCHECK(SolutionIsFeasible( - model_proto, - std::vector(truncated_response.solution().begin(), - truncated_response.solution().end()))); - for (const auto& observer : observers) { - observer(truncated_response); - } - }, - model); - if (response.status() == CpSolverStatus::MODEL_SAT || - response.status() == CpSolverStatus::OPTIMAL) { - // Truncate the solution in case model expansion added more variables. - response.mutable_solution()->Truncate(num_original_variables); - CHECK(SolutionIsFeasible(model_proto, - std::vector(response.solution().begin(), - response.solution().end()))); - } - return response; + CpModelProto expanded_proto = ExpandCpModel(model_proto); + if (params.use_lns() && model_proto.has_objective()) { + return SolveCpModelWithLNS(model_proto, expanded_proto, model); } - - // Do the actual presolve. - CpModelProto mapping_proto; - std::vector postsolve_mapping; - PresolveCpModel(&presolved_proto, &mapping_proto, &postsolve_mapping); - VLOG(1) << CpModelStats(presolved_proto); - - // Note that it is okay to use the initial model_proto in the postsolve even - // though we called PresolveCpModel() on the expanded proto. This is because - // PostsolveResponse() only use the proto to known the number of variables to - // fill in the response and to check the solution feasibility of these - // variables. - CpSolverResponse response = SolveCpModelInternal( - presolved_proto, true, - [&](const CpSolverResponse& response) { - if (observers.empty()) return; - CpSolverResponse copy = response; - PostsolveResponse(model_proto, mapping_proto, postsolve_mapping, ©); - for (const auto& observer : observers) { - observer(copy); - } - }, - model); - PostsolveResponse(model_proto, mapping_proto, postsolve_mapping, &response); - return response; + return PresolveAndSolve(model_proto, &expanded_proto, model); } } // namespace sat diff --git a/ortools/sat/integer_search.cc b/ortools/sat/integer_search.cc index 41cb6b1b62..93d6c0b489 100644 --- a/ortools/sat/integer_search.cc +++ b/ortools/sat/integer_search.cc @@ -280,12 +280,18 @@ SatSolver::Status SolveIntegerProblemWithLazyEncoding( {SatSolverRestartPolicy(model)}, model); } case SatParameters::FIXED_SEARCH: { - auto no_restart = []() { return false; }; // Not all Boolean might appear in next_decision(), so once there is no // decision left, we fix all Booleans that are still undecided. - return SolveProblemWithPortfolioSearch( - {SequentialSearch({next_decision, SatSolverHeuristic(model)})}, - {no_restart}, model); + if (parameters.randomize_search()) { + return SolveProblemWithPortfolioSearch( + {SequentialSearch({next_decision, SatSolverHeuristic(model)})}, + {SatSolverRestartPolicy(model)}, model); + } else { + auto no_restart = []() { return false; }; + return SolveProblemWithPortfolioSearch( + {SequentialSearch({next_decision, SatSolverHeuristic(model)})}, + {no_restart}, model); + } } case SatParameters::PORTFOLIO_SEARCH: { auto incomplete_portfolio = AddModelHeuristics({next_decision}, model); diff --git a/ortools/sat/lns.h b/ortools/sat/lns.h new file mode 100644 index 0000000000..606396b594 --- /dev/null +++ b/ortools/sat/lns.h @@ -0,0 +1,97 @@ +// Copyright 2010-2017 Google +// Licensed under the Apache License, Version 2.0 (the "License"); +// you may not use this file except in compliance with the License. +// You may obtain a copy of the License at +// +// http://www.apache.org/licenses/LICENSE-2.0 +// +// Unless required by applicable law or agreed to in writing, software +// distributed under the License is distributed on an "AS IS" BASIS, +// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +// See the License for the specific language governing permissions and +// limitations under the License. + +#ifndef OR_TOOLS_SAT_LNS_H_ +#define OR_TOOLS_SAT_LNS_H_ + +#include +#include + +namespace operations_research { +namespace sat { + +// A simple deterministic and multithreaded LNS design. +// +// While !stop_function(), we call a batch of num_threads +// generate_and_solve_function() in parallel. Each such functions should return +// an update_global_solution() function. These update functions will only be +// called sequentially (and in seed order) once the batch is done. +// +// The seed starts at zero and will be increased one by one. Each +// generate_and_solve_function() will get a different seed. +// +// Only the generate_and_solve_function(int seed) need to be thread-safe. +// +// The two templated types should behave like: +// - StopFunction: std::function +// - SolveNeighborhoodFunction: std::function(int seed)> +template +void OptimizeWithLNS(int num_threads, StopFunction stop_function, + SolveNeighborhoodFunction generate_and_solve_function); + +// Basic adaptive [0.0, 1.0] parameter that can be increased or decreased with a +// step that get smaller and smaller with the number of updates. +// +// Note(user): The current logic work well in practice, but has no theoretical +// foundation. So it might be possible to use better formulas depending on the +// situation. +class AdaptiveParameterValue { + public: + // Initial value is in [0.0, 1.0], both 0.0 and 1.0 are valid. + explicit AdaptiveParameterValue(double initial_value) + : value_(initial_value) {} + + void Reset() { num_changes_ = 0; } + + void Increase() { + const double factor = IncreaseNumChangesAndGetFactor(); + value_ = std::min(1.0 - (1.0 - value_) / factor, value_ * factor); + } + + void Decrease() { + const double factor = IncreaseNumChangesAndGetFactor(); + value_ = std::max(value_ / factor, 1.0 - (1.0 - value_) * factor); + } + + double value() const { return value_; } + + private: + // We want to change the parameters more and more slowly. + double IncreaseNumChangesAndGetFactor() { + ++num_changes_; + return 1.0 + 1.0 / std::sqrt(num_changes_ + 1); + } + + double value_; + int64 num_changes_ = 0; +}; + +// ============================================================================ +// Implementation. +// ============================================================================ + +template +inline void OptimizeWithLNS( + int num_threads, StopFunction stop_function, + SolveNeighborhoodFunction generate_and_solve_function) { + int64 seed = 0; + while (!stop_function()) { + std::function update_function = generate_and_solve_function(seed++); + update_function(); + } +} + +} // namespace sat +} // namespace operations_research + +#endif // OR_TOOLS_SAT_LNS_H_ diff --git a/ortools/sat/restart.cc b/ortools/sat/restart.cc index 2e4e6f8cee..b243d3e2c7 100644 --- a/ortools/sat/restart.cc +++ b/ortools/sat/restart.cc @@ -27,7 +27,7 @@ void RestartPolicy::Reset() { conflicts_until_next_strategy_change_ = strategy_change_conflicts_; luby_count_ = 0; - conflicts_until_next_restart_ = parameters_.luby_restart_period(); + conflicts_until_next_restart_ = parameters_.restart_period(); dl_running_average_.Reset(parameters_.restart_running_window_size()); lbd_running_average_.Reset(parameters_.restart_running_window_size()); @@ -50,6 +50,8 @@ void RestartPolicy::Reset() { tmp = SatParameters::DL_MOVING_AVERAGE_RESTART; } else if (string_value == "LBD_MOVING_AVERAGE_RESTART") { tmp = SatParameters::LBD_MOVING_AVERAGE_RESTART; + } else if (string_value == "FIXED_RESTART") { + tmp = SatParameters::FIXED_RESTART; } else { LOG(WARNING) << "Couldn't parse the RestartAlgorithm name: '" << string_value << "'."; @@ -97,6 +99,11 @@ bool RestartPolicy::ShouldRestart() { should_restart = true; } break; + case SatParameters::FIXED_RESTART: + if (conflicts_until_next_restart_ == 0) { + should_restart = true; + } + break; } if (should_restart) { num_restarts_++; @@ -113,8 +120,11 @@ bool RestartPolicy::ShouldRestart() { // Reset the various restart strategies. dl_running_average_.ClearWindow(); lbd_running_average_.ClearWindow(); - conflicts_until_next_restart_ = - parameters_.luby_restart_period() * SUniv(luby_count_ + 1); + conflicts_until_next_restart_ = parameters_.restart_period(); + if (strategies_[strategy_counter_ % strategies_.size()] == + SatParameters::LUBY_RESTART) { + conflicts_until_next_restart_ *= SUniv(luby_count_ + 1); + } } return should_restart; } diff --git a/ortools/sat/sat_parameters.proto b/ortools/sat/sat_parameters.proto index dcd0198a10..d7846a6b89 100644 --- a/ortools/sat/sat_parameters.proto +++ b/ortools/sat/sat_parameters.proto @@ -18,7 +18,7 @@ package operations_research.sat; // Contains the definitions for all the sat algorithm parameters and their // default values. // -// NEXT TAG: 100 +// NEXT TAG: 105 message SatParameters { // ========================================================================== // Branching and polarity @@ -226,7 +226,7 @@ message SatParameters { enum RestartAlgorithm { NO_RESTART = 0; - // Just follow a Luby sequence times luby_restart_period. + // Just follow a Luby sequence times restart_period. LUBY_RESTART = 1; // Moving average restart based on the decision level of conflicts. @@ -234,6 +234,9 @@ message SatParameters { // Moving average restart based on the LBD of conflicts. LBD_MOVING_AVERAGE_RESTART = 3; + + // Fixed period restart every restart period. + FIXED_RESTART = 4; } // The restart strategies will change each time the strategy_counter is @@ -250,8 +253,9 @@ message SatParameters { [default = "LUBY_RESTART,LBD_MOVING_AVERAGE_RESTART,DL_MOVING_AVERAGE_RESTART"]; - // Multiplier for the LUBY_RESTART algorithm. - optional int32 luby_restart_period = 30 [default = 50]; + // Restart period for the FIXED_RESTART strategy. This is also the multiplier + // used by the LUBY_RESTART strategy. + optional int32 restart_period = 30 [default = 50]; // Size of the window for the moving average restarts. optional int32 restart_running_window_size = 62 [default = 50]; @@ -560,4 +564,24 @@ message SatParameters { // For an optimization problem, stop the solver as soon as we have a solution. optional bool stop_after_first_solution = 98 [default = false]; + + // Specify the number of parallel workers to use during search. + // A number <= 1 means no parallelism. + optional int32 num_search_workers = 100 [default = 0]; + + // LNS parameters. + optional bool use_lns = 101 [default = false]; + optional int32 lns_num_threads = 102 [default = 1]; + + // Randomize fixed search. + optional bool randomize_search = 103 [default = false]; + + // Search randomization will collect equivalent 'max valued' variables, and + // pick one randomly. For instance, if the variable strategy is CHOOSE_FIRST, + // all unassigned variables are equivalent. If the variable strategy is + // CHOOSE_LOWEST_MIN, and `lm` is the current lowest min of all unassigned + // variables, then the set of max valued variables will be all unassigned + // variables where + // lm <= variable min <= lm + search_randomization_tolerance + optional int64 search_randomization_tolerance = 104 [default = 0]; }