diff --git a/ortools/constraint_solver/routing.cc b/ortools/constraint_solver/routing.cc index 153a6d9988..a863578cb3 100644 --- a/ortools/constraint_solver/routing.cc +++ b/ortools/constraint_solver/routing.cc @@ -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 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)); } diff --git a/ortools/constraint_solver/routing.h b/ortools/constraint_solver/routing.h index 141c5d002a..dc6dde9038 100644 --- a/ortools/constraint_solver/routing.h +++ b/ortools/constraint_solver/routing.h @@ -1544,6 +1544,7 @@ class RoutingModel { std::vector > > temporal_required_type_alternatives_per_type_index_; + absl::flat_hash_set trivially_infeasible_visit_types_; bool has_temporal_type_requirements_; // clang-format on int num_visit_types_; diff --git a/ortools/linear_solver/cbc_interface.cc b/ortools/linear_solver/cbc_interface.cc index 50c32fc24e..4eab6c21bd 100644 --- a/ortools/linear_solver/cbc_interface.cc +++ b/ortools/linear_solver/cbc_interface.cc @@ -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 diff --git a/ortools/linear_solver/gurobi_interface.cc b/ortools/linear_solver/gurobi_interface.cc index 2918a88e38..dcc993e9fe 100644 --- a/ortools/linear_solver/gurobi_interface.cc +++ b/ortools/linear_solver/gurobi_interface.cc @@ -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."); diff --git a/ortools/sat/cp_model_solver.cc b/ortools/sat/cp_model_solver.cc index 1de74c8799..dc4f0a0e00 100644 --- a/ortools/sat/cp_model_solver.cc +++ b/ortools/sat/cp_model_solver.cc @@ -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(); - 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 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> incoming_arc_constraints; - std::map> 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& 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 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(); - 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) { diff --git a/ortools/sat/linear_relaxation.cc b/ortools/sat/linear_relaxation.cc index 4f9c570848..4a33cab790 100644 --- a/ortools/sat/linear_relaxation.cc +++ b/ortools/sat/linear_relaxation.cc @@ -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(); + 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 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> incoming_arc_constraints; + std::map> 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& 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 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(); + 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, diff --git a/ortools/sat/linear_relaxation.h b/ortools/sat/linear_relaxation.h index 5a23e7f72c..c5a4c0b508 100644 --- a/ortools/sat/linear_relaxation.h +++ b/ortools/sat/linear_relaxation.h @@ -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