always-cautious/ortools-clone
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

[SAT] new search options: fixed_restart, randomize_search, lns!

Browse Source
This commit is contained in:
Laurent Perron
2018-06-15 11:12:00 +02:00
parent 27286f42ae
commit 653de024ef
8 changed files with 778 additions and 153 deletions
Download Patch File Download Diff File Expand all files Collapse all files

201
ortools/sat/cp_model_lns.cc Normal file
View File

@@ -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 <unordered_set>
#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<int> 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<int> 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<int> random_var(0, num_variables - 1);
std::vector<bool> used_variables(num_variables, false);
std::vector<int> relaxed_variables;
relaxed_variables.push_back(random_var(random));
used_variables[relaxed_variables.back()] = true;
std::vector<int> 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<int> 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<bool> relaxed_variables(var_to_constraint_.size(), false);
std::vector<bool> added_constraints(constraint_to_var_.size(), false);
std::vector<int> next_constraints;
// Start by a random constraint.
{
std::uniform_int_distribution<int> random_start(
0, constraint_to_var_.size() - 1);
next_constraints.push_back(random_start(random));
added_constraints[next_constraints.back()] = true;
}
std::vector<int> 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<int> 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<int> 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

92
ortools/sat/cp_model_lns.h Normal file
View File

@@ -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 <vector>
#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<std::vector<int>> constraint_to_var_;
std::vector<std::vector<int>> 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_

210
ortools/sat/cp_model_search.cc
View File

@@ -13,9 +13,11 @@
#include "ortools/sat/cp_model_search.h"
#include <random>
#include <unordered_map>
#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<LiteralIndex()> ConstructSearchStrategyInternal(
const std::unordered_map<int, std::pair<int64, int64>>&
var_to_coeff_offset_pair,
@@ -38,87 +53,126 @@ const std::function<LiteralIndex()> 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<SatParameters>();
// 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<VarValue> 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<SearchStrategyRandomGenerator>() != nullptr;
SearchStrategyRandomGenerator* const rw =
model->GetOrCreate<SearchStrategyRandomGenerator>();
if (!seeded) {
rw->rand.seed(parameters->random_seed());
}
const int winner = std::uniform_int_distribution<int>(
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<LiteralIndex()> ConstructSearchStrategy(

269
ortools/sat/cp_model_solver.cc
View File

@@ -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<SatParameters>()->set_max_number_of_conflicts(10);
std::vector<BooleanOrIntegerVariable> 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<SatParameters>()->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<DratProofHandler>(
/*in_binary_format=*/false, output, FLAGS_drat_check);
} else {
drat_proof_handler.reset(new DratProofHandler());
drat_proof_handler = absl::make_unique<DratProofHandler>();
}
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<SatParameters>();
const auto& observers = model->GetOrCreate<SolutionObservers>()->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<int64>(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<int64>(response.solution().begin(),
response.solution().end())));
}
return response;
}
// Do the actual presolve.
CpModelProto mapping_proto;
std::vector<int> 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, &copy);
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<SatParameters>();
parameters->set_stop_after_first_solution(true);
CpSolverResponse response;
{
CpModelProto copy = expanded_model;
response = PresolveAndSolve(model_proto, &copy, 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<std::unique_ptr<NeighborhoodGenerator>> generators;
generators.push_back(
absl::make_unique<SimpleNeighborhoodGenerator>(&expanded_model));
generators.push_back(
absl::make_unique<VariableGraphNeighborhoodGenerator>(&expanded_model));
generators.push_back(
absl::make_unique<ConstraintGraphNeighborhoodGenerator>(&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<AdaptiveParameterValue> difficulties(generators.size(),
AdaptiveParameterValue(0.5));
TimeLimit* limit = model->GetOrCreate<TimeLimit>();
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<int64>(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<SatParameters>();
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<SolutionObservers>()->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<int64>(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<int64>(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<int> 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, &copy);
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

14
ortools/sat/integer_search.cc
View File

@@ -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);

97
ortools/sat/lns.h Normal file
View File

@@ -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 <functional>
#include <vector>
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<bool()>
// - SolveNeighborhoodFunction: std::function<std::function<void()>(int seed)>
template <class StopFunction, class SolveNeighborhoodFunction>
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 <class StopFunction, class SolveNeighborhoodFunction>
inline void OptimizeWithLNS(
int num_threads, StopFunction stop_function,
SolveNeighborhoodFunction generate_and_solve_function) {
int64 seed = 0;
while (!stop_function()) {
std::function<void()> update_function = generate_and_solve_function(seed++);
update_function();
}
}
} // namespace sat
} // namespace operations_research
#endif // OR_TOOLS_SAT_LNS_H_

16
ortools/sat/restart.cc
View File

@@ -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;
}

32
ortools/sat/sat_parameters.proto
View File

@@ -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];
}