From c36d1d7dc164a6ccdbe3e101bca185da2004e10c Mon Sep 17 00:00:00 2001 From: Laurent Perron Date: Tue, 6 Aug 2019 15:56:27 -0700 Subject: [PATCH] more type of constraints in the linear solver proto, support in SCIP --- ortools/linear_solver/linear_solver.proto | 29 ++ ortools/linear_solver/model_validator.cc | 79 +++++- ortools/linear_solver/scip_proto_solver.cc | 303 ++++++++++++++++++++- 3 files changed, 398 insertions(+), 13 deletions(-) diff --git a/ortools/linear_solver/linear_solver.proto b/ortools/linear_solver/linear_solver.proto index 1c5914ee1e..8a93c37e26 100644 --- a/ortools/linear_solver/linear_solver.proto +++ b/ortools/linear_solver/linear_solver.proto @@ -110,6 +110,9 @@ message MPGeneralConstraintProto { MPIndicatorConstraint indicator_constraint = 2; MPSosConstraint sos_constraint = 3; MPQuadraticConstraint quadratic_constraint = 4; + MPAbsConstraint abs_constraint = 5; + MPAndConstraint and_constraint = 6; + MPOrConstraint or_constraint = 7; } } @@ -194,6 +197,32 @@ message MPQuadraticConstraint { optional double upper_bound = 7 [default = inf]; } +// Sets a variable's value to the absolute value of another variable. +message MPAbsConstraint { + // Variable indices are relative to the "variable" field in MPModelProto. + // resultant_var = abs(var) + optional int32 var_index = 1; + optional int32 resultant_var_index = 2; +} + +// Sets a binary variable's value equal to one if and only if all variables in a +// set of binary variables are true. +message MPAndConstraint { + // Variable indices are relative to the "variable" field in MPModelProto. + // resultant_var = and(var_1, var_2... var_n) + repeated int32 var_index = 1; + optional int32 resultant_var_index = 2; +} + +// Sets a binary variable's value equal to one if and only if at least one +// variable in a set of binary variables is true. +message MPOrConstraint { + // Variable indices are relative to the "variable" field in MPModelProto. + // resultant_var = or(var_1, var_2... var_n) + repeated int32 var_index = 1; + optional int32 resultant_var_index = 2; +} + // Quadratic part of a model's objective. Added with other objectives (such as // linear), this creates the model's objective function to be optimized. // Note: the linear part of the objective currently needs to be specified in the diff --git a/ortools/linear_solver/model_validator.cc b/ortools/linear_solver/model_validator.cc index 05aae0cf29..4fe097250e 100644 --- a/ortools/linear_solver/model_validator.cc +++ b/ortools/linear_solver/model_validator.cc @@ -137,6 +137,12 @@ std::string CroppedConstraintDebugString(const MPConstraintProto& constraint) { ProtobufShortDebugString(constraint_light), suffix_str); } +bool IsBoolean(const MPVariableProto& variable) { + if (variable.lower_bound() < 0) return false; + if (variable.upper_bound() > 1) return false; + return variable.is_integer(); +} + std::string FindErrorInMPIndicatorConstraint( const MPModelProto& model, const MPIndicatorConstraint& indicator, std::vector* var_mask) { @@ -147,9 +153,7 @@ std::string FindErrorInMPIndicatorConstraint( if (var_index < 0 || var_index >= model.variable_size()) { return absl::StrCat("var_index=", var_index, " is out of bounds."); } - if (!model.variable(var_index).is_integer() || - model.variable(var_index).lower_bound() < 0 || - model.variable(var_index).upper_bound() > 1) { + if (!IsBoolean(model.variable(var_index))) { return absl::StrCat("var_index=", var_index, " is not Boolean."); } const int var_value = indicator.var_value(); @@ -203,7 +207,7 @@ std::string FindErrorInMPQuadraticConstraint(const MPModelProto& model, return "var_index_size() != coefficient_size()"; } for (int i = 0; i < qcst.var_index_size(); ++i) { - if (qcst.var_index(i) < 0 || qcst.var_index(i) >= model.variable_size()) { + if (qcst.var_index(i) < 0 || qcst.var_index(i) >= num_vars) { return absl::StrCat("var_index(", i, ")=", qcst.var_index(i), " is invalid.", " It must be in [0, ", num_vars, ")"); } @@ -238,6 +242,58 @@ std::string FindErrorInMPQuadraticConstraint(const MPModelProto& model, return ""; } +std::string FindErrorInMPAbsConstraint(const MPModelProto& model, + const MPAbsConstraint& abs) { + if (!abs.has_var_index()) { + return "var_index is required."; + } + if (!abs.has_resultant_var_index()) { + return "resultant_var_index is required."; + } + + const int num_vars = model.variable_size(); + if (abs.var_index() < 0 || abs.var_index() >= num_vars) { + return absl::StrCat("var_index=", abs.var_index(), " is invalid.", + " It must be in [0, ", num_vars, ")"); + } + if (abs.resultant_var_index() < 0 || abs.resultant_var_index() >= num_vars) { + return absl::StrCat("var_index=", abs.resultant_var_index(), " is invalid.", + " It must be in [0, ", num_vars, ")"); + } + return ""; +} + +template +std::string FindErrorInMPAndOrConstraint(const MPModelProto& model, + const MPAndOrConstraint& and_or) { + if (and_or.var_index_size() == 0) { + return "var_index cannot be empty."; + } + if (!and_or.has_resultant_var_index()) { + return "resultant_var_index is required."; + } + + const int num_vars = model.variable_size(); + for (int i = 0; i < and_or.var_index_size(); ++i) { + if (and_or.var_index(i) < 0 || and_or.var_index(i) >= num_vars) { + return absl::StrCat("var_index(", i, ")=", and_or.var_index(i), + " is invalid.", " It must be in [0, ", num_vars, ")"); + } + if (!IsBoolean(model.variable(and_or.var_index(i)))) { + return absl::StrCat("var_index=", i, " is not Boolean."); + } + } + if (and_or.resultant_var_index() < 0 || + and_or.resultant_var_index() >= num_vars) { + return absl::StrCat("resultant_var_index=", and_or.resultant_var_index(), + " is invalid.", " It must be in [0, ", num_vars, ")"); + } + if (!IsBoolean(model.variable(and_or.resultant_var_index()))) { + return absl::StrCat("resultant_var_index is not Boolean."); + } + return ""; +} + std::string FindErrorInQuadraticObjective(const MPQuadraticObjective& qobj, int num_vars) { if (qobj.qvar1_index_size() != qobj.qvar2_index_size() || @@ -346,6 +402,21 @@ std::string FindErrorInMPModelProto(const MPModelProto& model) { model, gen_constraint.quadratic_constraint(), &variable_appears); break; + case MPGeneralConstraintProto::kAbsConstraint: + error = + FindErrorInMPAbsConstraint(model, gen_constraint.abs_constraint()); + break; + + case MPGeneralConstraintProto::kAndConstraint: + error = FindErrorInMPAndOrConstraint(model, + gen_constraint.and_constraint()); + break; + + case MPGeneralConstraintProto::kOrConstraint: + error = + FindErrorInMPAndOrConstraint(model, gen_constraint.or_constraint()); + break; + default: return absl::StrCat("Unknown general constraint type ", gen_constraint.general_constraint_case()); diff --git a/ortools/linear_solver/scip_proto_solver.cc b/ortools/linear_solver/scip_proto_solver.cc index 5af76a4a84..121b0ed34c 100644 --- a/ortools/linear_solver/scip_proto_solver.cc +++ b/ortools/linear_solver/scip_proto_solver.cc @@ -34,13 +34,19 @@ #include "ortools/linear_solver/linear_solver.pb.h" #include "ortools/linear_solver/model_validator.h" #include "ortools/linear_solver/scip_helper_macros.h" +#include "scip/cons_disjunction.h" +#include "scip/cons_linear.h" +#include "scip/pub_var.h" #include "scip/scip.h" #include "scip/scip_param.h" +#include "scip/scip_prob.h" +#include "scip/scip_var.h" #include "scip/scipdefplugins.h" #include "scip/set.h" #include "scip/struct_paramset.h" #include "scip/type_cons.h" #include "scip/type_paramset.h" +#include "scip/type_var.h" namespace operations_research { @@ -128,6 +134,81 @@ util::Status ScipSetSolverSpecificParameters(const std::string& parameters, } namespace { +// This function will create a new constraint if the indicator constraint has +// both a lower bound and an upper bound. +util::Status AddIndicatorConstraint( + const MPGeneralConstraintProto& gen_cst, + const std::vector& scip_variables, SCIP* scip, + SCIP_CONS** scip_cst, std::vector* scip_constraints, + std::vector* tmp_variables, + std::vector* tmp_coefficients) { + CHECK(scip != nullptr); + CHECK(scip_cst != nullptr); + CHECK(scip_constraints != nullptr); + CHECK(tmp_variables != nullptr); + CHECK(tmp_coefficients != nullptr); + CHECK(gen_cst.has_indicator_constraint()); + constexpr double kInfinity = std::numeric_limits::infinity(); + + const auto& ind = gen_cst.indicator_constraint(); + if (!ind.has_constraint()) return util::OkStatus(); + + const MPConstraintProto& constraint = ind.constraint(); + const int size = constraint.var_index_size(); + tmp_variables->resize(size, nullptr); + tmp_coefficients->resize(size, 0); + for (int i = 0; i < size; ++i) { + (*tmp_variables)[i] = scip_variables[constraint.var_index(i)]; + (*tmp_coefficients)[i] = constraint.coefficient(i); + } + + SCIP_VAR* ind_var = scip_variables[ind.var_index()]; + if (ind.var_value() == 0) { + RETURN_IF_SCIP_ERROR( + SCIPgetNegatedVar(scip, scip_variables[ind.var_index()], &ind_var)); + } + + if (ind.constraint().upper_bound() < kInfinity) { + RETURN_IF_SCIP_ERROR(SCIPcreateConsIndicator( + scip, scip_cst, gen_cst.name().c_str(), ind_var, size, + tmp_variables->data(), tmp_coefficients->data(), + ind.constraint().upper_bound(), + /*initial=*/!ind.constraint().is_lazy(), + /*separate=*/true, + /*enforce=*/true, + /*check=*/true, + /*propagate=*/true, + /*local=*/false, + /*dynamic=*/false, + /*removable=*/ind.constraint().is_lazy(), + /*stickingatnode=*/false)); + RETURN_IF_SCIP_ERROR(SCIPaddCons(scip, *scip_cst)); + scip_constraints->push_back(nullptr); + scip_cst = &scip_constraints->back(); + } + if (ind.constraint().lower_bound() > -kInfinity) { + for (int i = 0; i < size; ++i) { + (*tmp_coefficients)[i] *= -1; + } + RETURN_IF_SCIP_ERROR(SCIPcreateConsIndicator( + scip, scip_cst, gen_cst.name().c_str(), ind_var, size, + tmp_variables->data(), tmp_coefficients->data(), + -ind.constraint().lower_bound(), + /*initial=*/!ind.constraint().is_lazy(), + /*separate=*/true, + /*enforce=*/true, + /*check=*/true, + /*propagate=*/true, + /*local=*/false, + /*dynamic=*/false, + /*removable=*/ind.constraint().is_lazy(), + /*stickingatnode=*/false)); + RETURN_IF_SCIP_ERROR(SCIPaddCons(scip, *scip_cst)); + } + + return util::OkStatus(); +} + util::Status AddSosConstraint(const MPGeneralConstraintProto& gen_cst, const std::vector& scip_variables, SCIP* scip, SCIP_CONS** scip_cst, @@ -247,6 +328,111 @@ util::Status AddQuadraticConstraint( return util::OkStatus(); } +// Models the constraint y = |x| as y >= 0 plus one disjunction constraint: +// y = x OR y = -x +util::Status AddAbsConstraint(const MPGeneralConstraintProto& gen_cst, + const std::vector& scip_variables, + SCIP* scip, + std::vector* scip_constraints) { + CHECK(scip != nullptr); + CHECK(scip_constraints != nullptr); + CHECK(gen_cst.has_abs_constraint()); + const auto& abs = gen_cst.abs_constraint(); + SCIP_VAR* scip_var = scip_variables[abs.var_index()]; + SCIP_VAR* scip_resultant_var = scip_variables[abs.resultant_var_index()]; + + // Set the resultant variable's lower bound to zero if it's negative. + if (SCIPvarGetLbLocal(scip_resultant_var) < 0.0) { + RETURN_IF_SCIP_ERROR(SCIPchgVarLb(scip, scip_resultant_var, 0.0)); + } + + std::vector vars; + std::vector vals; + std::vector cons; + auto add_abs_constraint = + [&](const std::string& name_prefix) -> util::Status { + scip_constraints->push_back(nullptr); + const std::string name = + gen_cst.has_name() ? absl::StrCat(gen_cst.name(), name_prefix) : ""; + RETURN_IF_SCIP_ERROR(SCIPcreateConsBasicLinear( + scip, /*cons=*/&scip_constraints->back(), + /*name=*/name.c_str(), /*nvars=*/2, /*vars=*/vars.data(), + /*vals=*/vals.data(), /*lhs=*/0.0, /*rhs=*/0.0)); + // Note that the constraints are, by design, not added into the model using + // SCIPaddCons. + return util::OkStatus(); + }; + + // Create an intermediary constraint such that y = -x + vars = {scip_resultant_var, scip_var}; + vals = {1, 1}; + RETURN_IF_ERROR(add_abs_constraint("_neg")); + cons.push_back(scip_constraints->back()); + + // Create an intermediary constraint such that y = x + vals = {1, -1}; + RETURN_IF_ERROR(add_abs_constraint("_pos")); + cons.push_back(scip_constraints->back()); + + // Activate at least one of the two above constraints. + const std::string name = + gen_cst.has_name() ? absl::StrCat(gen_cst.name(), "_disj") : ""; + RETURN_IF_SCIP_ERROR(SCIPcreateConsBasicDisjunction( + scip, /*cons=*/&scip_constraints->back(), /*name=*/name.c_str(), + /*nconss=*/2, /*conss=*/cons.data(), /*relaxcons=*/nullptr)); + RETURN_IF_SCIP_ERROR(SCIPaddCons(scip, scip_constraints->back())); + + return util::OkStatus(); +} + +util::Status AddAndConstraint(const MPGeneralConstraintProto& gen_cst, + const std::vector& scip_variables, + SCIP* scip, SCIP_CONS** scip_cst, + std::vector* tmp_variables) { + CHECK(scip != nullptr); + CHECK(scip_cst != nullptr); + CHECK(tmp_variables != nullptr); + CHECK(gen_cst.has_and_constraint()); + const auto& andcst = gen_cst.and_constraint(); + + tmp_variables->resize(andcst.var_index_size(), nullptr); + for (int i = 0; i < andcst.var_index_size(); ++i) { + (*tmp_variables)[i] = scip_variables[andcst.var_index(i)]; + } + RETURN_IF_SCIP_ERROR(SCIPcreateConsBasicAnd( + scip, /*cons=*/scip_cst, + /*name=*/gen_cst.name().c_str(), + /*resvar=*/scip_variables[andcst.resultant_var_index()], + /*nvars=*/andcst.var_index_size(), + /*vars=*/tmp_variables->data())); + RETURN_IF_SCIP_ERROR(SCIPaddCons(scip, *scip_cst)); + return util::OkStatus(); +} + +util::Status AddOrConstraint(const MPGeneralConstraintProto& gen_cst, + const std::vector& scip_variables, + SCIP* scip, SCIP_CONS** scip_cst, + std::vector* tmp_variables) { + CHECK(scip != nullptr); + CHECK(scip_cst != nullptr); + CHECK(tmp_variables != nullptr); + CHECK(gen_cst.has_or_constraint()); + const auto& orcst = gen_cst.or_constraint(); + + tmp_variables->resize(orcst.var_index_size(), nullptr); + for (int i = 0; i < orcst.var_index_size(); ++i) { + (*tmp_variables)[i] = scip_variables[orcst.var_index(i)]; + } + RETURN_IF_SCIP_ERROR(SCIPcreateConsBasicOr( + scip, /*cons=*/scip_cst, + /*name=*/gen_cst.name().c_str(), + /*resvar=*/scip_variables[orcst.resultant_var_index()], + /*nvars=*/orcst.var_index_size(), + /*vars=*/tmp_variables->data())); + RETURN_IF_SCIP_ERROR(SCIPaddCons(scip, *scip_cst)); + return util::OkStatus(); +} + util::Status AddQuadraticObjective(const MPQuadraticObjective& quadobj, SCIP* scip, std::vector* scip_variables, @@ -291,20 +477,96 @@ util::Status AddQuadraticObjective(const MPQuadraticObjective& quadobj, return util::OkStatus(); } + +util::Status AddSolutionHint(const MPModelProto& model, SCIP* scip, + const std::vector& scip_variables) { + CHECK(scip != nullptr); + if (!model.has_solution_hint()) return util::OkStatus(); + + const PartialVariableAssignment& solution_hint = model.solution_hint(); + SCIP_SOL* solution; + bool is_solution_partial = + solution_hint.var_index_size() != model.variable_size(); + if (is_solution_partial) { + RETURN_IF_SCIP_ERROR( + SCIPcreatePartialSol(scip, /*sol=*/&solution, /*heur=*/nullptr)); + } else { + RETURN_IF_SCIP_ERROR( + SCIPcreateSol(scip, /*sol=*/&solution, /*heur=*/nullptr)); + } + + for (int i = 0; i < solution_hint.var_index_size(); ++i) { + RETURN_IF_SCIP_ERROR(SCIPsetSolVal( + scip, solution, scip_variables[solution_hint.var_index(i)], + solution_hint.var_value(i))); + } + + SCIP_Bool is_stored; + RETURN_IF_SCIP_ERROR(SCIPaddSolFree(scip, &solution, &is_stored)); + + return util::OkStatus(); +} + +bool MPModelIsInvalidForScip(const MPModelProto& model, SCIP* scip, + MPSolutionResponse* response) { + CHECK(scip != nullptr); + CHECK(response != nullptr); + + const double infinity = SCIPinfinity(scip); + for (int v = 0; v < model.variable_size(); ++v) { + const MPVariableProto& variable = model.variable(v); + if (variable.lower_bound() >= infinity) { + response->set_status(MPSOLVER_MODEL_INVALID); + response->set_status_str(absl::StrFormat( + "Variable %i's lower bound is considered +infinity", v)); + return true; + } + if (variable.upper_bound() <= -infinity) { + response->set_status(MPSOLVER_MODEL_INVALID); + response->set_status_str(absl::StrFormat( + "Variable %i's lower bound is considered -infinity", v)); + return true; + } + const double coeff = variable.objective_coefficient(); + if (coeff >= infinity || coeff <= -infinity) { + response->set_status(MPSOLVER_MODEL_INVALID); + response->set_status_str(absl::StrFormat( + "Variable %i's objective coefficient is considered infinite", v)); + return true; + } + } + + if (model.has_solution_hint()) { + for (int i = 0; i < model.solution_hint().var_value_size(); ++i) { + const double value = model.solution_hint().var_value(i); + if (value >= infinity || value <= -infinity) { + response->set_status(MPSOLVER_MODEL_INVALID); + response->set_status_str(absl::StrFormat( + "Variable %i's solution hint is considered infinite", + model.solution_hint().var_index(i))); + return true; + } + } + } + + if (model.objective_offset() >= infinity || + model.objective_offset() <= -infinity) { + response->set_status(MPSOLVER_MODEL_INVALID); + response->set_status_str( + "Model's objective offset is considered infinite."); + return true; + } + + return false; +} } // namespace util::StatusOr ScipSolveProto( const MPModelRequest& request) { MPSolutionResponse response; - if (MPRequestIsEmptyOrInvalid(request, &response)) { - return response; - } + if (MPRequestIsEmptyOrInvalid(request, &response)) return response; const MPModelProto& model = request.model(); - if (model.has_solution_hint()) { - // TODO(user): Support solution hints. - return util::UnimplementedError("Solution hint not supported."); - } SCIP* scip = nullptr; std::vector scip_variables(model.variable_size(), nullptr); @@ -335,6 +597,7 @@ util::StatusOr ScipSolveProto( RETURN_IF_SCIP_ERROR(SCIPcreate(&scip)); RETURN_IF_SCIP_ERROR(SCIPincludeDefaultPlugins(scip)); + if (MPModelIsInvalidForScip(model, scip, &response)) return response; const auto parameters_status = ScipSetSolverSpecificParameters( request.solver_specific_parameters(), scip); @@ -404,9 +667,13 @@ util::StatusOr ScipSolveProto( const int lincst_size = model.constraint_size(); for (int c = 0; c < model.general_constraint_size(); ++c) { const MPGeneralConstraintProto& gen_cst = model.general_constraint(c); - // TODO(user): Move indicator constraint logic from linear_solver.cc - // to this file. switch (gen_cst.general_constraint_case()) { + case MPGeneralConstraintProto::kIndicatorConstraint: { + RETURN_IF_ERROR(AddIndicatorConstraint( + gen_cst, scip_variables, scip, &scip_constraints[lincst_size + c], + &scip_constraints, &ct_variables, &ct_coefficients)); + break; + } case MPGeneralConstraintProto::kSosConstraint: { RETURN_IF_ERROR(AddSosConstraint(gen_cst, scip_variables, scip, &scip_constraints[lincst_size + c], @@ -420,6 +687,23 @@ util::StatusOr ScipSolveProto( &ct_qcoefficients)); break; } + case MPGeneralConstraintProto::kAbsConstraint: { + RETURN_IF_ERROR(AddAbsConstraint(gen_cst, scip_variables, scip, + &scip_constraints)); + break; + } + case MPGeneralConstraintProto::kAndConstraint: { + RETURN_IF_ERROR(AddAndConstraint(gen_cst, scip_variables, scip, + &scip_constraints[lincst_size + c], + &ct_variables)); + break; + } + case MPGeneralConstraintProto::kOrConstraint: { + RETURN_IF_ERROR(AddOrConstraint(gen_cst, scip_variables, scip, + &scip_constraints[lincst_size + c], + &ct_variables)); + break; + } default: return util::UnimplementedError( absl::StrFormat("General constraints of type %i not supported.", @@ -433,6 +717,7 @@ util::StatusOr ScipSolveProto( &scip_variables, &scip_constraints)); } RETURN_IF_SCIP_ERROR(SCIPaddOrigObjoffset(scip, model.objective_offset())); + RETURN_IF_ERROR(AddSolutionHint(model, scip, scip_variables)); RETURN_IF_SCIP_ERROR(SCIPsolve(scip));