internal fixes; move CP-SAT code around

This commit is contained in:
Laurent Perron
2019-04-26 16:06:25 +02:00
parent f956f2ea1c
commit e49f50d4b9
7 changed files with 249 additions and 174 deletions

View File

@@ -691,7 +691,6 @@ RoutingModel::RoutingModel(const RoutingIndexManager& index_manager,
const int64 size = Size();
index_to_pickup_index_pairs_.resize(size);
index_to_delivery_index_pairs_.resize(size);
index_to_visit_type_.resize(index_manager.num_indices(), kUnassigned);
index_to_vehicle_.resize(index_manager.num_indices(), kUnassigned);
@@ -1847,11 +1846,17 @@ void RoutingModel::CloseModelWithParameters(
nexts_, active_, vehicle_vars_, zero_transit));
}
// Nodes which are not in a disjunction are mandatory.
// Nodes which are not in a disjunction are mandatory, and those with a
// trivially infeasible type are necessarily unperformed
for (int i = 0; i < size; ++i) {
if (GetDisjunctionIndices(i).empty() && active_[i]->Max() != 0) {
active_[i]->SetValue(1);
}
const int type = GetVisitType(i);
if (type != kUnassigned &&
gtl::ContainsKey(trivially_infeasible_visit_types_, type)) {
active_[i]->SetValue(0);
}
}
// Reduce domains of vehicle variables
@@ -3708,8 +3713,14 @@ void RoutingModel::AddTemporalRequiredTypeAlternatives(
int dependent_type, absl::flat_hash_set<int> required_type_alternatives) {
DCHECK_LT(dependent_type,
temporal_required_type_alternatives_per_type_index_.size());
has_temporal_type_requirements_ = true;
if (required_type_alternatives.empty()) {
// The dependent_type requires an infeasible (empty) set of types.
trivially_infeasible_visit_types_.insert(dependent_type);
return;
}
has_temporal_type_requirements_ = true;
temporal_required_type_alternatives_per_type_index_[dependent_type].push_back(
std::move(required_type_alternatives));
}

View File

