add branching priorities to MPSolver (supported by SCIP and Gurobi)

This commit is contained in:
Laurent Perron
2019-06-17 11:35:17 +02:00
parent f5396c9d1f
commit 7487f03c65
6 changed files with 215 additions and 76 deletions

View File

@@ -84,6 +84,7 @@ class GurobiInterface : public MPSolverInterface {
// Clears the objective from all its terms.
void ClearObjective() override;
bool CheckBestObjectiveBoundExists() const override;
void BranchingPriorityChangedForVariable(int var_index) override;
// ------ Query statistics on the solution and the solve ------
// Number of simplex or interior-point iterations
@@ -175,6 +176,7 @@ class GurobiInterface : public MPSolverInterface {
GRBenv* env_;
bool mip_;
int current_solution_index_;
bool update_branching_priorities_ = false;
};
namespace {
@@ -320,6 +322,10 @@ void GurobiInterface::SetObjectiveOffset(double value) {
void GurobiInterface::ClearObjective() { sync_status_ = MUST_RELOAD; }
void GurobiInterface::BranchingPriorityChangedForVariable(int var_index) {
update_branching_priorities_ = true;
}
// ------ Query statistics on the solution and the solve ------
int64 GurobiInterface::iterations() const {
@@ -705,6 +711,16 @@ MPSolver::ResultStatus GurobiInterface::Solve(const MPSolverParameters& param) {
GRBsetdblattrelement(model_, "Start", p.first->index(), p.second));
}
// Pass branching priority annotations if at least one has been updated.
if (update_branching_priorities_) {
for (const MPVariable* var : solver_->variables_) {
CheckedGurobiCall(
GRBsetintattrelement(model_, GRB_INT_ATTR_BRANCHPRIORITY,
var->index(), var->branching_priority()));
}
update_branching_priorities_ = false;
}
// Time limit.
if (solver_->time_limit() != 0) {
VLOG(1) << "Setting time limit = " << solver_->time_limit() << " ms.";

View File

@@ -285,6 +285,12 @@ void MPVariable::SetInteger(bool integer) {
}
}
void MPVariable::SetBranchingPriority(int priority) {
if (priority == branching_priority_) return;
branching_priority_ = priority;
interface_->BranchingPriorityChangedForVariable(index_);
}
// ----- Interface shortcuts -----
bool MPSolver::IsMIP() const { return interface_->IsMIP(); }
@@ -601,6 +607,9 @@ MPSolverResponseStatus MPSolver::LoadModelFromProtoInternal(
MakeNumVar(var_proto.lower_bound(), var_proto.upper_bound(),
clear_names ? empty : var_proto.name());
variable->SetInteger(var_proto.is_integer());
if (var_proto.branching_priority() != 0) {
variable->SetBranchingPriority(var_proto.branching_priority());
}
objective->SetCoefficient(variable, var_proto.objective_coefficient());
}
@@ -774,6 +783,9 @@ void MPSolver::ExportModelToProto(MPModelProto* output_model) const {
variable_proto->set_objective_coefficient(
objective_->GetCoefficient(var));
}
if (var->branching_priority() != 0) {
variable_proto->set_branching_priority(var->branching_priority());
}
}
// Map the variables to their indices. This is needed to output the

View File

@@ -393,7 +393,7 @@ class MPSolver {
// Loads model from protocol buffer. Returns MPSOLVER_MODEL_IS_VALID if the
// model is valid, and another status otherwise (currently only
// MPSOLVER_MODEL_INVALID and MPSOLVER_INFEASIBLE). If the model isn't
// valid, populate "error_message".
// valid, populates "error_message".
MPSolverResponseStatus LoadModelFromProto(const MPModelProto& input_model,
std::string* error_message);
// The same as above, except that the loading keeps original variable and
@@ -885,6 +885,16 @@ class MPVariable {
// @see MPSolver::BasisStatus.
MPSolver::BasisStatus basis_status() const;
// Advanced usage: Certain MIP solvers (e.g. Gurobi or SCIP) allow you to set
// a per-variable priority for determining which variable to branch on. A
// value of 0 is treated as default, and is equivalent to not setting the
// branching priority. The solver looks first to branch on fractional
// variables in higher priority levels. As of 2019-05, only Gurobi and SCIP
// support setting branching priority; all other solvers will simply ignore
// this annotation.
int branching_priority() const { return branching_priority_; }
void SetBranchingPriority(int priority);
protected:
friend class MPSolver;
friend class MPSolverInterface;
@@ -926,6 +936,7 @@ class MPVariable {
const std::string name_;
double solution_value_;
double reduced_cost_;
int branching_priority_ = 0;
MPSolverInterface* const interface_;
DISALLOW_COPY_AND_ASSIGN(MPVariable);
};
@@ -1303,6 +1314,7 @@ class MPSolverInterface {
// Clears the objective from all its terms.
virtual void ClearObjective() = 0;
virtual void BranchingPriorityChangedForVariable(int var_index) {}
// ------ Query statistics on the solution and the solve ------
// Returns the number of simplex iterations. The problem must be discrete,
// otherwise it crashes, or returns kUnknownNumberOfIterations in NDEBUG mode.

View File

@@ -64,6 +64,8 @@ message MPVariableProto {
// The name of the variable.
optional string name = 5 [default = ""];
optional int32 branching_priority = 6 [default = 0];
}
// A linear constraint is always of the form:
@@ -106,6 +108,7 @@ message MPGeneralConstraintProto {
oneof general_constraint {
MPIndicatorConstraint indicator_constraint = 2;
MPSosConstraint sos_constraint = 3;
}
}
@@ -115,6 +118,7 @@ message MPGeneralConstraintProto {
// y = 0 => 2 * x1 + 3 * x2 >= 42
// The 2 * x1 + 3 * x2 >= 42 constraint is only active if the variable y is
// equal to 0.
// As of 2019/04, only SCIP, CP-SAT and Gurobi support this constraint type.
message MPIndicatorConstraint {
// Variable index (w.r.t. the "variable" field of MPModelProto) of the Boolean
// variable used as indicator.
@@ -127,6 +131,36 @@ message MPIndicatorConstraint {
optional MPConstraintProto constraint = 3;
}
// Special Ordered Set (SOS) constraints of type 1 or 2.
// See https://en.wikipedia.org/wiki/Special_ordered_set
// As of 2019/04, only SCIP and Gurobi support this constraint type.
message MPSosConstraint {
enum Type {
// At most one variable in `var_index` must be non-zero.
SOS1_DEFAULT = 0;
// At most two consecutive variables from `var_index` must be non-zero (i.e.
// for some i, var_index[i] and var_index[i+1]). See
// http://www.eudoxus.com/lp-training/5/5-6-special-ordered-sets-of-type-2
SOS2 = 1;
}
optional Type type = 1 [default = SOS1_DEFAULT];
// Variable index (w.r.t. the "variable" field of MPModelProto) of the
// variables in the SOS.
repeated int32 var_index = 2;
// Optional: SOS weights. If non-empty, must be of the same size as
// "var_index", and strictly increasing.
// SUBTLE: The weights can help the solver make branch-and-bound decisions
// that fit the underlying optimization model: after each LP relaxation, it
// will compute the "average weight" of the SOS variables, weighted by value
// (this is confusing: here we're using the values as weights), and the binary
// branch decision will be: is the non-zero variable above or below that?
// (weights are strictly monotonous, so the "cutoff" average weight
// corresponds to a "cutoff" index in the var_index sequence).
repeated double weight = 3;
}
// This message encodes a partial (or full) assignment of the variables of a
// MPModelProto problem. The indices in var_index should be unique and valid
// variable indices of the associated problem.

View File

@@ -13,6 +13,7 @@
#include "ortools/linear_solver/model_validator.h"
#include <algorithm>
#include <cmath>
#include <limits>
@@ -50,6 +51,26 @@ std::string FindErrorInMPVariable(const MPVariableProto& variable) {
return std::string();
}
// Returns an error message if 'var_indices' contains a duplicate index.
template <typename Iterable>
std::string FindDuplicateVarIndex(const Iterable& var_indices,
std::vector<bool>* var_mask) {
int duplicate_var_index = -1;
for (const int var_index : var_indices) {
if ((*var_mask)[var_index]) duplicate_var_index = var_index;
(*var_mask)[var_index] = true;
}
// Reset "var_mask" to all false, sparsely.
for (const int var_index : var_indices) {
(*var_mask)[var_index] = false;
}
if (duplicate_var_index >= 0) {
return absl::StrCat("var_index #", duplicate_var_index,
" appears several times");
}
return "";
}
// Internal method to detect errors in a single constraint.
// "var_mask" is a std::vector<bool> whose size is the number of variables in
// the model, and it will be all set to false before and after the call.
@@ -86,25 +107,91 @@ std::string FindErrorInMPConstraint(const MPConstraintProto& constraint,
}
}
// Verify that the var_index don't have duplicates. We use "var_mask".
int duplicate_var_index = -1;
for (const int var_index : constraint.var_index()) {
if ((*var_mask)[var_index]) duplicate_var_index = var_index;
(*var_mask)[var_index] = true;
}
// Reset "var_mask" to all false, sparsely.
for (const int var_index : constraint.var_index()) {
(*var_mask)[var_index] = false;
}
if (duplicate_var_index >= 0) {
return absl::StrCat("var_index #", duplicate_var_index,
" appears several times");
}
const std::string error =
FindDuplicateVarIndex(constraint.var_index(), var_mask);
if (!error.empty()) return error;
// We found no error, all is fine.
return std::string();
}
std::string CroppedConstraintDebugString(const MPConstraintProto& constraint) {
const int kMaxPrintedVars = 10;
MPConstraintProto constraint_light = constraint;
std::string suffix_str;
if (constraint.var_index_size() > kMaxPrintedVars) {
constraint_light.mutable_var_index()->Truncate(kMaxPrintedVars);
absl::StrAppend(&suffix_str,
" (var_index cropped; size=", constraint.var_index_size(),
").");
}
if (constraint.coefficient_size() > kMaxPrintedVars) {
constraint_light.mutable_coefficient()->Truncate(kMaxPrintedVars);
absl::StrAppend(&suffix_str, " (coefficient cropped; size=",
constraint.coefficient_size(), ").");
}
return absl::StrCat("Constraint proto: ",
ProtobufShortDebugString(constraint_light), suffix_str);
}
std::string FindErrorInMPIndicatorConstraint(
const MPModelProto& model, const MPIndicatorConstraint& indicator,
std::vector<bool>* var_mask) {
if (!indicator.has_var_index()) {
return "var_index is required.";
}
const int var_index = indicator.var_index();
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) {
return absl::StrCat("var_index=", var_index, " is not Boolean.");
}
const int var_value = indicator.var_value();
if (var_value < 0 || var_value > 1) {
return absl::StrCat("var_value=", var_value, " must be 0 or 1.");
}
const MPConstraintProto& constraint = indicator.constraint();
std::string error = FindErrorInMPConstraint(constraint, var_mask);
if (!error.empty()) {
// Constraint protos can be huge, theoretically. So we guard against
// that.
return absl::StrCat(error, " in constraint ",
CroppedConstraintDebugString(constraint));
}
return "";
}
std::string FindErrorInMPSosConstraint(const MPModelProto& model,
const MPSosConstraint& sos,
std::vector<bool>* var_mask) {
if (sos.weight_size() != 0 && sos.weight_size() != sos.var_index_size()) {
return "weight_size() > 0 and var_index_size() != weight_size()";
}
for (const int var_index : sos.var_index()) {
if (var_index < 0 || var_index >= model.variable_size()) {
return absl::StrCat("var_index=", var_index, " is out of bounds.");
}
}
for (int i = 0; i < sos.weight_size(); ++i) {
if (!std::isfinite(sos.weight(i))) {
return absl::StrCat("Invalid weight: ", sos.weight(i));
}
if (i == 0) continue;
if (sos.weight(i - 1) >= sos.weight(i)) {
return "SOS weights must be strictly increasing";
}
}
const std::string error = FindDuplicateVarIndex(sos.var_index(), var_mask);
if (!error.empty()) return error;
return "";
}
std::string FindErrorInSolutionHint(
const PartialVariableAssignment& solution_hint, int num_vars) {
if (solution_hint.var_index_size() != solution_hint.var_value_size()) {
@@ -130,28 +217,6 @@ std::string FindErrorInSolutionHint(
}
return std::string();
}
std::string GetCroppedConstraintError(const MPConstraintProto& constraint,
int constraint_index,
const std::string& error,
int max_printed_vars) {
MPConstraintProto constraint_light = constraint;
std::string suffix_str;
if (constraint.var_index_size() > max_printed_vars) {
constraint_light.mutable_var_index()->Truncate(max_printed_vars);
absl::StrAppend(&suffix_str,
" (var_index cropped; size=", constraint.var_index_size(),
").");
}
if (constraint.coefficient_size() > max_printed_vars) {
constraint_light.mutable_coefficient()->Truncate(max_printed_vars);
absl::StrAppend(&suffix_str, " (coefficient cropped; size=",
constraint.coefficient_size(), ").");
}
return absl::StrCat("In constraint #", constraint_index, ": ", error,
". Constraint proto: ",
ProtobufShortDebugString(constraint_light), suffix_str);
}
} // namespace
std::string FindErrorInMPModelProto(const MPModelProto& model) {
@@ -180,14 +245,13 @@ std::string FindErrorInMPModelProto(const MPModelProto& model) {
// Validate constraints.
std::vector<bool> variable_appears(num_vars, false);
const int kMaxNumVarsInPrintedConstraint = 10;
for (int i = 0; i < num_cts; ++i) {
const MPConstraintProto& constraint = model.constraint(i);
error = FindErrorInMPConstraint(constraint, &variable_appears);
if (!error.empty()) {
// Constraint protos can be huge, theoretically. So we guard against that.
return GetCroppedConstraintError(constraint, i, error,
kMaxNumVarsInPrintedConstraint);
return absl::StrCat("In constraint #", i, ": ", error, ". ",
CroppedConstraintDebugString(constraint));
}
}
@@ -195,45 +259,24 @@ std::string FindErrorInMPModelProto(const MPModelProto& model) {
for (int i = 0; i < model.general_constraint_size(); ++i) {
const MPGeneralConstraintProto& gen_constraint =
model.general_constraint(i);
std::string error;
switch (gen_constraint.general_constraint_case()) {
case MPGeneralConstraintProto::kIndicatorConstraint: {
if (!gen_constraint.indicator_constraint().has_var_index()) {
return absl::StrCat("In general constraint #", i,
": var_index is required.");
}
const int var_index = gen_constraint.indicator_constraint().var_index();
if (var_index < 0 || var_index >= num_vars) {
return absl::StrCat("In general constraint #", i,
": 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) {
return absl::StrCat("In general constraint #", i,
": var_index=", var_index, " is not Boolean.");
}
const int var_value = gen_constraint.indicator_constraint().var_value();
if (var_value < 0 || var_value > 1) {
return absl::StrCat("In general constraint #", i,
": var_value=", var_value, " is invalid.");
}
const MPConstraintProto& constraint =
gen_constraint.indicator_constraint().constraint();
error = FindErrorInMPConstraint(constraint, &variable_appears);
if (!error.empty()) {
// Constraint protos can be huge, theoretically. So we guard against
// that.
return absl::StrCat(
"In general constraint #", i, ": ",
GetCroppedConstraintError(constraint, i, error,
kMaxNumVarsInPrintedConstraint));
}
case MPGeneralConstraintProto::kIndicatorConstraint:
error = FindErrorInMPIndicatorConstraint(
model, gen_constraint.indicator_constraint(), &variable_appears);
break;
}
default: {
case MPGeneralConstraintProto::kSosConstraint:
error = FindErrorInMPSosConstraint(
model, gen_constraint.sos_constraint(), &variable_appears);
break;
default:
return absl::StrCat("Unknown general constraint type ",
gen_constraint.general_constraint_case());
}
}
if (!error.empty()) {
return absl::StrCat("In general constraint #", i, ": ", error);
}
}

View File

@@ -62,6 +62,7 @@ class SCIPInterface : public MPSolverInterface {
double coefficient) override;
void SetObjectiveOffset(double value) override;
void ClearObjective() override;
void BranchingPriorityChangedForVariable(int var_index) override;
int64 iterations() const override;
int64 nodes() const override;
@@ -138,6 +139,7 @@ class SCIPInterface : public MPSolverInterface {
SCIP_VAR* objective_offset_variable_;
std::vector<SCIP_VAR*> scip_variables_;
std::vector<SCIP_CONS*> scip_constraints_;
bool branching_priority_reset_ = false;
};
// Our own version of SCIP_CALL to do error management.
@@ -368,6 +370,18 @@ void SCIPInterface::ClearObjective() {
RETURN_IF_SCIP_ERROR(SCIPchgVarObj(scip_, objective_offset_variable_, 0.0));
}
void SCIPInterface::BranchingPriorityChangedForVariable(int var_index) {
// As of 2019-05, SCIP does not support setting branching priority for
// variables in models that have already been solved. Therefore, we force
// reset the model when setting the priority on an already extracted variable.
// Note that this is a more drastic step than merely changing the sync_status.
// This may be slightly conservative, as it is technically possible that
// the extraction has occurred without a call to Solve().
if (variable_is_extracted(var_index)) {
branching_priority_reset_ = true;
}
}
void SCIPInterface::AddRowConstraint(MPConstraint* ct) {
sync_status_ = MUST_RELOAD;
}
@@ -399,6 +413,12 @@ void SCIPInterface::ExtractNewVariables() {
false, nullptr, nullptr, nullptr, nullptr, nullptr));
RETURN_IF_SCIP_ERROR(SCIPaddVar(scip_, scip_var));
scip_variables_.push_back(scip_var);
const int branching_priority = var->branching_priority();
if (branching_priority != 0) {
const int index = var->index();
RETURN_IF_SCIP_ERROR(SCIPchgVarBranchPriority(
scip_, scip_variables_[index], branching_priority));
}
}
// Add new variables to existing constraints.
for (int i = 0; i < last_constraint_index_; i++) {
@@ -564,8 +584,10 @@ MPSolver::ResultStatus SCIPInterface::Solve(const MPSolverParameters& param) {
// Note that SCIP does not provide any incrementality.
// TODO(user): Is that still true now (2018) ?
if (param.GetIntegerParam(MPSolverParameters::INCREMENTALITY) ==
MPSolverParameters::INCREMENTALITY_OFF) {
MPSolverParameters::INCREMENTALITY_OFF ||
branching_priority_reset_) {
Reset();
branching_priority_reset_ = false;
}
// Set log level.