more type of constraints in the linear solver proto, support in SCIP

This commit is contained in:
Laurent Perron
2019-08-06 15:56:27 -07:00
parent 747ee3d562
commit c36d1d7dc1
3 changed files with 398 additions and 13 deletions

View File

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

View File

@@ -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<bool>* 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 <typename MPAndOrConstraint>
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());

View File

@@ -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_VAR*>& scip_variables, SCIP* scip,
SCIP_CONS** scip_cst, std::vector<SCIP_CONS*>* scip_constraints,
std::vector<SCIP_VAR*>* tmp_variables,
std::vector<double>* 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<double>::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_VAR*>& 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_VAR*>& scip_variables,
SCIP* scip,
std::vector<SCIP_CONS*>* 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<SCIP_VAR*> vars;
std::vector<double> vals;
std::vector<SCIP_CONS*> 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_VAR*>& scip_variables,
SCIP* scip, SCIP_CONS** scip_cst,
std::vector<SCIP_VAR*>* 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_VAR*>& scip_variables,
SCIP* scip, SCIP_CONS** scip_cst,
std::vector<SCIP_VAR*>* 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_VAR*>* 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_VAR*>& 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<MPSolutionResponse> 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_VAR*> scip_variables(model.variable_size(), nullptr);
@@ -335,6 +597,7 @@ util::StatusOr<MPSolutionResponse> 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<MPSolutionResponse> 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<MPSolutionResponse> 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<MPSolutionResponse> 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));