@@ -1544,6 +1544,7 @@ class RoutingModel {
std::vector<std::vector<absl::flat_hash_set<int> > >
temporal_required_type_alternatives_per_type_index_;
absl::flat_hash_set<int> trivially_infeasible_visit_types_;
bool has_temporal_type_requirements_;
// clang-format on
int num_visit_types_;

View File

@@ -55,11 +55,11 @@ class CBCInterface : public MPSolverInterface {
// ----- Parameters -----
util::Status SetNumThreads(int num_threads) override {
CHECK_GE(num_threads, 1);
num_threads_ = num_threads;
return util::OkStatus();
}
util::Status SetNumThreads(int num_threads) override {
CHECK_GE(num_threads, 1);
num_threads_ = num_threads;
return util::OkStatus();
}
// ----- Solve -----
// Solve the problem using the parameter values specified.
@@ -384,9 +384,10 @@ MPSolver::ResultStatus CBCInterface::Solve(const MPSolverParameters& param) {
// through callCbc.
model.setAllowableFractionGap(relative_mip_gap_);
// NOTE: Trailing space is required to avoid buffer overflow in cbc.
int return_status = num_threads_ == 1 ?
callCbc("-solve ", model) :
callCbc(absl::StrCat("-threads ", num_threads_, " -solve "), model);
int return_status =
num_threads_ == 1
? callCbc("-solve ", model)
: callCbc(absl::StrCat("-threads ", num_threads_, " -solve "), model);
const int kBadReturnStatus = 777;
CHECK_NE(kBadReturnStatus, return_status); // Should never happen according
// to the CBC source

View File

@@ -32,6 +32,8 @@
extern "C" {
#include "gurobi_c.h"
int __stdcall GRBisqp(GRBenv**, const char*, const char*, const char*, int,
const char*);
}
DEFINE_int32(num_gurobi_threads, 4, "Number of threads available for Gurobi.");

View File

@@ -646,169 +646,6 @@ IntegerVariable GetOrCreateVariableGreaterOrEqualToSumOf(
return new_var;
}
// Add a linear relaxation of the CP constraint to the set of linear
// constraints. The highest linearization_level is, the more types of constraint
// we encode. This method should be called only for linearization_level > 0.
//
// TODO(user): In full generality, we could encode all the constraint as an LP.
// TODO(user,user): Add unit tests for this method.
void TryToLinearizeConstraint(const CpModelProto& model_proto,
const ConstraintProto& ct, Model* m,
int linearization_level,
LinearRelaxation* relaxation) {
DCHECK_GT(linearization_level, 0);
auto* mapping = m->GetOrCreate<CpModelMapping>();
if (ct.constraint_case() == ConstraintProto::ConstraintCase::kBoolOr) {
if (linearization_level < 2) return;
LinearConstraintBuilder lc(m, IntegerValue(1), kMaxIntegerValue);
for (const int enforcement_ref : ct.enforcement_literal()) {
CHECK(lc.AddLiteralTerm(mapping->Literal(NegatedRef(enforcement_ref)),
IntegerValue(1)));
}
for (const int ref : ct.bool_or().literals()) {
CHECK(lc.AddLiteralTerm(mapping->Literal(ref), IntegerValue(1)));
}
relaxation->linear_constraints.push_back(lc.Build());
} else if (ct.constraint_case() ==
ConstraintProto::ConstraintCase::kBoolAnd) {
// TODO(user): These constraints can be many, and if they are not regrouped
// in big at most ones, then they should probably only added lazily as cuts.
// Regroup this with future clique-cut separation logic.
if (linearization_level < 2) return;
if (!HasEnforcementLiteral(ct)) return;
if (ct.enforcement_literal().size() == 1) {
const Literal enforcement = mapping->Literal(ct.enforcement_literal(0));
for (const int ref : ct.bool_and().literals()) {
relaxation->at_most_ones.push_back(
{enforcement, mapping->Literal(ref).Negated()});
}
return;
}
// Andi(e_i) => Andj(x_j)
// <=> num_rhs_terms <= Sum_j(x_j) + num_rhs_terms * Sum_i(~e_i)
int num_literals = ct.bool_and().literals_size();
LinearConstraintBuilder lc(m, IntegerValue(num_literals), kMaxIntegerValue);
for (const int ref : ct.bool_and().literals()) {
CHECK(lc.AddLiteralTerm(mapping->Literal(ref), IntegerValue(1)));
}
for (const int enforcement_ref : ct.enforcement_literal()) {
CHECK(lc.AddLiteralTerm(mapping->Literal(NegatedRef(enforcement_ref)),
IntegerValue(num_literals)));
}
relaxation->linear_constraints.push_back(lc.Build());
} else if (ct.constraint_case() ==
ConstraintProto::ConstraintCase::kAtMostOne) {
if (HasEnforcementLiteral(ct)) return;
std::vector<Literal> at_most_one;
for (const int ref : ct.at_most_one().literals()) {
at_most_one.push_back(mapping->Literal(ref));
}
relaxation->at_most_ones.push_back(at_most_one);
} else if (ct.constraint_case() == ConstraintProto::ConstraintCase::kIntMax) {
if (HasEnforcementLiteral(ct)) return;
const int target = ct.int_max().target();
for (const int var : ct.int_max().vars()) {
// This deal with the corner case X = max(X, Y, Z, ..) !
// Note that this can be presolved into X >= Y, X >= Z, ...
if (target == var) continue;
LinearConstraintBuilder lc(m, kMinIntegerValue, IntegerValue(0));
lc.AddTerm(mapping->Integer(var), IntegerValue(1));
lc.AddTerm(mapping->Integer(target), IntegerValue(-1));
relaxation->linear_constraints.push_back(lc.Build());
}
} else if (ct.constraint_case() == ConstraintProto::ConstraintCase::kIntMin) {
if (HasEnforcementLiteral(ct)) return;
const int target = ct.int_min().target();
for (const int var : ct.int_min().vars()) {
if (target == var) continue;
LinearConstraintBuilder lc(m, kMinIntegerValue, IntegerValue(0));
lc.AddTerm(mapping->Integer(target), IntegerValue(1));
lc.AddTerm(mapping->Integer(var), IntegerValue(-1));
relaxation->linear_constraints.push_back(lc.Build());
}
} else if (ct.constraint_case() ==
ConstraintProto::ConstraintCase::kIntProd) {
if (HasEnforcementLiteral(ct)) return;
const int target = ct.int_prod().target();
const int size = ct.int_prod().vars_size();
// We just linearize x = y^2 by x >= y which is far from ideal but at
// least pushes x when y moves away from zero. Note that if y is negative,
// we should probably also add x >= -y, but then this do not happen in
// our test set.
if (size == 2 && ct.int_prod().vars(0) == ct.int_prod().vars(1)) {
LinearConstraintBuilder lc(m, kMinIntegerValue, IntegerValue(0));
lc.AddTerm(mapping->Integer(ct.int_prod().vars(0)), IntegerValue(1));
lc.AddTerm(mapping->Integer(target), IntegerValue(-1));
relaxation->linear_constraints.push_back(lc.Build());
}
} else if (ct.constraint_case() == ConstraintProto::ConstraintCase::kLinear) {
AppendLinearConstraintRelaxation(ct, linearization_level, *m, relaxation);
} else if (ct.constraint_case() ==
ConstraintProto::ConstraintCase::kCircuit) {
if (HasEnforcementLiteral(ct)) return;
const int num_arcs = ct.circuit().literals_size();
CHECK_EQ(num_arcs, ct.circuit().tails_size());
CHECK_EQ(num_arcs, ct.circuit().heads_size());
// Each node must have exactly one incoming and one outgoing arc (note that
// it can be the unique self-arc of this node too).
std::map<int, std::vector<Literal>> incoming_arc_constraints;
std::map<int, std::vector<Literal>> outgoing_arc_constraints;
for (int i = 0; i < num_arcs; i++) {
const Literal arc = mapping->Literal(ct.circuit().literals(i));
const int tail = ct.circuit().tails(i);
const int head = ct.circuit().heads(i);
// Make sure this literal has a view.
m->Add(NewIntegerVariableFromLiteral(arc));
outgoing_arc_constraints[tail].push_back(arc);
incoming_arc_constraints[head].push_back(arc);
}
for (const auto* node_map :
{&outgoing_arc_constraints, &incoming_arc_constraints}) {
for (const auto& entry : *node_map) {
const std::vector<Literal>& exactly_one = entry.second;
if (exactly_one.size() > 1) {
LinearConstraintBuilder at_least_one_lc(m, IntegerValue(1),
kMaxIntegerValue);
for (const Literal l : exactly_one) {
CHECK(at_least_one_lc.AddLiteralTerm(l, IntegerValue(1)));
}
// We separate the two constraints.
relaxation->at_most_ones.push_back(exactly_one);
relaxation->linear_constraints.push_back(at_least_one_lc.Build());
}
}
}
} else if (ct.constraint_case() ==
ConstraintProto::ConstraintCase::kElement) {
const IntegerVariable index = mapping->Integer(ct.element().index());
const IntegerVariable target = mapping->Integer(ct.element().target());
const std::vector<IntegerVariable> vars =
mapping->Integers(ct.element().vars());
// We only relax the case where all the vars are constant.
// target = sum (index == i) * fixed_vars[i].
LinearConstraintBuilder constraint(m, IntegerValue(0), IntegerValue(0));
constraint.AddTerm(target, IntegerValue(-1));
IntegerTrail* integer_trail = m->GetOrCreate<IntegerTrail>();
for (const auto literal_value : m->Add(FullyEncodeVariable((index)))) {
const IntegerVariable var = vars[literal_value.value.value()];
if (!m->Get(IsFixed(var))) return;
// Make sure this literal has a view.
m->Add(NewIntegerVariableFromLiteral(literal_value.literal));
CHECK(constraint.AddLiteralTerm(literal_value.literal,
integer_trail->LowerBound(var)));
}
relaxation->linear_constraints.push_back(constraint.Build());
}
}
void TryToAddCutGenerators(const CpModelProto& model_proto,
const ConstraintProto& ct, Model* m,
LinearRelaxation* relaxation) {

View File

@@ -224,6 +224,223 @@ void AppendPartialGreaterThanEncodingRelaxation(IntegerVariable var,
}
}
// Add a linear relaxation of the CP constraint to the set of linear
// constraints. The highest linearization_level is, the more types of constraint
// we encode. This method should be called only for linearization_level > 0.
//
// TODO(user): In full generality, we could encode all the constraint as an LP.
// TODO(user,user): Add unit tests for this method.
void TryToLinearizeConstraint(const CpModelProto& model_proto,
const ConstraintProto& ct, Model* model,
int linearization_level,
LinearRelaxation* relaxation) {
DCHECK_GT(linearization_level, 0);
auto* mapping = model->GetOrCreate<CpModelMapping>();
if (ct.constraint_case() == ConstraintProto::ConstraintCase::kBoolOr) {
if (linearization_level < 2) return;
LinearConstraintBuilder lc(model, IntegerValue(1), kMaxIntegerValue);
for (const int enforcement_ref : ct.enforcement_literal()) {
CHECK(lc.AddLiteralTerm(mapping->Literal(NegatedRef(enforcement_ref)),
IntegerValue(1)));
}
for (const int ref : ct.bool_or().literals()) {
CHECK(lc.AddLiteralTerm(mapping->Literal(ref), IntegerValue(1)));
}
relaxation->linear_constraints.push_back(lc.Build());
} else if (ct.constraint_case() ==
ConstraintProto::ConstraintCase::kBoolAnd) {
// TODO(user): These constraints can be many, and if they are not regrouped
// in big at most ones, then they should probably only added lazily as cuts.
// Regroup this with future clique-cut separation logic.
if (linearization_level < 2) return;
if (!HasEnforcementLiteral(ct)) return;
if (ct.enforcement_literal().size() == 1) {
const Literal enforcement = mapping->Literal(ct.enforcement_literal(0));
for (const int ref : ct.bool_and().literals()) {
relaxation->at_most_ones.push_back(
{enforcement, mapping->Literal(ref).Negated()});
}
return;
}
// Andi(e_i) => Andj(x_j)
// <=> num_rhs_terms <= Sum_j(x_j) + num_rhs_terms * Sum_i(~e_i)
int num_literals = ct.bool_and().literals_size();
LinearConstraintBuilder lc(model, IntegerValue(num_literals),
kMaxIntegerValue);
for (const int ref : ct.bool_and().literals()) {
CHECK(lc.AddLiteralTerm(mapping->Literal(ref), IntegerValue(1)));
}
for (const int enforcement_ref : ct.enforcement_literal()) {
CHECK(lc.AddLiteralTerm(mapping->Literal(NegatedRef(enforcement_ref)),
IntegerValue(num_literals)));
}
relaxation->linear_constraints.push_back(lc.Build());
} else if (ct.constraint_case() ==
ConstraintProto::ConstraintCase::kAtMostOne) {
if (HasEnforcementLiteral(ct)) return;
std::vector<Literal> at_most_one;
for (const int ref : ct.at_most_one().literals()) {
at_most_one.push_back(mapping->Literal(ref));
}
relaxation->at_most_ones.push_back(at_most_one);
} else if (ct.constraint_case() == ConstraintProto::ConstraintCase::kIntMax) {
if (HasEnforcementLiteral(ct)) return;
// Case X = max(X_1, X_2, ..., X_N)
// Part 1: Encode X >= max(X_1, X_2, ..., X_N)
const int target = ct.int_max().target();
for (const int var : ct.int_max().vars()) {
// This deal with the corner case X = max(X, Y, Z, ..) !
// Note that this can be presolved into X >= Y, X >= Z, ...
if (target == var) continue;
LinearConstraintBuilder lc(model, kMinIntegerValue, IntegerValue(0));
lc.AddTerm(mapping->Integer(var), IntegerValue(1));
lc.AddTerm(mapping->Integer(target), IntegerValue(-1));
relaxation->linear_constraints.push_back(lc.Build());
}
// TODO(user): Check the coefficient of target in the objective to avoid
// the second part.
// Part 2: Encode upper bound on X.
// NOTE(user) for size = 2, we can do this with 1 less variable.
if (ct.int_max().vars_size() == 2 || linearization_level >= 2) {
const int positive_target = PositiveRef(target);
// TODO(user): We can get tighter bound on max target value by going
// through all the constraint vars.
const int64 max_target_value =
RefIsPositive(target)
? model_proto.variables(positive_target)
.domain(
model_proto.variables(positive_target).domain_size() -
1)
: -model_proto.variables(positive_target).domain(0);
// For each X_i, we encode y_i => X <= X_i. And at least one of the y_i
// is true. Note that the correct y_i will be chosen because of the first
// part in linearlization (X >= X_i).
// TODO(user): Only lower bound is needed, experiment.
LinearConstraintBuilder lc_exactly_one(model, IntegerValue(1),
IntegerValue(1));
for (const int var : ct.int_max().vars()) {
if (target == var) continue;
// y => X <= X_i.
// <=> max_term_value * y + X - X_i <= max_term_value.
// where max_tern_value is X_ub - X_i_lb.
IntegerVariable y = model->Add(NewIntegerVariable(0, 1));
const int positive_var = PositiveRef(var);
const int64 min_var_value =
RefIsPositive(var)
? model_proto.variables(positive_var).domain(0)
: -model_proto.variables(positive_var)
.domain(
model_proto.variables(positive_var).domain_size() -
1);
int64 max_term_value = max_target_value - min_var_value;
LinearConstraintBuilder lc(model, kMinIntegerValue,
IntegerValue(max_term_value));
lc.AddTerm(mapping->Integer(target), IntegerValue(1));
lc.AddTerm(mapping->Integer(var), IntegerValue(-1));
lc.AddTerm(y, IntegerValue(max_term_value));
relaxation->linear_constraints.push_back(lc.Build());
lc_exactly_one.AddTerm(y, IntegerValue(1));
}
relaxation->linear_constraints.push_back(lc_exactly_one.Build());
}
} else if (ct.constraint_case() == ConstraintProto::ConstraintCase::kIntMin) {
if (HasEnforcementLiteral(ct)) return;
const int target = ct.int_min().target();
for (const int var : ct.int_min().vars()) {
if (target == var) continue;
LinearConstraintBuilder lc(model, kMinIntegerValue, IntegerValue(0));
lc.AddTerm(mapping->Integer(target), IntegerValue(1));
lc.AddTerm(mapping->Integer(var), IntegerValue(-1));
relaxation->linear_constraints.push_back(lc.Build());
}
} else if (ct.constraint_case() ==
ConstraintProto::ConstraintCase::kIntProd) {
if (HasEnforcementLiteral(ct)) return;
const int target = ct.int_prod().target();
const int size = ct.int_prod().vars_size();
// We just linearize x = y^2 by x >= y which is far from ideal but at
// least pushes x when y moves away from zero. Note that if y is negative,
// we should probably also add x >= -y, but then this do not happen in
// our test set.
if (size == 2 && ct.int_prod().vars(0) == ct.int_prod().vars(1)) {
LinearConstraintBuilder lc(model, kMinIntegerValue, IntegerValue(0));
lc.AddTerm(mapping->Integer(ct.int_prod().vars(0)), IntegerValue(1));
lc.AddTerm(mapping->Integer(target), IntegerValue(-1));
relaxation->linear_constraints.push_back(lc.Build());
}
} else if (ct.constraint_case() == ConstraintProto::ConstraintCase::kLinear) {
AppendLinearConstraintRelaxation(ct, linearization_level, *model,
relaxation);
} else if (ct.constraint_case() ==
ConstraintProto::ConstraintCase::kCircuit) {
if (HasEnforcementLiteral(ct)) return;
const int num_arcs = ct.circuit().literals_size();
CHECK_EQ(num_arcs, ct.circuit().tails_size());
CHECK_EQ(num_arcs, ct.circuit().heads_size());
// Each node must have exactly one incoming and one outgoing arc (note that
// it can be the unique self-arc of this node too).
std::map<int, std::vector<Literal>> incoming_arc_constraints;
std::map<int, std::vector<Literal>> outgoing_arc_constraints;
for (int i = 0; i < num_arcs; i++) {
const Literal arc = mapping->Literal(ct.circuit().literals(i));
const int tail = ct.circuit().tails(i);
const int head = ct.circuit().heads(i);
// Make sure this literal has a view.
model->Add(NewIntegerVariableFromLiteral(arc));
outgoing_arc_constraints[tail].push_back(arc);
incoming_arc_constraints[head].push_back(arc);
}
for (const auto* node_map :
{&outgoing_arc_constraints, &incoming_arc_constraints}) {
for (const auto& entry : *node_map) {
const std::vector<Literal>& exactly_one = entry.second;
if (exactly_one.size() > 1) {
LinearConstraintBuilder at_least_one_lc(model, IntegerValue(1),
kMaxIntegerValue);
for (const Literal l : exactly_one) {
CHECK(at_least_one_lc.AddLiteralTerm(l, IntegerValue(1)));
}
// We separate the two constraints.
relaxation->at_most_ones.push_back(exactly_one);
relaxation->linear_constraints.push_back(at_least_one_lc.Build());
}
}
}
} else if (ct.constraint_case() ==
ConstraintProto::ConstraintCase::kElement) {
const IntegerVariable index = mapping->Integer(ct.element().index());
const IntegerVariable target = mapping->Integer(ct.element().target());
const std::vector<IntegerVariable> vars =
mapping->Integers(ct.element().vars());
// We only relax the case where all the vars are constant.
// target = sum (index == i) * fixed_vars[i].
LinearConstraintBuilder constraint(model, IntegerValue(0), IntegerValue(0));
constraint.AddTerm(target, IntegerValue(-1));
IntegerTrail* integer_trail = model->GetOrCreate<IntegerTrail>();
for (const auto literal_value : model->Add(FullyEncodeVariable((index)))) {
const IntegerVariable var = vars[literal_value.value.value()];
if (!model->Get(IsFixed(var))) return;
// Make sure this literal has a view.
model->Add(NewIntegerVariableFromLiteral(literal_value.literal));
CHECK(constraint.AddLiteralTerm(literal_value.literal,
integer_trail->LowerBound(var)));
}
relaxation->linear_constraints.push_back(constraint.Build());
}
}
void AppendLinearConstraintRelaxation(const ConstraintProto& constraint_proto,
const int linearization_level,
const Model& model,

View File

@@ -69,6 +69,12 @@ void AppendPartialGreaterThanEncodingRelaxation(IntegerVariable var,
const Model& model,
LinearRelaxation* relaxation);
// Adds linearization of different types of constraints.
void TryToLinearizeConstraint(const CpModelProto& model_proto,
const ConstraintProto& ct, Model* model,
int linearization_level,
LinearRelaxation* relaxation);
// Appends linear constraints to the relaxation. This also handles the
// relaxation of linear constraints with enforcement literals.
// A linear constraint lb <= ax <= ub with enforcement literals {ei} is relaxed