always-cautious/ortools-clone
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

Update math_opt

Browse Source
This commit is contained in:
Corentin Le Molgat
2022-03-02 22:10:23 +01:00
parent 1eba970bf3
commit 156190eea8
63 changed files with 4358 additions and 742 deletions
Download Patch File Download Diff File Expand all files Collapse all files

12
ortools/math_opt/core/BUILD.bazel
View File

@@ -193,3 +193,15 @@ cc_library(
"@com_google_absl//absl/types:span",
],
)
cc_library(
name = "inverted_bounds",
srcs = ["inverted_bounds.cc"],
hdrs = ["inverted_bounds.h"],
deps = [
"//ortools/base:status_builder",
"@com_google_absl//absl/status",
"@com_google_absl//absl/strings",
"@com_google_absl//absl/types:span",
],
)

58
ortools/math_opt/core/inverted_bounds.cc Normal file
View File

@@ -0,0 +1,58 @@
// Copyright 2010-2021 Google LLC
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#include "ortools/math_opt/core/inverted_bounds.h"
#include <cstdint>
#include <string>
#include <string_view>
#include <vector>
#include "absl/status/status.h"
#include "absl/strings/str_join.h"
#include "absl/types/span.h"
#include "ortools/base/status_builder.h"
namespace operations_research::math_opt {
absl::Status InvertedBounds::ToStatus() const {
if (empty()) {
return absl::OkStatus();
}
auto builder = util::InvalidArgumentErrorBuilder();
const auto format_bounds_ids = [&builder](const std::string_view name,
const std::vector<int64_t>& ids) {
if (ids.empty()) {
return;
}
builder << name << " with ids "
<< absl::StrJoin(absl::MakeSpan(ids).subspan(0, kMaxInvertedBounds),
",");
if (ids.size() > kMaxInvertedBounds) {
builder << "...";
}
};
format_bounds_ids("variables", variables);
if (!variables.empty() && !linear_constraints.empty()) {
builder << " and ";
}
format_bounds_ids("linear constraints", linear_constraints);
builder << " have lower_bound > upper_bound";
return builder;
}
} // namespace operations_research::math_opt

50
ortools/math_opt/core/inverted_bounds.h Normal file
View File

@@ -0,0 +1,50 @@
// Copyright 2010-2021 Google LLC
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#ifndef OR_TOOLS_MATH_OPT_CORE_INVERTED_BOUNDS_H_
#define OR_TOOLS_MATH_OPT_CORE_INVERTED_BOUNDS_H_
#include <cstdint>
#include <vector>
#include "absl/status/status.h"
namespace operations_research::math_opt {
// The maximum number of variables/constraints with inverted bounds to report.
constexpr std::size_t kMaxInvertedBounds = 10;
// The ids of the variables and linear constraints with inverted bounds
// (lower_bounds > upper_bounds).
//
// This is used internally by solvers to return an error on Solve() when bounds
// are inverted.
struct InvertedBounds {
// Returns true if this object contains no variable/constraint ids.
bool empty() const { return variables.empty() && linear_constraints.empty(); }
// Returns an error listing at most kMaxInvertedBounds variables and linear
// constraints ids (kMaxInvertedBounds of each). Returns OK status if this
// object is empty.
absl::Status ToStatus() const;
// Ids of the variables with inverted bounds.
std::vector<int64_t> variables;
// Ids of the linear constraints with inverted bounds.
std::vector<int64_t> linear_constraints;
};
} // namespace operations_research::math_opt
#endif // OR_TOOLS_MATH_OPT_CORE_INVERTED_BOUNDS_H_

2
ortools/math_opt/core/math_opt_proto_utils.h
View File

@@ -99,7 +99,7 @@ class SparseVectorFilterPredicate {
// Invariant: next input id must be >= next_input_id_lower_bound_.
//
// The initial value is 0 since all ids are expected to be non-negative.
int next_input_id_lower_bound_ = 0;
int64_t next_input_id_lower_bound_ = 0;
#endif // NDEBUG
};

4
ortools/math_opt/core/model_storage.cc
View File

@@ -748,7 +748,7 @@ std::optional<ModelUpdateProto> ModelStorage::ExportModelUpdate(
// * Use n-way merge here if the performances justify it.
const auto merge = std::make_shared<ModelUpdateProto>();
for (const auto& update : data->updates) {
MergeIntoUpdate(/*from=*/*update, /*into=*/*merge);
MergeIntoUpdate(/*from_new=*/*update, /*into_old=*/*merge);
}
// Push the merge to all trackers that have the same checkpoint (including
@@ -762,7 +762,7 @@ std::optional<ModelUpdateProto> ModelStorage::ExportModelUpdate(
const std::optional<ModelUpdateProto> pending_update =
ExportSharedModelUpdate();
if (pending_update) {
MergeIntoUpdate(/*from=*/*pending_update, /*into=*/update);
MergeIntoUpdate(/*from_new=*/*pending_update, /*into_old=*/update);
}
// Named returned value optimization (NRVO) does not apply here since the

2
ortools/math_opt/core/model_storage.h
View File

@@ -140,7 +140,7 @@ DEFINE_STRONG_INT_TYPE(UpdateTrackerId, int64_t);
// the modifications since the previous call to ModelStorage::Checkpoint(). Note
// that, for newly initialized models, before the first checkpoint, there is no
// additional memory overhead from tracking changes. See
// g3doc/ortools/math_opt/g3doc/model_building_complexity.md
// docs/ortools/math_opt/docs/model_building_complexity.md
// for details.
//
// On bad input:

2
ortools/math_opt/core/model_update_merge.h
View File

@@ -180,7 +180,7 @@ void UpdateNewElementProperty(
int updates_i = 0;
for (int i = 0; i < ids.size(); ++i) {
const int id = ids[i];
const int64_t id = ids[i];
while (deleted_i < deleted.size() && deleted[deleted_i] < id) {
++deleted_i;

16
ortools/math_opt/core/non_streamable_solver_init_arguments.h
View File

@@ -21,12 +21,12 @@
namespace operations_research {
namespace math_opt {
struct NonStreamableBiscoInitArguments;
struct NonStreamableCpSatInitArguments;
struct NonStreamableGScipInitArguments;
struct NonStreamableGlopInitArguments;
struct NonStreamableGlpkInitArguments;
struct NonStreamableGurobiInitArguments;
struct NonStreamablePdlpInitArguments;
// Interface for solver specific parameters used at the solver instantiation
// that can't be streamed (for example instances of C/C++ types that only exist
@@ -51,13 +51,6 @@ struct NonStreamableSolverInitArguments {
// Returns the type of solver that the implementation is for.
virtual SolverTypeProto solver_type() const = 0;
// Returns this for the NonStreamableBiscoInitArguments class, nullptr for
// other classes.
virtual const NonStreamableBiscoInitArguments*
ToNonStreamableBiscoInitArguments() const {
return nullptr;
}
// Returns this for the NonStreamableCpSatInitArguments class, nullptr for
// other classes.
virtual const NonStreamableCpSatInitArguments*
@@ -93,6 +86,13 @@ struct NonStreamableSolverInitArguments {
return nullptr;
}
// Returns this for the NonStreamablePdlpInitArguments class, nullptr for
// other classes.
virtual const NonStreamablePdlpInitArguments*
ToNonStreamablePdlpInitArguments() const {
return nullptr;
}
// Return a copy of this.
//
// The NonStreamableSolverInitArgumentsHelper implements this automatically

2
ortools/math_opt/cpp/BUILD.bazel
View File

@@ -242,6 +242,7 @@ cc_library(
deps = [
":streamable_solver_init_arguments",
"//ortools/math_opt/core:non_streamable_solver_init_arguments",
"@com_google_absl//absl/status:statusor",
],
)
@@ -289,6 +290,7 @@ cc_library(
"//ortools/math_opt:parameters_cc_proto",
"//ortools/math_opt/cpp:enums",
"//ortools/math_opt/solvers:gurobi_cc_proto",
"@com_google_absl//absl/status:statusor",
],
)

6
ortools/math_opt/cpp/id_map.h
View File

@@ -163,7 +163,7 @@ class IdMap {
inline std::pair<iterator, bool> try_emplace(const K& k, Args&&... args);
// Returns the number of elements erased (zero or one).
inline int erase(const K& k);
inline size_type erase(const K& k);
// In STL erase(const_iterator) and erase(iterator) both return an
// iterator. But flat_hash_map instead has void return types. So here we also
// use void.
@@ -442,9 +442,9 @@ std::pair<typename IdMap<K, V>::iterator, bool> IdMap<K, V>::try_emplace(
}
template <typename K, typename V>
int IdMap<K, V>::erase(const K& k) {
typename IdMap<K, V>::size_type IdMap<K, V>::erase(const K& k) {
CheckModel(k);
const int ret = map_.erase(k.typed_id());
const size_type ret = map_.erase(k.typed_id());
if (map_.empty()) {
storage_ = nullptr;
}

6
ortools/math_opt/cpp/id_set.h
View File

@@ -125,7 +125,7 @@ class IdSet {
inline std::pair<const_iterator, bool> emplace(const K& k);
// Returns the number of elements erased (zero or one).
inline int erase(const K& k);
inline size_type erase(const K& k);
// In STL erase(const_iterator) returns an iterator. But flat_hash_set instead
// has void return types. So here we also use void.
inline void erase(const_iterator pos);
@@ -289,9 +289,9 @@ std::pair<typename IdSet<K>::const_iterator, bool> IdSet<K>::emplace(
}
template <typename K>
int IdSet<K>::erase(const K& k) {
typename IdSet<K>::size_type IdSet<K>::erase(const K& k) {
CheckModel(k);
const int ret = set_.erase(k.typed_id());
const size_type ret = set_.erase(k.typed_id());
if (set_.empty()) {
storage_ = nullptr;
}

1
ortools/math_opt/cpp/matchers.cc
View File

@@ -40,7 +40,6 @@ using ::testing::AnyOf;
using ::testing::AnyOfArray;
using ::testing::Contains;
using ::testing::DoubleNear;
using ::testing::Eq;
using ::testing::ExplainMatchResult;
using ::testing::Field;
using ::testing::IsEmpty;

2
ortools/math_opt/cpp/model.cc
View File

@@ -209,7 +209,7 @@ QuadraticExpression Model::ObjectiveAsQuadraticExpression() const {
ModelProto Model::ExportModel() const { return storage()->ExportModel(); }
std::unique_ptr<UpdateTracker> Model::NewUpdateTracker() const {
std::unique_ptr<UpdateTracker> Model::NewUpdateTracker() {
return std::make_unique<UpdateTracker>(storage_);
}

159
ortools/math_opt/cpp/model.h
View File

@@ -46,65 +46,29 @@ namespace math_opt {
// x in {0.0, 1.0}
// y in [0.0, 2.5]
//
// using ::operations_research::math_opt::LinearConstraint;
// using ::operations_research::math_opt::Model;
// using ::operations_research::math_opt::SolveResult;
// using ::operations_research::math_opt::SolveParameters;
// using ::operations_research::math_opt::SolveResultProto;
// using ::operations_research::math_opt::Variable;
// using ::operations_research::math_opt::SolverType;
//
// Version 1:
//
// Model model("my_model");
// const Variable x = model.AddBinaryVariable("x");
// const Variable y = model.AddContinuousVariable(0.0, 2.5, "y");
// const LinearConstraint c = model.AddLinearConstraint(
// -std::numeric_limits<double>::infinity(), 1.5, "c");
// model.set_coefficient(c, x, 1.0);
// model.set_coefficient(c, y, 1.0);
// model.set_objective_coefficient(x, 2.0);
// model.set_objective_coefficient(y, 1.0);
// model.set_maximize();
// const SolveResult result = Solve(
// model, SolverType::kGscip, SolveParametersProto()).value();
// for (const auto& warning : result.warnings) {
// std::cerr << "Solver warning: " << warning << std::endl;
// }
// CHECK_EQ(result.termination.reason, TerminationReason::kOptimal)
// << result.termination_detail;
// // The following code will print:
// // objective value: 2.5
// // value for variable x: 1
// std::cout << "objective value: " << result.objective_value()
// << "\nvalue for variable x: " << result.variable_values().at(x)
// << std::endl;
//
// Version 2 (with linear expressions):
//
// Model model("my_model");
// const Variable x = model.AddBinaryVariable("x");
// const Variable y = model.AddContinuousVariable(0.0, 2.5, "y");
// // We can directly use linear combinations of variables ...
// model.AddLinearConstraint(x + y <= 1.5, "c");
// // ... or build them incrementally.
// LinearExpression objective_expression;
// objective_expression += 2*x;
// objective_expression += y;
// model.Maximize(objective_expression);
// const SolveResult result = Solve(
// model, SolverType::kGscip, SolveParametersProto()).value();
// for (const auto& warning : result.warnings) {
// std::cerr << "Solver warning: " << warning << std::endl;
// }
// CHECK_EQ(result.termination.reason, TerminationReason::kOptimal)
// << result.termination_detail;
// // The following code will print:
// // objective value: 2.5
// // value for variable x: 1
// std::cout << "objective value: " << result.objective_value()
// << "\nvalue for variable x: " << result.variable_values().at(x)
// << std::endl;
// math_opt::Model model("my_model");
// const math_opt::Variable x = model.AddBinaryVariable("x");
// const math_opt::Variable y = model.AddContinuousVariable(0.0, 2.5, "y");
// // We can directly use linear combinations of variables ...
// model.AddLinearConstraint(x + y <= 1.5, "c");
// // ... or build them incrementally.
// math_opt::LinearExpression objective_expression;
// objective_expression += 2 * x;
// objective_expression += y;
// model.Maximize(objective_expression);
// ASSIGN_OR_RETURN(const math_opt::SolveResult result,
// Solve(model, math_opt::SolverType::kGscip));
// switch (result.termination.reason) {
// case math_opt::TerminationReason::kOptimal:
// case math_opt::TerminationReason::kFeasible:
// std::cout << "objective value: " << result.objective_value() << std::endl
// << "value for variable x: " << result.variable_values().at(x)
// << std::endl;
// return absl::OkStatus();
// default:
// return util::InternalErrorBuilder()
// << "model failed to solve: " << result.termination;
// }
//
// Memory model:
//
@@ -165,10 +129,12 @@ class Model {
// That said, the Variable and LinearConstraint reference objects are model
// specific. Hence the ones linked to the original model must NOT be used with
// the clone. The Variable and LinearConstraint reference objects for the
// clone can be obtained via Variables() and LinearConstraints(). One can also
// use SortedVariables() and SortedLinearConstraints() that will return (until
// one of the two models is modified) the variables and constraints in the
// same order for the two models and provide a one-to-one correspondence.
// clone can be obtained using:
// * the variable() and linear_constraint() methods on the ids from the old
// Variable and LinearConstraint objects.
// * in increasing id order using SortedVariables() and
// SortedLinearConstraints()
// * in an arbitrary order using Variables() and LinearConstraints().
//
// Note that the returned model does not have any update tracker.
std::unique_ptr<Model> Clone() const;
@@ -209,10 +175,19 @@ class Model {
// The returned id of the next call to AddVariable.
//
// Equal to the number of variables created.
inline int next_variable_id() const;
inline int64_t next_variable_id() const;
// Returns true if this id has been created and not yet deleted.
inline bool has_variable(int id) const;
inline bool has_variable(int64_t id) const;
// Returns true if this id has been created and not yet deleted.
inline bool has_variable(VariableId id) const;
// Will CHECK if has_variable(id) is false.
inline Variable variable(int64_t id) const;
// Will CHECK if has_variable(id) is false.
inline Variable variable(VariableId id) const;
// Returns the variable name.
inline const std::string& name(Variable variable) const;
@@ -292,10 +267,19 @@ class Model {
// The returned id of the next call to AddLinearConstraint.
//
// Equal to the number of linear constraints created.
inline int next_linear_constraint_id() const;
inline int64_t next_linear_constraint_id() const;
// Returns true if this id has been created and not yet deleted.
inline bool has_linear_constraint(int id) const;
inline bool has_linear_constraint(int64_t id) const;
// Returns true if this id has been created and not yet deleted.
inline bool has_linear_constraint(LinearConstraintId id) const;
// Will CHECK if has_linear_constraint(id) is false.
inline LinearConstraint linear_constraint(int64_t id) const;
// Will CHECK if has_linear_constraint(id) is false.
inline LinearConstraint linear_constraint(LinearConstraintId id) const;
// Returns the linear constraint name.
inline const std::string& name(LinearConstraint constraint) const;
@@ -447,7 +431,7 @@ class Model {
// (variables, constraints, ...). The user is expected to use proper
// synchronization primitive to serialize changes to the model and the use of
// this method.
std::unique_ptr<UpdateTracker> NewUpdateTracker() const;
std::unique_ptr<UpdateTracker> NewUpdateTracker();
// Apply the provided update to this model. Returns a failure if the update is
// not valid.
@@ -534,12 +518,25 @@ void Model::DeleteVariable(const Variable variable) {
int Model::num_variables() const { return storage()->num_variables(); }
int Model::next_variable_id() const {
int64_t Model::next_variable_id() const {
return storage()->next_variable_id().value();
}
bool Model::has_variable(const int id) const {
return storage()->has_variable(VariableId(id));
bool Model::has_variable(const int64_t id) const {
return has_variable(VariableId(id));
}
bool Model::has_variable(const VariableId id) const {
return storage()->has_variable(id);
}
Variable Model::variable(const int64_t id) const {
return variable(VariableId(id));
}
Variable Model::variable(const VariableId id) const {
CHECK(has_variable(id)) << "No variable with id: " << id.value();
return Variable(storage(), id);
}
const std::string& Model::name(const Variable variable) const {
@@ -604,12 +601,26 @@ int Model::num_linear_constraints() const {
return storage()->num_linear_constraints();
}
int Model::next_linear_constraint_id() const {
int64_t Model::next_linear_constraint_id() const {
return storage()->next_linear_constraint_id().value();
}
bool Model::has_linear_constraint(const int id) const {
return storage()->has_linear_constraint(LinearConstraintId(id));
bool Model::has_linear_constraint(const int64_t id) const {
return has_linear_constraint(LinearConstraintId(id));
}
bool Model::has_linear_constraint(const LinearConstraintId id) const {
return storage()->has_linear_constraint(id);
}
LinearConstraint Model::linear_constraint(const int64_t id) const {
return linear_constraint(LinearConstraintId(id));
}
LinearConstraint Model::linear_constraint(const LinearConstraintId id) const {
CHECK(has_linear_constraint(id))
<< "No linear constraint with id: " << id.value();
return LinearConstraint(storage(), id);
}
const std::string& Model::name(const LinearConstraint constraint) const {

22
ortools/math_opt/cpp/parameters.cc
View File

@@ -151,6 +151,9 @@ SolveParametersProto SolveParameters::Proto() const {
if (iteration_limit.has_value()) {
result.set_iteration_limit(*iteration_limit);
}
if (node_limit.has_value()) {
result.set_node_limit(*node_limit);
}
if (cutoff_limit.has_value()) {
result.set_cutoff_limit(*cutoff_limit);
}
@@ -169,11 +172,11 @@ SolveParametersProto SolveParameters::Proto() const {
if (random_seed.has_value()) {
result.set_random_seed(*random_seed);
}
if (relative_gap_limit.has_value()) {
result.set_relative_gap_limit(*relative_gap_limit);
if (relative_gap_tolerance.has_value()) {
result.set_relative_gap_tolerance(*relative_gap_tolerance);
}
if (absolute_gap_limit.has_value()) {
result.set_absolute_gap_limit(*absolute_gap_limit);
if (absolute_gap_tolerance.has_value()) {
result.set_absolute_gap_tolerance(*absolute_gap_tolerance);
}
result.set_lp_algorithm(EnumToProto(lp_algorithm));
result.set_presolve(EnumToProto(presolve));
@@ -200,6 +203,9 @@ absl::StatusOr<SolveParameters> SolveParameters::FromProto(
if (proto.has_iteration_limit()) {
result.iteration_limit = proto.iteration_limit();
}
if (proto.has_node_limit()) {
result.node_limit = proto.node_limit();
}
if (proto.has_cutoff_limit()) {
result.cutoff_limit = proto.cutoff_limit();
}
@@ -218,11 +224,11 @@ absl::StatusOr<SolveParameters> SolveParameters::FromProto(
if (proto.has_random_seed()) {
result.random_seed = proto.random_seed();
}
if (proto.has_absolute_gap_limit()) {
result.absolute_gap_limit = proto.absolute_gap_limit();
if (proto.has_absolute_gap_tolerance()) {
result.absolute_gap_tolerance = proto.absolute_gap_tolerance();
}
if (proto.has_relative_gap_limit()) {
result.relative_gap_limit = proto.relative_gap_limit();
if (proto.has_relative_gap_tolerance()) {
result.relative_gap_tolerance = proto.relative_gap_tolerance();
}
result.lp_algorithm = EnumFromProto(proto.lp_algorithm());
result.presolve = EnumFromProto(proto.presolve());

53
ortools/math_opt/cpp/parameters.h
View File

@@ -189,19 +189,19 @@ struct SolveParameters {
// Limit on the iterations of the underlying algorithm (e.g. simplex pivots).
// The specific behavior is dependent on the solver and algorithm used, but
// should result in a deterministic solve limit.
// TODO(b/195295177): suggest node_limit as an alternative when it's added
// often can give a deterministic solve limit (further configuration may be
// needed, e.g. one thread).
//
// Typically supported by LP, QP, and MIP solvers, but for MIP solvers see
// also node_limit.
std::optional<int64_t> iteration_limit;
// Optimality tolerances (primarily) for MIP solvers. The absolute GAP of a
// feasible solution is the distance between its objective value and a dual
// bound (e.g. an upper bound on the optimal value for maximization problems).
// The relative GAP is a solver-dependent scaled version of the absolute GAP
// (e.g. it could be the relative GAP divided by the objective value of the
// feasible solution if this is non-zero). Solvers consider a solution optimal
// if its GAPs are below these limits (most solvers use both versions).
std::optional<double> relative_gap_limit;
std::optional<double> absolute_gap_limit;
// Limit on the number of subproblems solved in enumerative search (e.g.
// branch and bound). For many solvers this can be used to deterministically
// limit computation (further configuration may be needed, e.g. one thread).
//
// Typically for MIP solvers, see also iteration_limit.
std::optional<int64_t> node_limit;
// The solver stops early if it can prove there are no primal solutions at
// least as good as cutoff.
@@ -217,8 +217,7 @@ struct SolveParameters {
std::optional<double> cutoff_limit;
// The solver stops early as soon as it finds a solution at least this good,
// with termination reason kFeasible or kNoSolutionFound and limit kObjective.
// TODO(b/214567536): maybe it should only be kFeasible.
// with termination reason kFeasible and limit kObjective.
std::optional<double> objective_limit;
// The solver stops early as soon as it proves the best bound is at least this
@@ -258,6 +257,34 @@ struct SolveParameters {
// MAX(0, MIN(MAX_VALID_VALUE_FOR_SOLVER, random_seed)).
std::optional<int32_t> random_seed;
// An absolute optimality tolerance (primarily) for MIP solvers.
//
// The absolute GAP is the absolute value of the difference between:
// * the objective value of the best feasible solution found,
// * the dual bound produced by the search.
// The solver can stop once the absolute GAP is at most absolute_gap_tolerance
// (when set), and return TerminationReason::kOptimal.
//
// Must be >= 0 if set.
//
// See also relative_gap_tolerance.
std::optional<double> absolute_gap_tolerance;
// A relative optimality tolerance (primarily) for MIP solvers.
//
// The relative GAP is a normalized version of the absolute GAP (defined on
// absolute_gap_tolerance), where the normalization is solver-dependent, e.g.
// the absolute GAP divided by the objective value of the best feasible
// solution found.
//
// The solver can stop once the relative GAP is at most relative_gap_tolerance
// (when set), and return TerminationReason::kOptimal.
//
// Must be >= 0 if set.
//
// See also absolute_gap_tolerance.
std::optional<double> relative_gap_tolerance;
// The algorithm for solving a linear program. If nullopt, use the solver
// default algorithm.
//

29
ortools/math_opt/cpp/solve.cc
View File

@@ -102,14 +102,18 @@ absl::StatusOr<SolveResult> Solve(const Model& model,
}
absl::StatusOr<std::unique_ptr<IncrementalSolver>> IncrementalSolver::New(
Model& model, const SolverType solver_type, SolverInitArguments arguments) {
std::unique_ptr<UpdateTracker> update_tracker = model.NewUpdateTracker();
ASSIGN_OR_RETURN(
std::unique_ptr<Solver> solver,
Solver::New(EnumToProto(solver_type), update_tracker->ExportModel(),
ToSolverInitArgs(arguments)));
Model* const model, const SolverType solver_type,
SolverInitArguments arguments) {
if (model == nullptr) {
return absl::InvalidArgumentError("input model can't be null");
}
std::unique_ptr<UpdateTracker> update_tracker = model->NewUpdateTracker();
ASSIGN_OR_RETURN(const ModelProto model_proto, update_tracker->ExportModel());
ASSIGN_OR_RETURN(std::unique_ptr<Solver> solver,
Solver::New(EnumToProto(solver_type), model_proto,
ToSolverInitArgs(arguments)));
return absl::WrapUnique<IncrementalSolver>(
new IncrementalSolver(solver_type, std::move(arguments), model.storage(),
new IncrementalSolver(solver_type, std::move(arguments), model->storage(),
std::move(update_tracker), std::move(solver)));
}
@@ -131,21 +135,22 @@ absl::StatusOr<SolveResult> IncrementalSolver::Solve(
}
absl::StatusOr<IncrementalSolver::UpdateResult> IncrementalSolver::Update() {
std::optional<ModelUpdateProto> model_update =
update_tracker_->ExportModelUpdate();
ASSIGN_OR_RETURN(std::optional<ModelUpdateProto> model_update,
update_tracker_->ExportModelUpdate());
if (!model_update) {
return UpdateResult(true, std::move(model_update));
}
ASSIGN_OR_RETURN(const bool did_update, solver_->Update(*model_update));
update_tracker_->Checkpoint();
RETURN_IF_ERROR(update_tracker_->Checkpoint());
if (did_update) {
return UpdateResult(true, std::move(model_update));
}
ASSIGN_OR_RETURN(solver_, Solver::New(EnumToProto(solver_type_),
update_tracker_->ExportModel(),
ASSIGN_OR_RETURN(const ModelProto model_proto,
update_tracker_->ExportModel());
ASSIGN_OR_RETURN(solver_, Solver::New(EnumToProto(solver_type_), model_proto,
ToSolverInitArgs(init_args_)));
return UpdateResult(false, std::move(model_update));

18
ortools/math_opt/cpp/solve.h
View File

@@ -75,15 +75,15 @@ absl::StatusOr<SolveResult> Solve(const Model& model, SolverType solver_type,
// MIPs. In both cases, even if the solver supports the model changes,
// incremental solve may actually be slower.
//
// The New() function instantiates the solver and setup it from the current
// state of the Model. Calling Solve() will update the underlying solver with
// latest model changes and solve this model.
// The New() function instantiates the solver, setup it from the current state
// of the Model and register on it to listen to changes. Calling Solve() will
// update the underlying solver with latest model changes and solve this model.
//
// Usage:
// Model model = ...;
// ASSIGN_OR_RETURN(
// const std::unique_ptr<IncrementalSolver> incremental_solve,
// IncrementalSolver::New(model, SolverType::kXxx));
// IncrementalSolver::New(&model, SolverType::kXxx));
//
// ASSIGN_OR_RETURN(const SolveResult result1, incremental_solve->Solve());
//
@@ -94,6 +94,11 @@ absl::StatusOr<SolveResult> Solve(const Model& model, SolverType solver_type,
//
// ...
//
// Lifecycle: The IncrementalSolver is only keeping a std::weak_ref on Model's
// internal data and thus it returns an error if Update() or Solve() are called
// after the Model has been destroyed. It is fine to destroy the
// IncrementalSolver after the associated Model though.
//
// Thread-safety: The New(), Solve() and Update() methods must not be called
// while modifying the Model() (adding variables...). The user is expected to
// use proper synchronization primitives to serialize changes to the model and
@@ -134,9 +139,10 @@ class IncrementalSolver {
// The returned IncrementalSolver keeps a copy of `arguments`. Thus the
// content of arguments.non_streamable (for example pointers to solver
// specific struct) must be valid until the destruction of the
// IncrementalSolver.
// IncrementalSolver. It also registers on the Model to keep track of updates
// (see class documentation for details).
static absl::StatusOr<std::unique_ptr<IncrementalSolver>> New(
Model& model, SolverType solver_type, SolverInitArguments arguments = {});
Model* model, SolverType solver_type, SolverInitArguments arguments = {});
// Updates the underlying solver with latest model changes and runs the solve.
//

30
ortools/math_opt/cpp/streamable_solver_init_arguments.cc
View File

@@ -16,6 +16,7 @@
#include <optional>
#include <type_traits>
#include "absl/status/statusor.h"
#include "ortools/math_opt/parameters.pb.h"
#include "ortools/math_opt/solvers/gurobi.pb.h"
@@ -31,6 +32,16 @@ GurobiInitializerProto::ISVKey GurobiISVKey::Proto() const {
return isv_key_proto;
}
GurobiISVKey GurobiISVKey::FromProto(
const GurobiInitializerProto::ISVKey& key_proto) {
return GurobiISVKey{
.name = key_proto.name(),
.application_name = key_proto.application_name(),
.expiration = key_proto.expiration(),
.key = key_proto.key(),
};
}
GurobiInitializerProto StreamableGurobiInitArguments::Proto() const {
GurobiInitializerProto params_proto;
@@ -41,6 +52,15 @@ GurobiInitializerProto StreamableGurobiInitArguments::Proto() const {
return params_proto;
}
StreamableGurobiInitArguments StreamableGurobiInitArguments::FromProto(
const GurobiInitializerProto& args_proto) {
StreamableGurobiInitArguments args;
if (args_proto.has_isv_key()) {
args.isv_key = GurobiISVKey::FromProto(args_proto.isv_key());
}
return args;
}
SolverInitializerProto StreamableSolverInitArguments::Proto() const {
SolverInitializerProto params_proto;
@@ -51,5 +71,15 @@ SolverInitializerProto StreamableSolverInitArguments::Proto() const {
return params_proto;
}
absl::StatusOr<StreamableSolverInitArguments>
StreamableSolverInitArguments::FromProto(
const SolverInitializerProto& args_proto) {
StreamableSolverInitArguments args;
if (args_proto.has_gurobi()) {
args.gurobi = StreamableGurobiInitArguments::FromProto(args_proto.gurobi());
}
return args;
}
} // namespace math_opt
} // namespace operations_research

14
ortools/math_opt/cpp/streamable_solver_init_arguments.h
View File

@@ -24,7 +24,7 @@
#include <optional>
#include <string>
#include "ortools/math_opt/cpp/enums.h"
#include "absl/status/statusor.h"
#include "ortools/math_opt/parameters.pb.h"
#include "ortools/math_opt/solvers/gurobi.pb.h"
@@ -52,10 +52,12 @@ struct StreamableGlpkInitArguments {};
struct GurobiISVKey {
std::string name;
std::string application_name;
int64_t expiration = 0;
int32_t expiration = 0;
std::string key;
GurobiInitializerProto::ISVKey Proto() const;
static GurobiISVKey FromProto(
const GurobiInitializerProto::ISVKey& key_proto);
};
// Streamable Gurobi specific parameters for solver instantiation.
@@ -66,6 +68,10 @@ struct StreamableGurobiInitArguments {
// Returns the proto corresponding to these parameters.
GurobiInitializerProto Proto() const;
// Parses the proto corresponding to these parameters.
static StreamableGurobiInitArguments FromProto(
const GurobiInitializerProto& args_proto);
};
// Solver initialization parameters that can be streamed to be exchanged with
@@ -83,6 +89,10 @@ struct StreamableSolverInitArguments {
// Returns the proto corresponding to these parameters.
SolverInitializerProto Proto() const;
// Parses the proto corresponding to these parameters.
static absl::StatusOr<StreamableSolverInitArguments> FromProto(
const SolverInitializerProto& args_proto);
};
} // namespace math_opt

22
ortools/math_opt/cpp/update_tracker.cc
View File

@@ -16,6 +16,8 @@
#include <memory>
#include <optional>
#include "absl/status/status.h"
#include "absl/status/statusor.h"
#include "ortools/base/logging.h"
#include "ortools/math_opt/core/model_storage.h"
#include "ortools/math_opt/model.pb.h"
@@ -37,21 +39,29 @@ UpdateTracker::UpdateTracker(const std::shared_ptr<ModelStorage>& storage)
: storage_(ABSL_DIE_IF_NULL(storage)),
update_tracker_(storage->NewUpdateTracker()) {}
std::optional<ModelUpdateProto> UpdateTracker::ExportModelUpdate() {
absl::StatusOr<std::optional<ModelUpdateProto>>
UpdateTracker::ExportModelUpdate() {
const std::shared_ptr<ModelStorage> storage = storage_.lock();
CHECK(storage != nullptr) << internal::kModelIsDestroyed;
if (storage == nullptr) {
return absl::InvalidArgumentError(internal::kModelIsDestroyed);
}
return storage->ExportModelUpdate(update_tracker_);
}
void UpdateTracker::Checkpoint() {
absl::Status UpdateTracker::Checkpoint() {
const std::shared_ptr<ModelStorage> storage = storage_.lock();
CHECK(storage != nullptr) << internal::kModelIsDestroyed;
if (storage == nullptr) {
return absl::InvalidArgumentError(internal::kModelIsDestroyed);
}
storage->Checkpoint(update_tracker_);
return absl::OkStatus();
}
ModelProto UpdateTracker::ExportModel() const {
absl::StatusOr<ModelProto> UpdateTracker::ExportModel() const {
const std::shared_ptr<ModelStorage> storage = storage_.lock();
CHECK(storage != nullptr) << internal::kModelIsDestroyed;
if (storage == nullptr) {
return absl::InvalidArgumentError(internal::kModelIsDestroyed);
}
return storage->ExportModel();
}

24
ortools/math_opt/cpp/update_tracker.h
View File

@@ -17,6 +17,8 @@
#include <memory>
#include <optional>
#include "absl/status/status.h"
#include "absl/status/statusor.h"
#include "absl/strings/string_view.h"
#include "ortools/math_opt/core/model_storage.h"
#include "ortools/math_opt/model.pb.h"
@@ -55,9 +57,9 @@ namespace math_opt {
// model.AddVariable(0.0, 1.0, true, "y");
// model.set_maximize(true);
//
// const std::optional<ModelUpdateProto> update_proto =
// update_tracker.ExportModelUpdate();
// update_tracker.Checkpoint();
// ASSIGN_OR_RETURN(const std::optional<ModelUpdateProto> update_proto,
// update_tracker.ExportModelUpdate());
// RETURN_IF_ERROR(update_tracker.Checkpoint());
//
// if (update_proto) {
// ... use *update_proto here ...
@@ -73,18 +75,24 @@ class UpdateTracker {
// Returns a proto representation of the changes to the model since the most
// recent checkpoint (i.e. last time Checkpoint() was called); nullopt if
// the update would have been empty.
std::optional<ModelUpdateProto> ExportModelUpdate();
//
// If fails if the Model has been destroyed.
absl::StatusOr<std::optional<ModelUpdateProto>> ExportModelUpdate();
// Uses the current model state as the starting point to calculate the
// ModelUpdateProto next time ExportModelUpdate() is called.
void Checkpoint();
//
// If fails if the Model has been destroyed.
absl::Status Checkpoint();
// Returns a proto representation of the whole model.
//
// This is a shortcut method that is equivalent to calling
// Model::ExportModel(). It is there so that users of the UpdateTracker
// can avoid having to keep a reference to the Model model.
ModelProto ExportModel() const;
//
// If fails if the Model has been destroyed.
absl::StatusOr<ModelProto> ExportModel() const;
private:
const std::weak_ptr<ModelStorage> storage_;
@@ -93,10 +101,10 @@ class UpdateTracker {
namespace internal {
// The CHECK message used when a function of UpdateTracker is called after the
// The failure message used when a function of UpdateTracker is called after the
// destruction of the model..
constexpr absl::string_view kModelIsDestroyed =
"Can't call this function after the associated model has been destroyed.";
"can't call this function after the associated model has been destroyed";
} // namespace internal

16
ortools/math_opt/io/BUILD.bazel
View File

@@ -22,3 +22,19 @@ cc_library(
"@com_google_absl//absl/status:statusor",
],
)
cc_library(
name = "mps_converter",
srcs = ["mps_converter.cc"],
hdrs = ["mps_converter.h"],
deps = [
":proto_converter",
"//ortools/linear_solver:linear_solver_cc_proto",
"//ortools/lp_data:mps_reader",
"//ortools/math_opt:model_cc_proto",
"//ortools/util:file_util",
"@com_google_absl//absl/status",
"@com_google_absl//absl/status:statusor",
"@com_google_absl//absl/strings",
],
)

35
ortools/math_opt/io/mps_converter.cc Normal file
View File

@@ -0,0 +1,35 @@
// Copyright 2010-2021 Google LLC
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#include "ortools/math_opt/io/mps_converter.h"
#include <string>
#include "absl/status/statusor.h"
#include "absl/strings/string_view.h"
#include "ortools/linear_solver/linear_solver.pb.h"
#include "ortools/lp_data/mps_reader.h"
#include "ortools/math_opt/io/proto_converter.h"
#include "ortools/math_opt/model.pb.h"
#include "ortools/util/file_util.h"
namespace operations_research::math_opt {
absl::StatusOr<ModelProto> ReadMpsFile(const absl::string_view filename) {
glop::MPSReader mps_reader;
MPModelProto mp_model;
RETURN_IF_ERROR(mps_reader.ParseFile(std::string(filename), &mp_model));
return MPModelProtoToMathOptModel(mp_model);
}
} // namespace operations_research::math_opt

35
ortools/math_opt/io/mps_converter.h Normal file
View File

@@ -0,0 +1,35 @@
// Copyright 2010-2021 Google LLC
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#ifndef OR_TOOLS_MATH_OPT_IO_MPS_CONVERTER_H_
#define OR_TOOLS_MATH_OPT_IO_MPS_CONVERTER_H_
#include <string>
#include "absl/status/statusor.h"
#include "absl/strings/string_view.h"
#include "ortools/math_opt/model.pb.h"
namespace operations_research::math_opt {
// Reads an MPS file and converts it to a ModelProto (like MpsToModelProto
// above, but takes a file name instead of the file contents and reads the file.
//
//
// The file can be stored as plain text or gzipped (with the .gz extension).
//
absl::StatusOr<ModelProto> ReadMpsFile(absl::string_view filename);
} // namespace operations_research::math_opt
#endif // OR_TOOLS_MATH_OPT_IO_MPS_CONVERTER_H_

75
ortools/math_opt/parameters.proto
View File

@@ -152,29 +152,26 @@ message SolveParametersProto {
// Limit on the iterations of the underlying algorithm (e.g. simplex pivots).
// The specific behavior is dependent on the solver and algorithm used, but
// should result in a deterministic solve limit.
// TODO(b/195295177): suggest node_limit as an alternative when it's added
// often can give a deterministic solve limit (further configuration may be
// needed, e.g. one thread).
//
// Typically supported by LP, QP, and MIP solvers, but for MIP solvers see
// also node_limit.
optional int64 iteration_limit = 2;
// Optimality tolerances (primarily) for MIP solvers. The absolute GAP of a
// feasible solution is the distance between its objective value and a dual
// bound (e.g. an upper bound on the optimal value for maximization problems).
// The relative GAP is a solver-dependent scaled version of the absolute GAP
// (e.g. it could be the relative GAP divided by the objective value of the
// feasible solution if this is non-zero). Solvers consider a solution optimal
// if its GAPs are below these limits (most solvers use both versions).
// TODO(b/213603287): rename as relative_gap_tolerance and
// absolute_gap_tolerance.
optional double relative_gap_limit = 17;
optional double absolute_gap_limit = 18;
// Limit on the number of subproblems solved in enumerative search (e.g.
// branch and bound). For many solvers this can be used to deterministically
// limit computation (further configuration may be needed, e.g. one thread).
//
// Typically for MIP solvers, see also iteration_limit.
optional int64 node_limit = 24;
// The solver stops early if it can prove there are no primal solutions at
// least as good as cutoff.
//
// On an early stop, the solver returns termination reason
// LIMIT_NO_SOLUTION_FOUND and with limit CUTOFF and is not required to give
// any extra solution information. Has no effect on the return value if there
// is no early stop.
// On an early stop, the solver returns termination reason NO_SOLUTION_FOUND
// and with limit CUTOFF and is not required to give any extra solution
// information. Has no effect on the return value if there is no early stop.
//
// It is recommended that you use a tolerance if you want solutions with
// objective exactly equal to cutoff to be returned.
@@ -183,23 +180,21 @@ message SolveParametersProto {
optional double cutoff_limit = 20;
// The solver stops early as soon as it finds a solution at least this good,
// with termination reason LIMIT_FEASIBLE or LIMIT_NO_SOLUTION_FOUND and
// limit LIMIT_OBJECTIVE.
// TODO(b/214567536): maybe it should only be LIMIT_FEASIBLE.
// with termination reason FEASIBLE and limit OBJECTIVE.
optional double objective_limit = 21;
// The solver stops early as soon as it proves the best bound is at least this
// good, with termination reason LIMIT_FEASIBLE or LIMIT_NO_SOLUTION_FOUND and
// limit LIMIT_OBJECTIVE.
// good, with termination reason FEASIBLE or NO_SOLUTION_FOUND and limit
// OBJECTIVE.
//
// See the user guide for more details and a comparison with cutoff_limit.
optional double best_bound_limit = 22;
// The solver stops early after finding this many feasible solutions, with
// termination reason LIMIT_FEASIBLE and limit LIMIT_SOLUTION. Must be greater
// than zero if set. It is often used get the solver to stop on the first
// feasible solution found. Note that there is no guarantee on the objective
// value for any of the returned solutions.
// termination reason FEASIBLE and limit SOLUTION. Must be greater than zero
// if set. It is often used get the solver to stop on the first feasible
// solution found. Note that there is no guarantee on the objective value for
// any of the returned solutions.
//
// Solvers will typically not return more solutions than the solution limit,
// but this is not enforced by MathOpt, see also b/214041169.
@@ -234,6 +229,34 @@ message SolveParametersProto {
// MAX(0, MIN(MAX_VALID_VALUE_FOR_SOLVER, random_seed)).
optional int32 random_seed = 5;
// An absolute optimality tolerance (primarily) for MIP solvers.
//
// The absolute GAP is the absolute value of the difference between:
// * the objective value of the best feasible solution found,
// * the dual bound produced by the search.
// The solver can stop once the absolute GAP is at most absolute_gap_tolerance
// (when set), and return TERMINATION_REASON_OPTIMAL.
//
// Must be >= 0 if set.
//
// See also relative_gap_tolerance.
optional double absolute_gap_tolerance = 18;
// A relative optimality tolerance (primarily) for MIP solvers.
//
// The relative GAP is a normalized version of the absolute GAP (defined on
// absolute_gap_tolerance), where the normalization is solver-dependent, e.g.
// the absolute GAP divided by the objective value of the best feasible
// solution found.
//
// The solver can stop once the relative GAP is at most relative_gap_tolerance
// (when set), and return TERMINATION_REASON_OPTIMAL.
//
// Must be >= 0 if set.
//
// See also absolute_gap_tolerance.
optional double relative_gap_tolerance = 17;
// The algorithm for solving a linear program. If LP_ALGORITHM_UNSPECIFIED,
// use the solver default algorithm.
//

2
ortools/math_opt/result.proto
View File

@@ -83,7 +83,7 @@ message SolveStatsProto {
// minimization and larger for maximization) than best_primal_bound:
// * best_primal_bound is trivial (+inf for minimization and -inf
// maximization) when the solver does not claim to have such bound. This
// may happen for some solvers (e.g. bisco, typically continuous solvers)
// may happen for some solvers (e.g. PDLP, typically continuous solvers)
// even when returning optimal (solver could terminate with slightly
// infeasible primal solutions).
// * best_primal_bound can be closer to the optimal value than the objective

80
ortools/math_opt/samples/basic_example.cc
View File

@@ -16,82 +16,56 @@
#include <iostream>
#include <limits>
#include "absl/status/status.h"
#include "absl/status/statusor.h"
#include "ortools/base/init_google.h"
#include "ortools/base/logging.h"
#include "ortools/base/status_builder.h"
#include "ortools/base/status_macros.h"
#include "ortools/math_opt/cpp/math_opt.h"
namespace {
namespace math_opt = ::operations_research::math_opt;
// Model the problem:
// max 2.0 * x + y
// s.t. x + y <= 1.5
// x in {0.0, 1.0}
// y in [0.0, 2.5]
//
void SolveVersion1() {
using ::operations_research::math_opt::LinearConstraint;
using ::operations_research::math_opt::Model;
using ::operations_research::math_opt::SolveResult;
using ::operations_research::math_opt::SolverType;
using ::operations_research::math_opt::TerminationReason;
using ::operations_research::math_opt::Variable;
Model model("my_model");
const Variable x = model.AddBinaryVariable("x");
const Variable y = model.AddContinuousVariable(0.0, 2.5, "y");
const LinearConstraint c = model.AddLinearConstraint(
-std::numeric_limits<double>::infinity(), 1.5, "c");
model.set_coefficient(c, x, 1.0);
model.set_coefficient(c, y, 1.0);
model.set_objective_coefficient(x, 2.0);
model.set_objective_coefficient(y, 1.0);
model.set_maximize();
const SolveResult result = Solve(model, SolverType::kGscip).value();
CHECK_EQ(result.termination.reason, TerminationReason::kOptimal)
<< result.termination;
// The following code will print:
// objective value: 2.5
// value for variable x: 1
std::cout << "objective value: " << result.objective_value()
<< "\nvalue for variable x: " << result.variable_values().at(x)
<< std::endl;
}
void SolveVersion2() {
using ::operations_research::math_opt::LinearExpression;
using ::operations_research::math_opt::Model;
using ::operations_research::math_opt::SolveResult;
using ::operations_research::math_opt::SolverType;
using ::operations_research::math_opt::TerminationReason;
using ::operations_research::math_opt::Variable;
Model model("my_model");
const Variable x = model.AddBinaryVariable("x");
const Variable y = model.AddContinuousVariable(0.0, 2.5, "y");
absl::Status Main() {
math_opt::Model model("my_model");
const math_opt::Variable x = model.AddBinaryVariable("x");
const math_opt::Variable y = model.AddContinuousVariable(0.0, 2.5, "y");
// We can directly use linear combinations of variables ...
model.AddLinearConstraint(x + y <= 1.5, "c");
// ... or build them incrementally.
LinearExpression objective_expression;
math_opt::LinearExpression objective_expression;
objective_expression += 2 * x;
objective_expression += y;
model.Maximize(objective_expression);
const SolveResult result = Solve(model, SolverType::kGscip).value();
CHECK_EQ(result.termination.reason, TerminationReason::kOptimal)
<< result.termination;
// The following code will print:
// objective value: 2.5
// value for variable x: 1
std::cout << "objective value: " << result.objective_value()
<< "\nvalue for variable x: " << result.variable_values().at(x)
<< std::endl;
ASSIGN_OR_RETURN(const math_opt::SolveResult result,
Solve(model, math_opt::SolverType::kGscip));
switch (result.termination.reason) {
case math_opt::TerminationReason::kOptimal:
case math_opt::TerminationReason::kFeasible:
std::cout << "objective value: " << result.objective_value() << std::endl
<< "value for variable x: " << result.variable_values().at(x)
<< std::endl;
return absl::OkStatus();
default:
return util::InternalErrorBuilder()
<< "model failed to solve: " << result.termination;
}
}
} // namespace
int main(int argc, char** argv) {
InitGoogle(argv[0], &argc, &argv, true);
SolveVersion1();
SolveVersion2();
const absl::Status status = Main();
if (!status.ok()) {
LOG(QFATAL) << status;
}
return 0;
}

51
ortools/math_opt/samples/cocktail_hour.cc
View File

@@ -45,6 +45,7 @@
#include "ortools/base/map_util.h"
#include "ortools/base/status_macros.h"
#include "ortools/math_opt/cpp/math_opt.h"
#include "ortools/util/status_macros.h"
ABSL_FLAG(std::string, mode, "text",
"One of \"text\", \"latex\", or \"analysis\".");
@@ -251,10 +252,14 @@ absl::StatusOr<Menu> SolveForMenu(
ASSIGN_OR_RETURN(const math_opt::SolveResult result,
math_opt::Solve(model, math_opt::SolverType::kGscip, args));
// Check that the problem has an optimal solution.
QCHECK_EQ(result.termination.reason, math_opt::TerminationReason::kOptimal)
<< "Failed to find an optimal solution: " << result.termination;
switch (result.termination.reason) {
case math_opt::TerminationReason::kOptimal:
case math_opt::TerminationReason::kFeasible:
break;
default:
return util::InternalErrorBuilder()
<< "failed to find a solution: " << result.termination;
}
Menu menu;
for (const absl::string_view ingredient : kIngredients) {
if (result.variable_values().at(ingredient_vars.at(ingredient)) > 0.5) {
@@ -315,7 +320,7 @@ std::string ExportToLaTeX(const std::vector<Cocktail>& cocktails,
return absl::StrReplaceAll(absl::StrJoin(lines, "\n"), {{"#", "\\#"}});
}
void RealMain() {
absl::Status Main() {
const std::string mode = absl::GetFlag(FLAGS_mode);
CHECK(absl::flat_hash_set<std::string>({"text", "latex", "analysis"})
.contains(mode))
@@ -323,52 +328,50 @@ void RealMain() {
// We are in analysis mode.
if (mode == "analysis") {
const absl::Status status = AnalysisMode();
if (!status.ok()) {
LOG(QFATAL) << status;
}
return;
return AnalysisMode();
}
absl::StatusOr<Menu> menu =
OR_ASSIGN_OR_RETURN3(
Menu menu,
SolveForMenu(absl::GetFlag(FLAGS_num_ingredients), mode == "text",
SetFromVec(absl::GetFlag(FLAGS_existing_ingredients)),
SetFromVec(absl::GetFlag(FLAGS_unavailable_ingredients)),
SetFromVec(absl::GetFlag(FLAGS_required_cocktails)),
SetFromVec(absl::GetFlag(FLAGS_blocked_cocktails)));
if (!menu.ok()) {
LOG(QFATAL) << "Error when solving for optimal set of ingredients: "
<< menu.status();
}
SetFromVec(absl::GetFlag(FLAGS_blocked_cocktails))),
_ << "error when solving for optimal set of ingredients");
// We are in latex mode.
if (mode == "latex") {
std::cout << ExportToLaTeX(menu->cocktails) << std::endl;
return;
std::cout << ExportToLaTeX(menu.cocktails) << std::endl;
return absl::OkStatus();
}
// We are in text mode
std::cout << "Considered " << AllCocktails().size() << " cocktails and "
<< kIngredientsSize << " ingredients." << std::endl;
std::cout << "Solution has " << menu->ingredients.size()
<< " ingredients to make " << menu->cocktails.size()
<< " cocktails." << std::endl
std::cout << "Solution has " << menu.ingredients.size()
<< " ingredients to make " << menu.cocktails.size() << " cocktails."
<< std::endl
<< std::endl;
std::cout << "Ingredients:" << std::endl;
for (const std::string& ingredient : menu->ingredients) {
for (const std::string& ingredient : menu.ingredients) {
std::cout << " " << ingredient << std::endl;
}
std::cout << "Cocktails:" << std::endl;
for (const Cocktail& cocktail : menu->cocktails) {
for (const Cocktail& cocktail : menu.cocktails) {
std::cout << " " << cocktail.name << std::endl;
}
return absl::OkStatus();
}
} // namespace
int main(int argc, char** argv) {
InitGoogle(argv[0], &argc, &argv, true);
RealMain();
const absl::Status status = Main();
if (!status.ok()) {
LOG(QFATAL) << status;
}
return 0;
}

33
ortools/math_opt/samples/cutting_stock.cc
View File

@@ -137,8 +137,8 @@ absl::StatusOr<std::pair<Configuration, double>> BestConfiguration(
math_opt::Solve(model, math_opt::SolverType::kCpSat));
if (solve_result.termination.reason !=
math_opt::TerminationReason::kOptimal) {
return util::InvalidArgumentErrorBuilder()
<< "Failed to solve knapsack pricing problem: "
return util::InternalErrorBuilder()
<< "Failed to solve knapsack pricing problem to optimality: "
<< solve_result.termination;
}
Configuration config;
@@ -185,18 +185,26 @@ absl::StatusOr<CuttingStockSolution> SolveCuttingStock(
}
ASSIGN_OR_RETURN(auto solver, math_opt::IncrementalSolver::New(
model, math_opt::SolverType::kGlop));
&model, math_opt::SolverType::kGlop));
int pricing_round = 0;
while (true) {
ASSIGN_OR_RETURN(math_opt::SolveResult solve_result, solver->Solve());
if (solve_result.termination.reason !=
math_opt::TerminationReason::kOptimal) {
return util::InternalErrorBuilder()
<< "Failed to solve leader LP problem at iteration "
<< "Failed to solve leader LP problem to optimality at "
"iteration "
<< pricing_round << " termination: " << solve_result.termination;
}
// GLOP always returns a dual solution on optimal
CHECK(solve_result.has_dual_feasible_solution());
if (!solve_result.has_dual_feasible_solution()) {
// MathOpt does not require solvers to return a dual solution on optimal,
// but most LP solvers always will, see go/mathopt-solver-contracts for
// details.
return util::InternalErrorBuilder()
<< "no dual solution was returned with optimal solution at "
"iteration "
<< pricing_round;
}
std::vector<double> prices;
for (const math_opt::LinearConstraint d : demand_met) {
prices.push_back(solve_result.dual_values().at(d));
@@ -219,11 +227,14 @@ absl::StatusOr<CuttingStockSolution> SolveCuttingStock(
}
ASSIGN_OR_RETURN(const math_opt::SolveResult solve_result,
math_opt::Solve(model, math_opt::SolverType::kCpSat));
if (solve_result.termination.reason !=
math_opt::TerminationReason::kOptimal) {
return util::InternalErrorBuilder()
<< "Failed to solve final cutting stock MIP, termination: "
<< solve_result.termination;
switch (solve_result.termination.reason) {
case math_opt::TerminationReason::kOptimal:
case math_opt::TerminationReason::kFeasible:
break;
default:
return util::InternalErrorBuilder()
<< "Failed to solve final cutting stock MIP, termination: "
<< solve_result.termination;
}
CuttingStockSolution solution;
for (const auto& [config, var] : configs) {

725
ortools/math_opt/samples/facility_lp_benders.cc
View File

@@ -11,7 +11,40 @@
// See the License for the specific language governing permissions and
// limitations under the License.
// Simple linear programming example
// An advanced benders decomposition example
//
// We consider a network design problem where each location has a demand that
// must be met by its neighboring facilities, and each facility can control
// its total capacity. In this version we also require that locations cannot
// use more that a specified fraction of a facilities capacity.
//
// Problem data:
// * F: set of facilities.
// * L: set of locations.
// * E: subset of {(f,l) : f in F, l in L} that describes the network between
// facilities and locations.
// * d: demand at location (all demands are equal for simplicity).
// * c: cost per unit of capacity at a facility (all facilities are have the
// same cost for simplicity).
// * h: cost per unit transported through an edge.
// * a: fraction of a facility's capacity that can be used by each location.
//
// Decision variables:
// * z_f: capacity at facility f in F.
// * x_(f,l): flow from facility f to location l for all (f,l) in E.
//
// Formulation:
//
// min c * sum(z_f : f in F) + sum(h_e * x_e : e in E)
// s.t.
// x_(f,l) <= a * z_f for all (f,l) in E
// sum(x_(f,l) : l such that (f,l) in E) <= z_f for all f in F
// sum(x_(f,l) : f such that (f,l) in E) >= d for all l in L
// x_e >= 0 for all e in E
// z_f >= 0 for all f in F
//
// Below we solve this problem directly and using a benders decompostion
// approach.
#include <algorithm>
#include <iostream>
@@ -22,7 +55,9 @@
#include "absl/container/flat_hash_map.h"
#include "absl/flags/flag.h"
#include "absl/random/random.h"
#include "absl/random/seed_sequences.h"
#include "absl/random/uniform_int_distribution.h"
#include "absl/status/status.h"
#include "absl/status/statusor.h"
#include "absl/strings/str_format.h"
#include "absl/strings/string_view.h"
@@ -30,7 +65,9 @@
#include "absl/time/time.h"
#include "ortools/base/init_google.h"
#include "ortools/base/logging.h"
#include "ortools/base/status_macros.h"
#include "ortools/math_opt/cpp/math_opt.h"
#include "ortools/util/status_macros.h"
ABSL_FLAG(int, num_facilities, 3000, "Number of facilities.");
ABSL_FLAG(int, num_locations, 50, "Number of locations.");
@@ -41,18 +78,19 @@ ABSL_FLAG(double, facility_cost, 100, "Facility capacity cost.");
ABSL_FLAG(
double, location_fraction, 0.001,
"Fraction of a facility's capacity that can be used by each location.");
ABSL_FLAG(operations_research::math_opt::SolverType, solver_type,
operations_research::math_opt::SolverType::kGlop,
"the LP solver to use, possible values: glop, gurobi, glpk, pdlp");
namespace {
using ::operations_research::math_opt::IncrementalSolver;
using ::operations_research::math_opt::LinearConstraint;
using ::operations_research::math_opt::LinearExpression;
using ::operations_research::math_opt::Model;
using ::operations_research::math_opt::SolveArguments;
using ::operations_research::math_opt::SolveResult;
using ::operations_research::math_opt::SolverType;
using ::operations_research::math_opt::Sum;
using ::operations_research::math_opt::TerminationReason;
using ::operations_research::math_opt::Variable;
namespace math_opt = operations_research::math_opt;
constexpr double kInf = std::numeric_limits<double>::infinity();
constexpr double kZeroTol = 1.0e-3;
////////////////////////////////////////////////////////////////////////////////
// Facility location instance representation and generation
////////////////////////////////////////////////////////////////////////////////
// First element is a facility and second is a location.
using Edge = std::pair<int, int>;
@@ -91,7 +129,8 @@ class Network {
Network::Network(const int num_facilities, const int num_locations,
const double edge_probability) {
absl::BitGen bitgen;
absl::SeedSeq seq({1, 2, 3});
absl::BitGen bitgen(seq);
num_facilities_ = num_facilities;
num_locations_ = num_locations;
@@ -138,287 +177,460 @@ Network::Network(const int num_facilities, const int num_locations,
}
}
constexpr double kInf = std::numeric_limits<double>::infinity();
struct FacilityLocationInstance {
Network network;
double location_demand;
double facility_cost;
double location_fraction;
};
// We consider a network design problem where each location has a demand that
// must be met by its neighboring facilities, and each facility can control
// its total capacity. In this version we also require that locations cannot
// use more that a specified fraction of a facilities capacity.
//
// Problem data:
// * F: set of facilities.
// * L: set of locations.
// * E: subset of {(f,l) : f in F, l in L} that describes the network between
// facilities and locations.
// * d: demand at location (all demands are equal for simplicity).
// * c: cost per unit of capacity at a facility (all facilities are have the
// same cost for simplicity).
// * h: cost per unit transported through an edge.
// * a: fraction of a facility's capacity that can be used by each location.
//
// Decision variables:
// * z_f: capacity at facility f in F.
// * x_(f,l): flow from facility f to location l for all (f,l) in E.
//
// Formulation:
// min c * sum(z_f : f in F) + sum(h_e * x_e : e in E)
// s.t.
// x_(f,l) <= a * z_f for all (f,l) in E
// sum(x_(f,l) : l such that (f,l) in E) <= z_f for all f in F
// sum(x_(f,l) : f such that (f,l) in E) >= d for all l in L
// x_e >= 0 for all e in E
// z_f >= 0 for all f in F
//
void FullProblem(const Network& network, const double location_demand,
const double facility_cost, const double location_fraction) {
const int num_facilities = network.num_facilities();
const int num_locations = network.num_locations();
////////////////////////////////////////////////////////////////////////////////
// Direct solve
////////////////////////////////////////////////////////////////////////////////
Model model("Full network design problem");
model.set_minimize();
// See file level comment for problem description and formulation.
absl::Status FullProblem(const FacilityLocationInstance& instance,
const math_opt::SolverType solver_type) {
const int num_facilities = instance.network.num_facilities();
const int num_locations = instance.network.num_locations();
math_opt::Model model("Full network design problem");
// Capacity variables
std::vector<Variable> z;
std::vector<math_opt::Variable> z;
for (int j = 0; j < num_facilities; j++) {
const Variable z_j = model.AddContinuousVariable(0.0, kInf);
z.push_back(z_j);
model.set_objective_coefficient(z_j, facility_cost);
z.push_back(model.AddContinuousVariable(0.0, kInf));
}
// Flow variables
absl::flat_hash_map<Edge, Variable> x;
for (const auto& edge : network.edges()) {
const Variable x_edge = model.AddContinuousVariable(0.0, kInf);
absl::flat_hash_map<Edge, math_opt::Variable> x;
for (const auto& edge : instance.network.edges()) {
const math_opt::Variable x_edge = model.AddContinuousVariable(0.0, kInf);
x.insert({edge, x_edge});
model.set_objective_coefficient(x_edge, network.edge_cost(edge));
}
// Objective Function
math_opt::LinearExpression objective_for_edges;
for (const auto& [edge, x_edge] : x) {
objective_for_edges += instance.network.edge_cost(edge) * x_edge;
}
model.Minimize(objective_for_edges +
instance.facility_cost * math_opt::Sum(z));
// Demand constraints
for (int location = 0; location < num_locations; ++location) {
LinearExpression linear_expression;
for (const auto& edge : network.edges_incident_to_location(location)) {
linear_expression += x.at(edge);
math_opt::LinearExpression incomming_supply;
for (const auto& edge :
instance.network.edges_incident_to_location(location)) {
incomming_supply += x.at(edge);
}
model.AddLinearConstraint(linear_expression >= location_demand);
model.AddLinearConstraint(incomming_supply >= instance.location_demand);
}
// Supply constraints
for (int facility = 0; facility < num_facilities; ++facility) {
LinearExpression linear_expression;
for (const auto& edge : network.edges_incident_to_facility(facility)) {
linear_expression += x.at(edge);
math_opt::LinearExpression outgoing_supply;
for (const auto& edge :
instance.network.edges_incident_to_facility(facility)) {
outgoing_supply += x.at(edge);
}
model.AddLinearConstraint(linear_expression <= z[facility]);
model.AddLinearConstraint(outgoing_supply <= z[facility]);
}
// Arc constraints
for (int facility = 0; facility < num_facilities; ++facility) {
for (const auto& edge : network.edges_incident_to_facility(facility)) {
model.AddLinearConstraint(x.at(edge) <= location_fraction * z[facility]);
for (const auto& edge :
instance.network.edges_incident_to_facility(facility)) {
model.AddLinearConstraint(x.at(edge) <=
instance.location_fraction * z[facility]);
}
}
const SolveResult result = Solve(model, SolverType::kGurobi).value();
QCHECK_EQ(result.termination.reason, TerminationReason::kOptimal)
<< "Failed to find an optimal solution: " << result.termination;
ASSIGN_OR_RETURN(const math_opt::SolveResult result,
math_opt::Solve(model, solver_type));
if (result.termination.reason != math_opt::TerminationReason::kOptimal) {
return util::InternalErrorBuilder()
<< "failed to find an optimal solution: " << result.termination;
}
std::cout << "Full problem optimal objective: "
<< absl::StrFormat("%.9f", result.objective_value()) << std::endl;
return absl::OkStatus();
}
void Benders(const Network network, const double location_demand,
const double facility_cost, const double location_fraction,
const double target_precission,
const int maximum_iterations = 30000) {
////////////////////////////////////////////////////////////////////////////////
// Benders solver
////////////////////////////////////////////////////////////////////////////////
// Setup first stage model:
//
// min c * sum(z_f : f in F) + w
// s.t.
// z_f >= 0 for all f in F
// sum(fcut_f^i z_f) + fcut_const^i <= 0 for i = 1,...
// sum(ocut_f^j z_f) + ocut_const^j <= w for j = 1,...
struct FirstStageProblem {
math_opt::Model model;
std::vector<math_opt::Variable> z;
math_opt::Variable w;
FirstStageProblem(const Network& network, const double facility_cost);
};
FirstStageProblem::FirstStageProblem(const Network& network,
const double facility_cost)
: model("First stage problem"), w(model.AddContinuousVariable(0.0, kInf)) {
const int num_facilities = network.num_facilities();
const int num_locations = network.num_locations();
// Setup first stage model.
//
// min c * sum(z_f : f in F) + w
// s.t.
// z_f >= 0 for all f in F
// sum(fcut_f^i z_f) + fcut_const^i <= 0 for i = 1,...
// sum(ocut_f^j z_f) + ocut_const^j <= w for j = 1,...
Model first_stage_model("First stage problem");
std::vector<Variable> z;
// Capacity variables
for (int j = 0; j < num_facilities; j++) {
z.push_back(first_stage_model.AddContinuousVariable(0.0, kInf));
z.push_back(model.AddContinuousVariable(0.0, kInf));
}
const Variable w = first_stage_model.AddContinuousVariable(0.0, kInf);
// First stage objective
model.Minimize(w + facility_cost * Sum(z));
}
first_stage_model.Minimize(facility_cost * Sum(z) + w);
// Represents a cut if the form:
//
// z_coefficients^T z + constant <= w_coefficient * w
//
// This will be a feasibility cut if w_coefficient = 0.0 and an optimality
// cut if w_coefficient = 1.
struct Cut {
std::vector<double> z_coefficients;
double constant;
double w_coefficient;
};
// Setup second stage model.
// min sum(h_e * x_e : e in E)
// s.t.
// x_(f,l) <= a * zz_f for all (f,l) in E
// sum(x_(f,l) : l such that (f,l) in E) <= zz_f for all f in F
// sum(x_(f,l) : f such that (f,l) in E) >= d for all l in L
// x_e >= 0 for all e in E
//
// where zz_f are fixed values for z_f from the first stage model.
Model second_stage_model("Second stage model");
second_stage_model.set_minimize();
absl::flat_hash_map<Edge, Variable> x;
for (const auto& edge : network.edges()) {
const Variable x_edge = second_stage_model.AddContinuousVariable(0.0, kInf);
x.insert({edge, x_edge});
second_stage_model.set_objective_coefficient(x_edge,
network.edge_cost(edge));
// Solves the second stage model:
//
// min sum(h_e * x_e : e in E)
// s.t.
// x_(f,l) <= a * zz_f for all (f,l) in E
// sum(x_(f,l) : l such that (f,l) in E) <= zz_f for all f in F
// sum(x_(f,l) : f such that (f,l) in E) >= d for all l in L
// x_e >= 0 for all e in E
//
// where zz_f are fixed values for z_f from the first stage model, and generates
// an infeasibility or optimality cut as needed.
class SecondStageSolver {
public:
static absl::StatusOr<std::unique_ptr<SecondStageSolver>> New(
FacilityLocationInstance instance, math_opt::SolverType solver_type);
absl::StatusOr<std::pair<double, Cut>> Solve(
const std::vector<double>& z_values, double w_value,
double fist_stage_objective);
private:
SecondStageSolver(FacilityLocationInstance instance,
math_opt::SolveParameters second_stage_params);
absl::StatusOr<Cut> OptimalityCut(
const math_opt::SolveResult& second_stage_result);
absl::StatusOr<Cut> FeasibilityCut(
const math_opt::SolveResult& second_stage_result);
math_opt::Model second_stage_model_;
const Network network_;
const double location_fraction_;
math_opt::SolveParameters second_stage_params_;
absl::flat_hash_map<Edge, math_opt::Variable> x_;
std::vector<math_opt::LinearConstraint> supply_constraints_;
std::vector<math_opt::LinearConstraint> demand_constraints_;
std::unique_ptr<math_opt::IncrementalSolver> solver_;
};
absl::StatusOr<math_opt::SolveParameters> EnsureDualRaySolveParameters(
const math_opt::SolverType solver_type) {
math_opt::SolveParameters parameters;
switch (solver_type) {
case math_opt::SolverType::kGurobi:
parameters.gurobi.param_values["InfUnbdInfo"] = "1";
break;
case math_opt::SolverType::kGlop:
parameters.presolve = math_opt::Emphasis::kOff;
parameters.scaling = math_opt::Emphasis::kOff;
parameters.lp_algorithm = math_opt::LPAlgorithm::kDualSimplex;
break;
case math_opt::SolverType::kGlpk:
parameters.presolve = math_opt::Emphasis::kOff;
parameters.lp_algorithm = math_opt::LPAlgorithm::kDualSimplex;
break;
default:
return util::InternalErrorBuilder()
<< "unsupported solver: " << solver_type;
}
std::vector<LinearConstraint> demand_constraints;
return parameters;
}
absl::StatusOr<std::unique_ptr<SecondStageSolver>> SecondStageSolver::New(
FacilityLocationInstance instance, const math_opt::SolverType solver_type) {
// Set solver arguments to ensure a dual ray is returned.
ASSIGN_OR_RETURN(math_opt::SolveParameters parameters,
EnsureDualRaySolveParameters(solver_type));
std::unique_ptr<SecondStageSolver> second_stage_solver =
absl::WrapUnique<SecondStageSolver>(
new SecondStageSolver(std::move(instance), parameters));
ASSIGN_OR_RETURN(std::unique_ptr<math_opt::IncrementalSolver> solver,
math_opt::IncrementalSolver::New(
&second_stage_solver->second_stage_model_, solver_type));
second_stage_solver->solver_ = std::move(solver);
return std::move(second_stage_solver);
}
SecondStageSolver::SecondStageSolver(
FacilityLocationInstance instance,
math_opt::SolveParameters second_stage_params)
: second_stage_model_("Second stage model"),
network_(std::move(instance.network)),
location_fraction_(instance.location_fraction),
second_stage_params_(second_stage_params) {
const int num_facilities = network_.num_facilities();
const int num_locations = network_.num_locations();
// Flow variables
for (const auto& edge : network_.edges()) {
const math_opt::Variable x_edge =
second_stage_model_.AddContinuousVariable(0.0, kInf);
x_.insert({edge, x_edge});
}
// Objective Function
math_opt::LinearExpression objective_for_edges;
for (const auto& [edge, x_edge] : x_) {
objective_for_edges += network_.edge_cost(edge) * x_edge;
}
second_stage_model_.Minimize(objective_for_edges);
// Demand constraints
for (int location = 0; location < num_locations; ++location) {
LinearExpression linear_expression;
for (const auto& edge : network.edges_incident_to_location(location)) {
linear_expression += x.at(edge);
math_opt::LinearExpression incomming_supply;
for (const auto& edge : network_.edges_incident_to_location(location)) {
incomming_supply += x_.at(edge);
}
demand_constraints.push_back(second_stage_model.AddLinearConstraint(
linear_expression >= location_demand));
demand_constraints_.push_back(second_stage_model_.AddLinearConstraint(
incomming_supply >= instance.location_demand));
}
std::vector<LinearConstraint> supply_constraints;
// Supply constraints
for (int facility = 0; facility < num_facilities; ++facility) {
LinearExpression linear_expression;
for (const auto& edge : network.edges_incident_to_facility(facility)) {
linear_expression += x.at(edge);
math_opt::LinearExpression outgoing_supply;
for (const auto& edge : network_.edges_incident_to_facility(facility)) {
outgoing_supply += x_.at(edge);
}
supply_constraints.push_back(
second_stage_model.AddLinearConstraint(linear_expression <= kInf));
// Set supply constraint with trivial upper bound to be updated with first
// stage information.
supply_constraints_.push_back(
second_stage_model_.AddLinearConstraint(outgoing_supply <= kInf));
}
}
absl::StatusOr<std::pair<double, Cut>> SecondStageSolver::Solve(
const std::vector<double>& z_values, const double w_value,
const double fist_stage_objective) {
const int num_facilities = network_.num_facilities();
// Update second stage with first stage solution.
for (int facility = 0; facility < num_facilities; ++facility) {
if (z_values[facility] < -kZeroTol) {
return util::InternalErrorBuilder()
<< "negative z_value in first stage: " << z_values[facility]
<< " for facility " << facility;
}
// Make sure variable bounds are valid (lb <= ub).
const double capacity_value = std::max(z_values[facility], 0.0);
for (const auto& edge : network_.edges_incident_to_facility(facility)) {
second_stage_model_.set_upper_bound(x_.at(edge),
location_fraction_ * capacity_value);
}
second_stage_model_.set_upper_bound(supply_constraints_[facility],
capacity_value);
}
// Solve and process second stage.
ASSIGN_OR_RETURN(const math_opt::SolveResult second_stage_result,
solver_->Solve(math_opt::SolveArguments{
.parameters = second_stage_params_}));
switch (second_stage_result.termination.reason) {
case math_opt::TerminationReason::kInfeasible:
// If the second stage problem is infeasible we can construct a
// feasibility cut from a returned dual ray.
{
OR_ASSIGN_OR_RETURN3(const Cut feasibility_cut,
FeasibilityCut(second_stage_result),
_ << "on infeasible for second stage solver");
return std::make_pair(kInf, feasibility_cut);
}
case math_opt::TerminationReason::kOptimal:
// If the second stage problem is optimal we can construct an optimality
// cut from a returned dual optimal solution. We can also update the upper
// bound.
{
// Upper bound is obtained by switching predicted second stage objective
// value w with the true second stage objective value.
const double upper_bound = fist_stage_objective - w_value +
second_stage_result.objective_value();
OR_ASSIGN_OR_RETURN3(const Cut optimality_cut,
OptimalityCut(second_stage_result),
_ << "on optimal for second stage solver");
return std::make_pair(upper_bound, optimality_cut);
}
default:
return util::InternalErrorBuilder()
<< "second stage was not solved to optimality or infeasibility: "
<< second_stage_result.termination;
}
}
// If the second stage problem is infeasible we get a dual ray (r, y) such that
//
// sum(r_(f,l)*a*zz_f : (f,l) in E, r_(f,l) < 0)
// + sum(y_f*zz_f : f in F, y_f < 0)
// + sum(y_l*d : l in L, y_l > 0) > 0.
//
// Then we get the feasibility cut (go/math_opt-advanced-dual-use#benders)
//
// sum(fcut_f*z_f) + fcut_const <= 0,
//
// where
//
// fcut_f = sum(r_(f,l)*a : (f,l) in E, r_(f,l) < 0)
// + min{y_f, 0}
// fcut_const = sum*(y_l*d : l in L, y_l > 0)
absl::StatusOr<Cut> SecondStageSolver::FeasibilityCut(
const math_opt::SolveResult& second_stage_result) {
const int num_facilities = network_.num_facilities();
Cut result;
if (!second_stage_result.has_dual_ray()) {
// MathOpt does not require solvers to return a dual solution on optimal,
// but most LP solvers always will, see go/mathopt-solver-contracts for
// details.
return util::InternalErrorBuilder()
<< "no dual ray available for feasibility cut";
}
SolveArguments second_stage_args;
second_stage_args.parameters.gurobi.param_values["InfUnbdInfo"] = "1";
for (int facility = 0; facility < num_facilities; ++facility) {
double coefficient = 0.0;
for (const auto& edge : network_.edges_incident_to_facility(facility)) {
const double reduced_cost =
second_stage_result.ray_reduced_costs().at(x_.at(edge));
coefficient += location_fraction_ * std::min(reduced_cost, 0.0);
}
const double dual_value =
second_stage_result.ray_dual_values().at(supply_constraints_[facility]);
coefficient += std::min(dual_value, 0.0);
result.z_coefficients.push_back(coefficient);
}
result.constant = 0.0;
for (const auto& constraint : demand_constraints_) {
const double dual_value =
second_stage_result.ray_dual_values().at(constraint);
result.constant += std::max(dual_value, 0.0);
}
result.w_coefficient = 0.0;
return result;
}
// If the second stage problem is optimal we get a dual solution (r, y) such
// that the optimal objective value is equal to
//
// sum(r_(f,l)*a*zz_f : (f,l) in E, r_(f,l) < 0)
// + sum(y_f*zz_f : f in F, y_f < 0)
// + sum*(y_l*d : l in L, y_l > 0) > 0.
//
// Then we get the optimality cut (go/math_opt-advanced-dual-use#benders)
//
// sum(ocut_f*z_f) + ocut_const <= w,
//
// where
//
// ocut_f = sum(r_(f,l)*a : (f,l) in E, r_(f,l) < 0)
// + min{y_f, 0}
// ocut_const = sum*(y_l*d : l in L, y_l > 0)
absl::StatusOr<Cut> SecondStageSolver::OptimalityCut(
const math_opt::SolveResult& second_stage_result) {
const int num_facilities = network_.num_facilities();
Cut result;
if (!second_stage_result.has_dual_feasible_solution()) {
return util::InternalErrorBuilder()
<< "no dual solution available for optimality cut";
}
for (int facility = 0; facility < num_facilities; ++facility) {
double coefficient = 0.0;
for (const auto& edge : network_.edges_incident_to_facility(facility)) {
const double reduced_cost =
second_stage_result.reduced_costs().at(x_.at(edge));
coefficient += location_fraction_ * std::min(reduced_cost, 0.0);
}
double dual_value =
second_stage_result.dual_values().at(supply_constraints_[facility]);
coefficient += std::min(dual_value, 0.0);
result.z_coefficients.push_back(coefficient);
}
result.constant = 0.0;
for (const auto& constraint : demand_constraints_) {
double dual_value = second_stage_result.dual_values().at(constraint);
result.constant += std::max(dual_value, 0.0);
}
result.w_coefficient = 1.0;
return result;
}
absl::Status Benders(const FacilityLocationInstance& instance,
const double target_precission,
const math_opt::SolverType solver_type,
const int maximum_iterations = 30000) {
const int num_facilities = instance.network.num_facilities();
// Setup first stage model and solver.
FirstStageProblem first_stage(instance.network, instance.facility_cost);
ASSIGN_OR_RETURN(
const std::unique_ptr<math_opt::IncrementalSolver> first_stage_solver,
math_opt::IncrementalSolver::New(&first_stage.model, solver_type));
// Setup second stage solver.
ASSIGN_OR_RETURN(std::unique_ptr<SecondStageSolver> second_stage_solver,
SecondStageSolver::New(instance, solver_type));
// Start Benders
int iteration = 0;
double best_upper_bound = kInf;
const std::unique_ptr<IncrementalSolver> first_stage_solver =
IncrementalSolver::New(first_stage_model, SolverType::kGurobi).value();
const std::unique_ptr<IncrementalSolver> second_stage_solver =
IncrementalSolver::New(second_stage_model, SolverType::kGurobi).value();
std::vector<double> z_values;
z_values.resize(num_facilities);
while (true) {
LOG(INFO) << "Iteration: " << iteration;
// Solve and process first stage.
const SolveResult first_stage_result = first_stage_solver->Solve().value();
QCHECK_EQ(first_stage_result.termination.reason,
TerminationReason::kOptimal)
<< first_stage_result.termination;
ASSIGN_OR_RETURN(const math_opt::SolveResult first_stage_result,
first_stage_solver->Solve());
if (first_stage_result.termination.reason !=
math_opt::TerminationReason::kOptimal) {
return util::InternalErrorBuilder()
<< "could not solve first stage problem to optimality:"
<< first_stage_result.termination;
}
for (int j = 0; j < num_facilities; j++) {
z_values[j] = first_stage_result.variable_values().at(first_stage.z[j]);
}
const double lower_bound = first_stage_result.objective_value();
LOG(INFO) << "LB = " << lower_bound;
// Setup second stage.
for (int facility = 0; facility < num_facilities; ++facility) {
const double capacity_value =
first_stage_result.variable_values().at(z[facility]);
for (const auto& edge : network.edges_incident_to_facility(facility)) {
second_stage_model.set_upper_bound(x.at(edge),
location_fraction * capacity_value);
}
second_stage_model.set_upper_bound(supply_constraints[facility],
capacity_value);
}
// Solve and process second stage.
const SolveResult second_stage_result =
second_stage_solver->Solve(second_stage_args).value();
if (second_stage_result.termination.reason ==
TerminationReason::kInfeasible) {
// If the second stage problem is infeasible we will get a dual ray
// (r, y) such that
//
// sum(r_(f,l)*a*zz_f : (f,l) in E, r_(f,l) < 0)
// + sum(y_f*zz_f : f in F, y_f < 0)
// + sum(y_l*d : l in L, y_l > 0) > 0.
//
// Then we get the feasibility cut (go/mathopt-advanced-dual-use#benders)
//
// sum(fcut_f*z_f) + fcut_const <= 0,
//
// where
//
// fcut_f = sum(r_(f,l)*a : (f,l) in E, r_(f,l) < 0)
// + min{y_f, 0}
// fcut_const = sum*(y_l*d : l in L, y_l > 0)
LOG(INFO) << "Adding feasibility cut...";
LinearExpression feasibility_cut_expression;
for (int facility = 0; facility < num_facilities; ++facility) {
double coefficient = 0.0;
for (const auto& edge : network.edges_incident_to_facility(facility)) {
const double reduced_cost =
second_stage_result.ray_reduced_costs().at(x.at(edge));
if (reduced_cost < 0) {
coefficient += location_fraction * reduced_cost;
}
}
double dual_value = second_stage_result.ray_dual_values().at(
supply_constraints[facility]);
coefficient += std::min(dual_value, 0.0);
feasibility_cut_expression += coefficient * z[facility];
}
double constant = 0.0;
for (const auto& constraint : demand_constraints) {
double dual_value =
second_stage_result.ray_dual_values().at(constraint);
if (dual_value > 0) {
constant += dual_value;
}
}
first_stage_model.AddLinearConstraint(
feasibility_cut_expression + constant <= 0.0);
} else {
// If the second stage problem is optimal we will get a dual solution
// (r, y) such that the optimal objective value is equal to
//
// sum(r_(f,l)*a*zz_f : (f,l) in E, r_(f,l) < 0)
// + sum(y_f*zz_f : f in F, y_f < 0)
// + sum*(y_l*d : l in L, y_l > 0) > 0.
//
// Then we get the optimality cut (go/mathopt-advanced-dual-use#benders)
//
// sum(ocut_f*z_f) + ocut_const <= w,
//
// where
//
// ocut_f = sum(r_(f,l)*a : (f,l) in E, r_(f,l) < 0)
// + min{y_f, 0}
// ocut_const = sum*(y_l*d : l in L, y_l > 0)
QCHECK_EQ(second_stage_result.termination.reason,
TerminationReason::kOptimal)
<< second_stage_result.termination;
LOG(INFO) << "Adding optimality cut...";
LinearExpression optimality_cut_expression;
double upper_bound = 0.0;
for (int facility = 0; facility < num_facilities; ++facility) {
double coefficient = 0.0;
for (const auto& edge : network.edges_incident_to_facility(facility)) {
upper_bound += network.edge_cost(edge) *
second_stage_result.variable_values().at(x.at(edge));
const double reduced_cost =
second_stage_result.reduced_costs().at(x.at(edge));
if (reduced_cost < 0) {
coefficient += location_fraction * reduced_cost;
}
}
double dual_value =
second_stage_result.dual_values().at(supply_constraints[facility]);
coefficient += std::min(dual_value, 0.0);
optimality_cut_expression += coefficient * z[facility];
}
double constant = 0.0;
for (const auto& constraint : demand_constraints) {
double dual_value = second_stage_result.dual_values().at(constraint);
if (dual_value > 0) {
constant += dual_value;
}
}
// Add first stage contribution to upper bound = facility_cost * Sum(z)
upper_bound += first_stage_result.objective_value() -
first_stage_result.variable_values().at(w);
best_upper_bound = std::min(best_upper_bound, upper_bound);
first_stage_model.AddLinearConstraint(
optimality_cut_expression + constant <= w);
ASSIGN_OR_RETURN(
(auto [upper_bound, cut]),
second_stage_solver->Solve(
z_values, first_stage_result.variable_values().at(first_stage.w),
first_stage_result.objective_value()));
math_opt::LinearExpression cut_expression;
for (int j = 0; j < num_facilities; j++) {
cut_expression += cut.z_coefficients[j] * first_stage.z[j];
}
cut_expression += cut.constant;
first_stage.model.AddLinearConstraint(cut_expression <=
cut.w_coefficient * first_stage.w);
best_upper_bound = std::min(upper_bound, best_upper_bound);
LOG(INFO) << "UB = " << best_upper_bound;
++iteration;
if (best_upper_bound - lower_bound < target_precission) {
@@ -433,25 +645,34 @@ void Benders(const Network network, const double location_demand,
break;
}
}
return absl::OkStatus();
}
absl::Status Main() {
const FacilityLocationInstance instance{
.network = Network(absl::GetFlag(FLAGS_num_facilities),
absl::GetFlag(FLAGS_num_locations),
absl::GetFlag(FLAGS_edge_probability)),
.location_demand = absl::GetFlag(FLAGS_location_demand),
.facility_cost = absl::GetFlag(FLAGS_facility_cost),
.location_fraction = absl::GetFlag(FLAGS_location_fraction)};
absl::Time start = absl::Now();
RETURN_IF_ERROR(FullProblem(instance, absl::GetFlag(FLAGS_solver_type)))
<< "full solve failed";
std::cout << "Full solve time : " << absl::Now() - start << std::endl;
start = absl::Now();
RETURN_IF_ERROR(Benders(instance, absl::GetFlag(FLAGS_benders_precission),
absl::GetFlag(FLAGS_solver_type)))
<< "Benders solve failed";
std::cout << "Benders solve time : " << absl::Now() - start << std::endl;
return absl::OkStatus();
}
} // namespace
int main(int argc, char** argv) {
InitGoogle(argv[0], &argc, &argv, true);
const Network network(absl::GetFlag(FLAGS_num_facilities),
absl::GetFlag(FLAGS_num_locations),
absl::GetFlag(FLAGS_edge_probability));
absl::Time start = absl::Now();
FullProblem(network, absl::GetFlag(FLAGS_location_demand),
absl::GetFlag(FLAGS_facility_cost),
absl::GetFlag(FLAGS_location_fraction));
std::cout << "Full solve time : " << absl::Now() - start << std::endl;
start = absl::Now();
Benders(network, absl::GetFlag(FLAGS_location_demand),
absl::GetFlag(FLAGS_facility_cost),
absl::GetFlag(FLAGS_location_fraction),
absl::GetFlag(FLAGS_benders_precission));
std::cout << "Benders solve time : " << absl::Now() - start << std::endl;
const absl::Status status = Main();
if (!status.ok()) {
LOG(QFATAL) << status;
}
return 0;
}

53
ortools/math_opt/samples/integer_programming.cc
View File

@@ -20,15 +20,13 @@
#include "absl/time/time.h"
#include "ortools/base/init_google.h"
#include "ortools/base/logging.h"
#include "ortools/base/status_builder.h"
#include "ortools/base/status_macros.h"
#include "ortools/math_opt/cpp/math_opt.h"
namespace {
using ::operations_research::math_opt::Model;
using ::operations_research::math_opt::SolveResult;
using ::operations_research::math_opt::SolverType;
using ::operations_research::math_opt::TerminationReason;
using ::operations_research::math_opt::Variable;
using ::operations_research::math_opt::VariableMap;
namespace math_opt = ::operations_research::math_opt;
constexpr double kInf = std::numeric_limits<double>::infinity();
@@ -39,12 +37,12 @@ constexpr double kInf = std::numeric_limits<double>::infinity();
// x in {0.0, 1.0, 2.0, ...,
// y in {0.0, 1.0, 2.0, ...,
//
void SolveSimpleMIP() {
Model model("Integer programming example");
absl::Status Main() {
math_opt::Model model("Integer programming example");
// Variables
const Variable x = model.AddIntegerVariable(0.0, kInf, "x");
const Variable y = model.AddIntegerVariable(0.0, kInf, "y");
const math_opt::Variable x = model.AddIntegerVariable(0.0, kInf, "x");
const math_opt::Variable y = model.AddIntegerVariable(0.0, kInf, "y");
// Constraints
model.AddLinearConstraint(x + 7 * y <= 17.5, "c1");
@@ -53,28 +51,29 @@ void SolveSimpleMIP() {
// Objective
model.Maximize(x + 10 * y);
std::cout << "Num variables: " << model.num_variables() << std::endl;
std::cout << "Num constraints: " << model.num_linear_constraints()
<< std::endl;
ASSIGN_OR_RETURN(const math_opt::SolveResult result,
Solve(model, math_opt::SolverType::kGscip));
const SolveResult result = Solve(model, SolverType::kGscip).value();
// Check that the problem has an optimal solution.
QCHECK_EQ(result.termination.reason, TerminationReason::kOptimal)
<< "Failed to find an optimal solution: " << result.termination;
std::cout << "Problem solved in " << result.solve_time() << std::endl;
std::cout << "Objective value: " << result.objective_value() << std::endl;
const double x_val = result.variable_values().at(x);
const double y_val = result.variable_values().at(y);
std::cout << "Variable values: [x=" << x_val << ", y=" << y_val << "]"
<< std::endl;
switch (result.termination.reason) {
case math_opt::TerminationReason::kOptimal:
case math_opt::TerminationReason::kFeasible:
std::cout << "Problem solved in " << result.solve_time() << std::endl;
std::cout << "objective value: " << result.objective_value() << std::endl;
std::cout << "Variable values: [x=" << result.variable_values().at(x)
<< ", y=" << result.variable_values().at(y) << "]" << std::endl;
return absl::OkStatus();
default:
return util::InternalErrorBuilder()
<< "model failed to solve: " << result.termination;
}
}
} // namespace
int main(int argc, char** argv) {
InitGoogle(argv[0], &argc, &argv, true);
SolveSimpleMIP();
const absl::Status status = Main();
if (!status.ok()) {
LOG(QFATAL) << status;
}
return 0;
}

192
ortools/math_opt/samples/lagrangian_relaxation.cc
View File

@@ -87,6 +87,7 @@
#include "absl/flags/flag.h"
#include "absl/memory/memory.h"
#include "absl/status/status.h"
#include "absl/status/statusor.h"
#include "absl/strings/str_format.h"
#include "absl/strings/string_view.h"
@@ -94,6 +95,8 @@
#include "ortools/base/init_google.h"
#include "ortools/base/logging.h"
#include "ortools/base/mathutil.h"
#include "ortools/base/status_builder.h"
#include "ortools/base/status_macros.h"
#include "ortools/math_opt/cpp/math_opt.h"
ABSL_FLAG(double, step_size, 0.95,
@@ -113,14 +116,8 @@ constexpr double kZeroTol = 1.0e-8;
namespace {
using ::operations_research::MathUtil;
using ::operations_research::math_opt::LinearExpression;
using ::operations_research::math_opt::Model;
using ::operations_research::math_opt::SolveArguments;
using ::operations_research::math_opt::SolveResult;
using ::operations_research::math_opt::SolverType;
using ::operations_research::math_opt::TerminationReason;
using ::operations_research::math_opt::Variable;
using ::operations_research::math_opt::VariableMap;
namespace math_opt = ::operations_research::math_opt;
struct Arc {
int i;
@@ -138,20 +135,20 @@ struct Graph {
};
struct FlowModel {
FlowModel() : model(std::make_unique<Model>("LagrangianProblem")) {}
std::unique_ptr<Model> model;
LinearExpression cost;
LinearExpression resource_1;
LinearExpression resource_2;
std::vector<Variable> flow_vars;
FlowModel() : model(std::make_unique<math_opt::Model>("LagrangianProblem")) {}
std::unique_ptr<math_opt::Model> model;
math_opt::LinearExpression cost;
math_opt::LinearExpression resource_1;
math_opt::LinearExpression resource_2;
std::vector<math_opt::Variable> flow_vars;
};
// Populates `model` with variables and constraints of a shortest path problem.
FlowModel CreateShortestPathModel(const Graph graph) {
FlowModel flow_model;
Model& model = *flow_model.model;
math_opt::Model& model = *flow_model.model;
for (const Arc& arc : graph.arcs) {
Variable var = model.AddContinuousVariable(
math_opt::Variable var = model.AddContinuousVariable(
/*lower_bound=*/0, /*upper_bound=*/1,
/*name=*/absl::StrFormat("x_%d_%d", arc.i, arc.j));
flow_model.cost += arc.cost * var;
@@ -161,8 +158,8 @@ FlowModel CreateShortestPathModel(const Graph graph) {
}
// Flow balance constraints
std::vector<LinearExpression> out_flow(graph.num_nodes);
std::vector<LinearExpression> in_flow(graph.num_nodes);
std::vector<math_opt::LinearExpression> out_flow(graph.num_nodes);
std::vector<math_opt::LinearExpression> in_flow(graph.num_nodes);
for (int arc_index = 0; arc_index < graph.arcs.size(); ++arc_index) {
out_flow[graph.arcs[arc_index].i] += flow_model.flow_vars[arc_index];
in_flow[graph.arcs[arc_index].j] += flow_model.flow_vars[arc_index];
@@ -206,12 +203,13 @@ Graph CreateSampleNetwork() {
}
// Solves the constrained shortest path as an MIP.
FlowModel SolveMip(const Graph graph, const double max_resource_1,
const double max_resource_2) {
absl::StatusOr<FlowModel> SolveMip(const Graph graph,
const double max_resource_1,
const double max_resource_2) {
FlowModel flow_model;
Model& model = *flow_model.model;
math_opt::Model& model = *flow_model.model;
for (const Arc& arc : graph.arcs) {
Variable var = model.AddBinaryVariable(
math_opt::Variable var = model.AddBinaryVariable(
/*name=*/absl::StrFormat("x_%d_%d", arc.i, arc.j));
flow_model.cost += arc.cost * var;
flow_model.resource_1 += +arc.resource_1 * var;
@@ -220,8 +218,8 @@ FlowModel SolveMip(const Graph graph, const double max_resource_1,
}
// Flow balance constraints
std::vector<LinearExpression> out_flow(graph.num_nodes);
std::vector<LinearExpression> in_flow(graph.num_nodes);
std::vector<math_opt::LinearExpression> out_flow(graph.num_nodes);
std::vector<math_opt::LinearExpression> in_flow(graph.num_nodes);
for (int arc_index = 0; arc_index < graph.arcs.size(); ++arc_index) {
out_flow[graph.arcs[arc_index].i] += flow_model.flow_vars[arc_index];
in_flow[graph.arcs[arc_index].j] += flow_model.flow_vars[arc_index];
@@ -239,51 +237,71 @@ FlowModel SolveMip(const Graph graph, const double max_resource_1,
model.AddLinearConstraint(flow_model.resource_2 <= max_resource_2,
"resource_ctr_2");
model.Minimize(flow_model.cost);
const SolveResult result = Solve(model, SolverType::kGscip).value();
const VariableMap<double>& variable_values = result.variable_values();
std::cout << "MIP Solution with 2 side constraints" << std::endl;
std::cout << absl::StrFormat("MIP objective value: %6.3f",
result.objective_value())
<< std::endl;
std::cout << "Resource 1: " << flow_model.resource_1.Evaluate(variable_values)
<< std::endl;
std::cout << "Resource 2: " << flow_model.resource_2.Evaluate(variable_values)
<< std::endl;
std::cout << "========================================" << std::endl;
return flow_model;
ASSIGN_OR_RETURN(const math_opt::SolveResult result,
Solve(model, math_opt::SolverType::kGscip));
switch (result.termination.reason) {
case math_opt::TerminationReason::kOptimal:
case math_opt::TerminationReason::kFeasible:
std::cout << "MIP Solution with 2 side constraints" << std::endl;
std::cout << absl::StrFormat("MIP objective value: %6.3f",
result.objective_value())
<< std::endl;
std::cout << "Resource 1: "
<< flow_model.resource_1.Evaluate(result.variable_values())
<< std::endl;
std::cout << "Resource 2: "
<< flow_model.resource_2.Evaluate(result.variable_values())
<< std::endl;
std::cout << "========================================" << std::endl;
return flow_model;
default:
return util::InternalErrorBuilder()
<< "model failed to solve: " << result.termination;
}
}
// Solves the linear relaxation of a constrained shortest path problem
// formulated as an MIP.
void SolveLinearRelaxation(FlowModel& flow_model, const Graph& graph,
const double max_resource_1,
const double max_resource_2) {
Model& model = *flow_model.model;
const SolveResult result = Solve(model, SolverType::kGscip).value();
const VariableMap<double>& variable_values = result.variable_values();
std::cout << "LP relaxation with 2 side constraints" << std::endl;
std::cout << absl::StrFormat("LP objective value: %6.3f",
result.objective_value())
<< std::endl;
std::cout << "Resource 1: " << flow_model.resource_1.Evaluate(variable_values)
<< std::endl;
std::cout << "Resource 2: " << flow_model.resource_2.Evaluate(variable_values)
<< std::endl;
std::cout << "========================================" << std::endl;
absl::Status SolveLinearRelaxation(FlowModel& flow_model, const Graph& graph,
const double max_resource_1,
const double max_resource_2) {
math_opt::Model& model = *flow_model.model;
ASSIGN_OR_RETURN(const math_opt::SolveResult result,
Solve(model, math_opt::SolverType::kGscip));
switch (result.termination.reason) {
case math_opt::TerminationReason::kOptimal:
case math_opt::TerminationReason::kFeasible:
std::cout << "LP relaxation with 2 side constraints" << std::endl;
std::cout << absl::StrFormat("LP objective value: %6.3f",
result.objective_value())
<< std::endl;
std::cout << "Resource 1: "
<< flow_model.resource_1.Evaluate(result.variable_values())
<< std::endl;
std::cout << "Resource 2: "
<< flow_model.resource_2.Evaluate(result.variable_values())
<< std::endl;
std::cout << "========================================" << std::endl;
return absl::OkStatus();
default:
return util::InternalErrorBuilder()
<< "model failed to solve: " << result.termination;
}
}
void SolveLagrangianRelaxation(const Graph graph, const double max_resource_1,
const double max_resource_2) {
absl::Status SolveLagrangianRelaxation(const Graph graph,
const double max_resource_1,
const double max_resource_2) {
// Model, variables, and linear expressions.
FlowModel flow_model = CreateShortestPathModel(graph);
Model& model = *flow_model.model;
LinearExpression& cost = flow_model.cost;
LinearExpression& resource_1 = flow_model.resource_1;
LinearExpression& resource_2 = flow_model.resource_2;
math_opt::Model& model = *flow_model.model;
math_opt::LinearExpression& cost = flow_model.cost;
math_opt::LinearExpression& resource_1 = flow_model.resource_1;
math_opt::LinearExpression& resource_2 = flow_model.resource_2;
// Dualized constraints and dual variable iterates.
std::vector<double> mu;
std::vector<LinearExpression> grad_mu;
std::vector<math_opt::LinearExpression> grad_mu;
const bool dualized_resource_1 = absl::GetFlag(FLAGS_dualize_resource_1);
const bool dualized_resource_2 = absl::GetFlag(FLAGS_dualize_resource_2);
const bool lagrangian_output = absl::GetFlag(FLAGS_lagrangian_output);
@@ -298,14 +316,14 @@ void SolveLagrangianRelaxation(const Graph graph, const double max_resource_1,
mu.push_back(initial_dual_value);
grad_mu.push_back(max_resource_1 - resource_1);
model.AddLinearConstraint(resource_2 <= max_resource_2);
for (Variable& var : flow_model.flow_vars) {
for (math_opt::Variable& var : flow_model.flow_vars) {
model.set_integer(var);
}
} else if (!dualized_resource_1 && dualized_resource_2) {
mu.push_back(initial_dual_value);
grad_mu.push_back(max_resource_2 - resource_2);
model.AddLinearConstraint(resource_1 <= max_resource_1);
for (Variable& var : flow_model.flow_vars) {
for (math_opt::Variable& var : flow_model.flow_vars) {
model.set_integer(var);
}
} else {
@@ -339,16 +357,25 @@ void SolveLagrangianRelaxation(const Graph graph, const double max_resource_1,
}
while (!termination) {
LinearExpression lagrangian_function;
math_opt::LinearExpression lagrangian_function;
lagrangian_function += cost;
for (int k = 0; k < mu.size(); ++k) {
lagrangian_function += mu[k] * grad_mu[k];
}
model.Minimize(lagrangian_function);
SolveResult result = Solve(model, SolverType::kGscip).value();
CHECK_EQ(result.termination.reason, TerminationReason::kOptimal)
<< result.termination;
const VariableMap<double>& vars_val = result.variable_values();
ASSIGN_OR_RETURN(math_opt::SolveResult result,
Solve(model, math_opt::SolverType::kGscip));
switch (result.termination.reason) {
case math_opt::TerminationReason::kOptimal:
case math_opt::TerminationReason::kFeasible:
break;
default:
return util::InternalErrorBuilder()
<< "failed to minimize lagrangian function: "
<< result.termination;
}
const math_opt::VariableMap<double>& vars_val = result.variable_values();
bool feasible = true;
// Iterate update. Takes a step in the direction of the gradient (since the
@@ -420,35 +447,46 @@ void SolveLagrangianRelaxation(const Graph graph, const double max_resource_1,
upper_bound)
<< std::endl;
std::cout << "========================================" << std::endl;
return absl::OkStatus();
}
void RelaxModel(FlowModel& flow_model) {
for (Variable& var : flow_model.flow_vars) {
for (math_opt::Variable& var : flow_model.flow_vars) {
flow_model.model->set_continuous(var);
flow_model.model->set_lower_bound(var, 0.0);
flow_model.model->set_upper_bound(var, 1.0);
}
}
void SolveFullModel(const Graph& graph, double max_resource_1,
double max_resource_2) {
FlowModel flow_model = SolveMip(graph, max_resource_1, max_resource_2);
absl::Status SolveFullModel(const Graph& graph, double max_resource_1,
double max_resource_2) {
ASSIGN_OR_RETURN(FlowModel flow_model,
SolveMip(graph, max_resource_1, max_resource_2));
RelaxModel(flow_model);
SolveLinearRelaxation(flow_model, graph, max_resource_1, max_resource_2);
return SolveLinearRelaxation(flow_model, graph, max_resource_1,
max_resource_2);
}
} // namespace
int main(int argc, char** argv) {
InitGoogle(argv[0], &argc, &argv, true);
absl::Status Main() {
// Problem data
const Graph graph = CreateSampleNetwork();
const double max_resource_1 = 10;
const double max_resource_2 = 4;
SolveFullModel(graph, max_resource_1, max_resource_2);
RETURN_IF_ERROR(SolveFullModel(graph, max_resource_1, max_resource_2))
<< "full solve failed";
RETURN_IF_ERROR(
SolveLagrangianRelaxation(graph, max_resource_1, max_resource_2))
<< "lagrangian solve failed";
return absl::OkStatus();
}
} // namespace
SolveLagrangianRelaxation(graph, max_resource_1, max_resource_2);
int main(int argc, char** argv) {
InitGoogle(argv[0], &argc, &argv, true);
const absl::Status status = Main();
if (!status.ok()) {
LOG(QFATAL) << status;
}
return 0;
}

50
ortools/math_opt/samples/linear_programming.cc
View File

@@ -24,22 +24,18 @@
#include "absl/time/time.h"
#include "ortools/base/init_google.h"
#include "ortools/base/logging.h"
#include "ortools/base/status_builder.h"
#include "ortools/base/status_macros.h"
#include "ortools/math_opt/cpp/math_opt.h"
namespace {
using ::operations_research::math_opt::LinearConstraint;
using ::operations_research::math_opt::LinearExpression;
using ::operations_research::math_opt::Model;
using ::operations_research::math_opt::SolveResult;
using ::operations_research::math_opt::SolverType;
using ::operations_research::math_opt::Sum;
using ::operations_research::math_opt::TerminationReason;
using ::operations_research::math_opt::Variable;
namespace math_opt = ::operations_research::math_opt;
constexpr double kInf = std::numeric_limits<double>::infinity();
// Model and solve the problem:
// max 10 * x0 + 6 * x1 + 4 * x2
// max 10 * x0 + 6 * x1 + 4 * x2
// s.t. 10 * x0 + 4 * x1 + 5 * x2 <= 600
// 2 * x0 + 2 * x1 + 6 * x2 <= 300
// x0 + x1 + x2 <= 100
@@ -47,17 +43,17 @@ constexpr double kInf = std::numeric_limits<double>::infinity();
// x1 in [0, infinity)
// x2 in [0, infinity)
//
void SolveSimpleLp() {
Model model("Linear programming example");
absl::Status Main() {
math_opt::Model model("Linear programming example");
// Variables
std::vector<Variable> x;
std::vector<math_opt::Variable> x;
for (int j = 0; j < 3; j++) {
x.push_back(model.AddContinuousVariable(0.0, kInf, absl::StrCat("x", j)));
}
// Constraints
std::vector<LinearConstraint> constraints;
std::vector<math_opt::LinearConstraint> constraints;
constraints.push_back(
model.AddLinearConstraint(10 * x[0] + 4 * x[1] + 5 * x[2] <= 600, "c1"));
constraints.push_back(
@@ -68,15 +64,12 @@ void SolveSimpleLp() {
// Objective
model.Maximize(10 * x[0] + 6 * x[1] + 4 * x[2]);
std::cout << "Num variables: " << model.num_variables() << std::endl;
std::cout << "Num constraints: " << model.num_linear_constraints()
<< std::endl;
const SolveResult result = Solve(model, SolverType::kGlop).value();
// Check that the problem has an optimal solution.
QCHECK_EQ(result.termination.reason, TerminationReason::kOptimal)
<< "Failed to find an optimal solution: " << result.termination;
ASSIGN_OR_RETURN(const math_opt::SolveResult result,
Solve(model, math_opt::SolverType::kGlop));
if (result.termination.reason != math_opt::TerminationReason::kOptimal) {
return util::InternalErrorBuilder()
<< "model failed to solve to optimality" << result.termination;
}
std::cout << "Problem solved in " << result.solve_time() << std::endl;
std::cout << "Objective value: " << result.objective_value() << std::endl;
@@ -84,19 +77,16 @@ void SolveSimpleLp() {
std::cout << "Variable values: ["
<< absl::StrJoin(result.variable_values().Values(x), ", ") << "]"
<< std::endl;
std::cout << "Constraint duals: ["
<< absl::StrJoin(result.dual_values().Values(constraints), ", ")
<< "]" << std::endl;
std::cout << "Reduced costs: ["
<< absl::StrJoin(result.reduced_costs().Values(x), ", ") << "]"
<< std::endl;
// TODO(user): add basis statuses when they are included in SolveResult
return absl::OkStatus();
}
} // namespace
int main(int argc, char** argv) {
InitGoogle(argv[0], &argc, &argv, true);
SolveSimpleLp();
const absl::Status status = Main();
if (!status.ok()) {
LOG(QFATAL) << status;
}
return 0;
}

49
ortools/math_opt/solvers/BUILD.bazel
View File

@@ -30,6 +30,7 @@ cc_library(
"//ortools/math_opt:result_cc_proto",
"//ortools/math_opt:solution_cc_proto",
"//ortools/math_opt:sparse_containers_cc_proto",
"//ortools/math_opt/core:inverted_bounds",
"//ortools/math_opt/core:math_opt_proto_utils",
"//ortools/math_opt/core:solve_interrupter",
"//ortools/math_opt/core:solver_interface",
@@ -108,6 +109,7 @@ cc_library(
"//ortools/math_opt:result_cc_proto",
"//ortools/math_opt:solution_cc_proto",
"//ortools/math_opt:sparse_containers_cc_proto",
"//ortools/math_opt/core:inverted_bounds",
"//ortools/math_opt/core:math_opt_proto_utils",
"//ortools/math_opt/core:solver_interface",
"//ortools/math_opt/core:sparse_vector_view",
@@ -153,6 +155,7 @@ cc_library(
"//ortools/math_opt:result_cc_proto",
"//ortools/math_opt:solution_cc_proto",
"//ortools/math_opt:sparse_containers_cc_proto",
"//ortools/math_opt/core:inverted_bounds",
"//ortools/math_opt/core:math_opt_proto_utils",
"//ortools/math_opt/core:solver_interface",
"//ortools/math_opt/core:sparse_vector_view",
@@ -191,6 +194,7 @@ cc_library(
"//ortools/math_opt:result_cc_proto",
"//ortools/math_opt:solution_cc_proto",
"//ortools/math_opt:sparse_containers_cc_proto",
"//ortools/math_opt/core:inverted_bounds",
"//ortools/math_opt/core:math_opt_proto_utils",
"//ortools/math_opt/core:solve_interrupter",
"//ortools/math_opt/core:solver_interface",
@@ -257,6 +261,51 @@ cc_library(
],
)
cc_library(
name = "glpk_solver",
srcs = [
"glpk_solver.cc",
"glpk_solver.h",
],
deps = [
":message_callback_data",
"//ortools/base",
"//ortools/base:cleanup",
"//ortools/base:protoutil",
"//ortools/base:status_macros",
"//ortools/glpk:glpk_env_deleter",
"//ortools/glpk:glpk_formatters",
"//ortools/math_opt:callback_cc_proto",
"//ortools/math_opt:model_cc_proto",
"//ortools/math_opt:model_parameters_cc_proto",
"//ortools/math_opt:model_update_cc_proto",
"//ortools/math_opt:parameters_cc_proto",
"//ortools/math_opt:result_cc_proto",
"//ortools/math_opt:solution_cc_proto",
"//ortools/math_opt:sparse_containers_cc_proto",
"//ortools/math_opt/core:inverted_bounds",
"//ortools/math_opt/core:math_opt_proto_utils",
"//ortools/math_opt/core:solve_interrupter",
"//ortools/math_opt/core:solver_interface",
"//ortools/math_opt/core:sparse_submatrix",
"//ortools/math_opt/core:sparse_vector_view",
"//ortools/math_opt/solvers/glpk:rays",
"//ortools/math_opt/validators:callback_validator",
"//ortools/port:proto_utils",
"@com_google_absl//absl/base:core_headers",
"@com_google_absl//absl/cleanup",
"@com_google_absl//absl/container:flat_hash_map",
"@com_google_absl//absl/memory",
"@com_google_absl//absl/status",
"@com_google_absl//absl/status:statusor",
"@com_google_absl//absl/strings",
"@com_google_absl//absl/synchronization",
"@com_google_absl//absl/time",
"@glpk",
],
alwayslink = 1,
)
proto_library(
name = "gurobi_proto",
srcs = ["gurobi.proto"],

48
ortools/math_opt/solvers/cp_sat_solver.cc
View File

@@ -100,6 +100,9 @@ std::vector<std::string> SetSolveParameters(
request.set_solver_time_limit_seconds(absl::ToDoubleSeconds(
util_time::DecodeGoogleApiProto(parameters.time_limit()).value()));
}
if (parameters.has_node_limit()) {
warnings.push_back("The node_limit parameter is not supported for CP-SAT.");
}
// Build CP SAT parameters by first initializing them from the common
// parameters, and then using the values in `solver_specific_parameters` to
@@ -120,12 +123,12 @@ std::vector<std::string> SetSolveParameters(
if (parameters.has_threads()) {
sat_parameters.set_num_search_workers(parameters.threads());
}
if (parameters.has_relative_gap_limit()) {
sat_parameters.set_relative_gap_limit(parameters.relative_gap_limit());
if (parameters.has_relative_gap_tolerance()) {
sat_parameters.set_relative_gap_limit(parameters.relative_gap_tolerance());
}
if (parameters.has_absolute_gap_limit()) {
sat_parameters.set_absolute_gap_limit(parameters.absolute_gap_limit());
if (parameters.has_absolute_gap_tolerance()) {
sat_parameters.set_absolute_gap_limit(parameters.absolute_gap_tolerance());
}
// cutoff_limit is handled outside this function as it modifies the model.
if (parameters.has_best_bound_limit()) {
@@ -315,6 +318,7 @@ GetTerminationAndStats(const bool is_interrupted, const bool maximize,
}
return std::make_pair(std::move(solve_stats), std::move(termination));
}
} // namespace
absl::StatusOr<std::unique_ptr<SolverInterface>> CpSatSolver::New(
@@ -323,14 +327,18 @@ absl::StatusOr<std::unique_ptr<SolverInterface>> CpSatSolver::New(
MathOptModelToMPModelProto(model));
std::vector variable_ids(model.variables().ids().begin(),
model.variables().ids().end());
std::vector linear_constraint_ids(model.linear_constraints().ids().begin(),
model.linear_constraints().ids().end());
// TODO(b/204083726): Remove this check if QP support is added
if (!model.objective().quadratic_coefficients().row_ids().empty()) {
return absl::InvalidArgumentError(
"MathOpt does not currently support CP-SAT models with quadratic "
"objectives");
}
return absl::WrapUnique(
new CpSatSolver(std::move(cp_sat_model), std::move(variable_ids)));
return absl::WrapUnique(new CpSatSolver(
std::move(cp_sat_model),
/*variable_ids=*/std::move(variable_ids),
/*linear_constraint_ids=*/std::move(linear_constraint_ids)));
}
absl::StatusOr<SolveResultProto> CpSatSolver::Solve(
@@ -435,6 +443,10 @@ absl::StatusOr<SolveResultProto> CpSatSolver::Solve(
// supported by CP-SAT and we have validated they are empty.
};
}
// CP-SAT returns "infeasible" for inverted bounds.
RETURN_IF_ERROR(ListInvertedBounds().ToStatus());
ASSIGN_OR_RETURN(const MPSolutionResponse response,
SatSolveProto(std::move(req), &interrupt_solve,
logging_callback, solution_callback));
@@ -473,9 +485,11 @@ absl::Status CpSatSolver::Update(const ModelUpdateProto& model_update) {
}
CpSatSolver::CpSatSolver(MPModelProto cp_sat_model,
std::vector<int64_t> variable_ids)
std::vector<int64_t> variable_ids,
std::vector<int64_t> linear_constraint_ids)
: cp_sat_model_(std::move(cp_sat_model)),
variable_ids_(std::move(variable_ids)) {}
variable_ids_(std::move(variable_ids)),
linear_constraint_ids_(std::move(linear_constraint_ids)) {}
SparseDoubleVectorProto CpSatSolver::ExtractSolution(
const absl::Span<const double> cp_sat_variable_values,
@@ -498,6 +512,24 @@ SparseDoubleVectorProto CpSatSolver::ExtractSolution(
return result;
}
InvertedBounds CpSatSolver::ListInvertedBounds() const {
InvertedBounds inverted_bounds;
for (int v = 0; v < cp_sat_model_.variable_size(); ++v) {
const MPVariableProto& var = cp_sat_model_.variable(v);
if (var.lower_bound() > var.upper_bound()) {
inverted_bounds.variables.push_back(variable_ids_[v]);
}
}
for (int c = 0; c < cp_sat_model_.constraint_size(); ++c) {
const MPConstraintProto& cstr = cp_sat_model_.constraint(c);
if (cstr.lower_bound() > cstr.upper_bound()) {
inverted_bounds.linear_constraints.push_back(linear_constraint_ids_[c]);
}
}
return inverted_bounds;
}
MATH_OPT_REGISTER_SOLVER(SOLVER_TYPE_CP_SAT, CpSatSolver::New);
} // namespace math_opt

12
ortools/math_opt/solvers/cp_sat_solver.h
View File

@@ -25,6 +25,7 @@
#include "absl/types/span.h"
#include "ortools/linear_solver/linear_solver.pb.h"
#include "ortools/math_opt/callback.pb.h"
#include "ortools/math_opt/core/inverted_bounds.h"
#include "ortools/math_opt/core/solve_interrupter.h"
#include "ortools/math_opt/core/solver_interface.h"
#include "ortools/math_opt/model.pb.h"
@@ -52,18 +53,27 @@ class CpSatSolver : public SolverInterface {
bool CanUpdate(const ModelUpdateProto& model_update) override;
private:
CpSatSolver(MPModelProto cp_sat_model, std::vector<int64_t> variable_ids);
CpSatSolver(MPModelProto cp_sat_model, std::vector<int64_t> variable_ids,
std::vector<int64_t> linear_constraint_ids);
// Extract the solution from CP-SAT's response.
SparseDoubleVectorProto ExtractSolution(
absl::Span<const double> cp_sat_variable_values,
const ModelSolveParametersProto& model_parameters) const;
// Returns the ids of variables and linear constraints with inverted bounds.
InvertedBounds ListInvertedBounds() const;
const MPModelProto cp_sat_model_;
// For the i-th variable in `cp_sat_model_`, `variable_ids_[i]` contains the
// corresponding id in the input `Model`.
const std::vector<int64_t> variable_ids_;
// For the i-th linear constraint in `cp_sat_model_`,
// `linear_constraint_ids_[i]` contains the corresponding id in the input
// `Model`.
const std::vector<int64_t> linear_constraint_ids_;
};
} // namespace math_opt

77
ortools/math_opt/solvers/glop_solver.cc
View File

@@ -67,6 +67,8 @@ namespace math_opt {
namespace {
constexpr double kInf = std::numeric_limits<double>::infinity();
absl::string_view SafeName(const VariablesProto& variables, int index) {
if (variables.names().empty()) {
return {};
@@ -299,6 +301,9 @@ absl::StatusOr<glop::GlopParameters> GlopSolver::MergeSolveParameters(
solver_parameters.iteration_limit()) {
result.set_max_number_of_iterations(solver_parameters.iteration_limit());
}
if (solver_parameters.has_node_limit()) {
warnings.emplace_back("GLOP does snot support 'node_limit' parameter");
}
if (!result.has_use_dual_simplex() &&
solver_parameters.lp_algorithm() != LP_ALGORITHM_UNSPECIFIED) {
switch (solver_parameters.lp_algorithm()) {
@@ -479,6 +484,72 @@ std::vector<int64_t> GetSortedIs(
return sorted;
}
// Returns a vector of containing the MathOpt id of each row or column. Here T
// is either (Col|Row)Index and id_map is expected to be
// GlopSolver::(linear_constraints_|variables_).
template <typename T>
glop::StrictITIVector<T, int64_t> IndexToId(
const absl::flat_hash_map<int64_t, T>& id_map) {
// Guard value used to identify not-yet-set elements of index_to_id.
constexpr int64_t kEmptyId = -1;
glop::StrictITIVector<T, int64_t> index_to_id(T(id_map.size()), kEmptyId);
for (const auto& [id, index] : id_map) {
CHECK(index >= 0 && index < index_to_id.size()) << index;
CHECK_EQ(index_to_id[index], kEmptyId);
index_to_id[index] = id;
}
// At this point, index_to_id can't contain any kEmptyId values since
// index_to_id.size() == id_map.size() and we modified id_map.size() elements
// in the loop, after checking that the modified element was changed by a
// previous iteration.
return index_to_id;
}
InvertedBounds GlopSolver::ListInvertedBounds() const {
// Identify rows and columns by index first.
std::vector<glop::ColIndex> inverted_columns;
const glop::ColIndex num_cols = linear_program_.num_variables();
for (glop::ColIndex col(0); col < num_cols; ++col) {
if (linear_program_.variable_lower_bounds()[col] >
linear_program_.variable_upper_bounds()[col]) {
inverted_columns.push_back(col);
}
}
std::vector<glop::RowIndex> inverted_rows;
const glop::RowIndex num_rows = linear_program_.num_constraints();
for (glop::RowIndex row(0); row < num_rows; ++row) {
if (linear_program_.constraint_lower_bounds()[row] >
linear_program_.constraint_upper_bounds()[row]) {
inverted_rows.push_back(row);
}
}
// Convert column/row indices into MathOpt ids. We avoid calling the expensive
// IndexToId() when not necessary.
InvertedBounds inverted_bounds;
if (!inverted_columns.empty()) {
const glop::StrictITIVector<glop::ColIndex, int64_t> ids =
IndexToId(variables_);
CHECK_EQ(ids.size(), num_cols);
inverted_bounds.variables.reserve(inverted_columns.size());
for (const glop::ColIndex col : inverted_columns) {
inverted_bounds.variables.push_back(ids[col]);
}
}
if (!inverted_rows.empty()) {
const glop::StrictITIVector<glop::RowIndex, int64_t> ids =
IndexToId(linear_constraints_);
CHECK_EQ(ids.size(), num_rows);
inverted_bounds.linear_constraints.reserve(inverted_rows.size());
for (const glop::RowIndex row : inverted_rows) {
inverted_bounds.linear_constraints.push_back(ids[row]);
}
}
return inverted_bounds;
}
void GlopSolver::FillSolution(const glop::ProblemStatus status,
const ModelSolveParametersProto& model_parameters,
SolveResultProto& solve_result) {
@@ -600,7 +671,6 @@ absl::Status GlopSolver::FillSolveStats(const glop::ProblemStatus status,
const absl::Duration solve_time,
SolveStatsProto& solve_stats) {
const bool is_maximize = linear_program_.IsMaximizationProblem();
constexpr double kInf = std::numeric_limits<double>::infinity();
// Set default status and bounds.
solve_stats.mutable_problem_status()->set_primal_status(
@@ -760,6 +830,11 @@ absl::StatusOr<SolveResultProto> GlopSolver::Solve(
}
});
// Glop returns an error when bounds are inverted and does not list the
// offending variables/constraints. Here we want to return a more detailed
// status.
RETURN_IF_ERROR(ListInvertedBounds().ToStatus());
const glop::ProblemStatus status =
lp_solver_.SolveWithTimeLimit(linear_program_, time_limit.get());
const absl::Duration solve_time = absl::Now() - start;

4
ortools/math_opt/solvers/glop_solver.h
View File

@@ -31,6 +31,7 @@
#include "ortools/lp_data/lp_data.h"
#include "ortools/lp_data/lp_types.h"
#include "ortools/math_opt/callback.pb.h"
#include "ortools/math_opt/core/inverted_bounds.h"
#include "ortools/math_opt/core/solve_interrupter.h"
#include "ortools/math_opt/core/solver_interface.h"
#include "ortools/math_opt/model.pb.h"
@@ -82,6 +83,9 @@ class GlopSolver : public SolverInterface {
void UpdateLinearConstraintBounds(
const LinearConstraintUpdatesProto& linear_constraint_updates);
// Returns the ids of variables and linear constraints with inverted bounds.
InvertedBounds ListInvertedBounds() const;
void FillSolution(glop::ProblemStatus status,
const ModelSolveParametersProto& model_parameters,
SolveResultProto& solve_result);

18
ortools/math_opt/solvers/glpk/BUILD.bazel Normal file
View File

@@ -0,0 +1,18 @@
# Code specific to GLPK used in glpk_solver.cc.
package(default_visibility = ["//ortools/math_opt/solvers:__subpackages__"])
cc_library(
name = "rays",
srcs = ["rays.cc"],
hdrs = ["rays.h"],
deps = [
"//ortools/base",
"//ortools/base:status_macros",
"//ortools/glpk:glpk_computational_form",
"//ortools/glpk:glpk_formatters",
"@com_google_absl//absl/status",
"@com_google_absl//absl/status:statusor",
"@com_google_absl//absl/strings",
"@glpk",
],
)

384
ortools/math_opt/solvers/glpk/rays.cc Normal file
View File

@@ -0,0 +1,384 @@
// Copyright 2010-2021 Google LLC
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#include "ortools/math_opt/solvers/glpk/rays.h"
#include <optional>
#include <string>
#include <utility>
#include <vector>
#include "absl/status/status.h"
#include "absl/status/statusor.h"
#include "absl/strings/str_cat.h"
#include "ortools/base/logging.h"
#include "ortools/base/status_macros.h"
#include "ortools/glpk/glpk_computational_form.h"
#include "ortools/glpk/glpk_formatters.h"
namespace operations_research::math_opt {
namespace {
absl::StatusOr<GlpkRay> ComputePrimalRay(glp_prob* const problem,
const int non_basic_variable) {
const int num_cstrs = glp_get_num_rows(problem);
// Check that the non_basic_variable is indeed non basic.
const int non_basic_variable_status =
ComputeFormVarStatus(problem,
/*num_cstrs=*/num_cstrs,
/*k=*/non_basic_variable);
CHECK_NE(non_basic_variable_status, GLP_BS);
// When we perform the (primal) simplex algorithm, we detect the primal
// unboundness when we have a non-basic variable (here variable can be a
// structural or an auxiliary variable) which contributes to increase (for
// maximization, decrease for minimization) the objective but none of the
// basic variables bounds are limiting its growth. GLPK returns the index of
// this non-basic tableau variable.
//
// To be more precise, here we will use the conventions used in
// glpk-5.0/doc/glpk.pdf available from glpk-5.0.tar.gz.
//
// From (glpk eq. 3.13) we know that the values of the basic variables are
// dependent on the values of the non-basic ones:
//
// x_B = 𝚵 x_N
//
// where 𝚵 is the tableau defined by (glpk eq. 3.12):
//
// 𝚵 = -B^-1 N
//
// Thus if the c-th non basic variable is changed:
//
// x'_N = x_N + t e_c , e_c ∈ R^n is the c-th standard unit vector
// t ∈ R is the change
//
// Then to keep the primal feasible we must have:
//
// x'_B = 𝚵 x'_N
// = 𝚵 x_N + t 𝚵 e_c
// = 𝚵 x_N + t 𝚵 e_c
// = x_B + t 𝚵 e_c
//
// We thus have the primal ray:
//
// x'_N - x_N = t e_c
// x'_B - x_B = t 𝚵 e_c
//
// From (glpk eq. 3.34) we know that the primal objective is:
//
// z = d^T x_N + c_0
//
// I.e. reduced cost d_j shows how the non-basic variable x_j influences the
// objective.
//
// Thus if the problem is a minimization we know that:
//
// t > 0 , if d_c < 0
// t < 0 , if d_c > 0
//
// Since if it was not the case, the primal simplex algorithm would not have
// picked this variable.
//
// The signs for a maximization are reversed:
//
// t < 0 , if d_c < 0
// t > 0 , if d_c > 0
const double reduced_cost = ComputeFormVarReducedCost(
problem, /*num_cstrs=*/num_cstrs, /*k=*/non_basic_variable);
const double t = (glp_get_obj_dir(problem) == GLP_MAX ? 1.0 : -1.0) *
(reduced_cost >= 0 ? 1.0 : -1.0);
// In case of bounded variables, we can check that the result agrees with the
// current active bound. We can't do so for free variables though.
switch (non_basic_variable_status) {
case GLP_NL: // At lower-bound.
if (t == -1.0) {
return absl::InternalError(
"a non-basic variable at its lower-bound is reported as cause of "
"unboundness but the reduced cost's sign indicates that the solver "
"considered making it smaller");
}
break;
case GLP_NU: // At upper-bound.
if (t == 1.0) {
return absl::InternalError(
"a non-basic variable at its upper-bound is reported as cause of "
"unboundness but the reduced cost's sign indicates that the solver "
"considered making it bigger");
}
break;
case GLP_NF: // Free (unbounded).
break;
default: // GLP_BS (basic), GLP_NS (fixed) or invalid value
return absl::InternalError(absl::StrCat(
"unexpected ", BasisStatusString(non_basic_variable_status),
" reported as cause of unboundness"));
}
GlpkRay::SparseVector ray_non_zeros;
// As seen in the maths above:
//
// x'_N - x_N = t e_c
//
ray_non_zeros.emplace_back(non_basic_variable, t);
// As seen in the maths above:
//
// x'_B - x_B = t 𝚵 e_c
//
// Here 𝚵 e_c is the c-th column of the tableau. We thus use the GLPK function
// that returns this column.
std::vector<int> inds(num_cstrs + 1);
std::vector<double> vals(num_cstrs + 1);
const int non_zeros =
glp_eval_tab_col(problem, non_basic_variable, inds.data(), vals.data());
for (int i = 1; i <= non_zeros; ++i) {
ray_non_zeros.emplace_back(inds[i], t * vals[i]);
}
return GlpkRay(GlpkRayType::kPrimal, std::move(ray_non_zeros));
}
absl::StatusOr<GlpkRay> ComputeDualRay(glp_prob* const problem,
const int basic_variable) {
const int num_cstrs = glp_get_num_rows(problem);
// Check that the basic_variable is indeed basic.
{
const int status = ComputeFormVarStatus(problem,
/*num_cstrs=*/num_cstrs,
/*k=*/basic_variable);
CHECK_EQ(status, GLP_BS) << BasisStatusString(status);
}
// The dual simplex proceeds by repeatedly finding basic variables (here
// variable includes structural and auxiliary variables) that are primal
// infeasible and replacing them in the basis with a non-basic variable whose
// growth is limited by their reduced cost.
//
// This algorithm detects dual unboundness when we have a basic variable is
// primal infeasible (out of its bounds) but there are no non-basic variable
// that would limit the growth of its reduced cost, and thus the growth of the
// dual objective.
//
// To be more precise, here we will use the conventions used in
// glpk-5.0/doc/glpk.pdf available from glpk-5.0.tar.gz. The dual simplex
// algorithm is defined by (https://d-nb.info/978580478/34): Koberstein,
// Achim. "The dual simplex method, techniques for a fast and stable
// implementation." Unpublished doctoral thesis, Universität Paderborn,
// Paderborn, Germany (2005).
//
// In the following reasoning, we will considering the dual after the
// permutation of the basis (glpk eq. 3.27):
//
// B^T π + λ_B = c_B
// N^T π + λ_N = c_N
//
// We will now see what happens when we relax a basic variable that would
// leave the base. See (Koberstein §3.1.2) for details.
//
// Let's assume we have (π, λ_B, λ_N) that is a basic dual feasible
// solution. By definition:
//
// λ_B = 0
//
// If we relax the equality constraint of the basic variable r that is primal
// infeasible, that is if we relax λ_B_r and get another solution (π', λ'_B,
// λ'_N). By definition, all other basic variables stays at equality and thus:
//
// λ'_B = t e_r , e_r ∈ R^m is the standard unit vector
// t ∈ R is the relaxation
//
// From (glpk eq. 3.30) we have:
//
// λ'_N = N^T B^-T λ'_B + (c_N - N^T B^-T c_B)
// λ'_N = t (B^-1 N)^T e_r + λ_N
//
// Using the (glpk eq. 3.12) definition of the tableau:
//
// 𝚵 = -B^-1 N
//
// We have:
//
// λ'_N = -t 𝚵^T e_r + λ_N
//
// That is that the change of the reduced cost of the basic variable r has to
// be compensated by the change of the reduced costs the non-basic variables.
//
// We can write the new dual objective:
//
// Z' = l^T λ'_l + u^T λ'_u
//
// If the problem is a minimization we have:
//
// Z' = sum_{j:λ'_N_j >= 0} l_N_j λ'_N_j +
// sum_{j:λ'_N_j <= 0} u_N_j λ'_N_j +
// {l_B_r, if t >= 0, u_B_r, else} t
//
// Here we assume the signs of λ'_N are identical to the ones of λ_N (this is
// not an issue with dual simplex since we want to make one non-basic tight to
// use it in the basis) we can replace λ'_N with the value computed above and
// considering the initial solution was basic which implied that non-basic
// where at their bound we can rewrite the objective as:
//
// Z' = Z - t e_r^T 𝚵 x_N + {l_B_r, if t >= 0, u_B_r, else} t
//
// We have, using (glpk eq. 3.13):
//
// e_r^T 𝚵 x_N = e_r^T x_B = x_B_r
//
// And thus, for a minimization we have:
//
// Z' - Z = t * {l_B_r - x_B_r, if t >= 0,
// u_B_r - x_B_r, if t <= 0}
//
// Depending on the type of constraint, i.e. depending on whether l_B_r and/or
// u_B_r are finite), we have constraints on the sign of `t`. But we can see
// that since we pick the basic variable r because it was primal infeasible,
// then it should break one of its finite bounds.
//
// either x_B_r < l_B_r
// or u_B_r < x_B_r
//
// If l_B_r is finite and x_B_r < l_B_r, then choosing:
//
// t >= 0
//
// leads to:
//
// Z' - Z >= 0
//
// and we see from (glpk eq. 3.17) and the "rule of signs" table (glpk page
// 101) that we keep the solution dual feasible by doing so.
//
// The same logic applies if x_B_r > u_B_r:
//
// t <= 0
//
// leads to:
//
// Z' - Z >= 0
//
// The dual objective increase in both cases; which is what we want for a
// minimization problem since the dual is a maximization.
//
// For a maximization problem the results are similar but the sign of t
// changes (which is expected since the dual is a minimization):
//
// Z' - Z = t * {l_B_r - x_B_r, if t <= 0,
// u_B_r - x_B_r, if t >= 0}
//
// If a problem is dual unbounded, this means that it is possible to grow t
// without limit. I.e. is possible to choose any value for t without making
// any λ'_N change sign.
//
// We can then express the changes of λ' from t:
//
// λ'_B = t e_r
// λ'_N = -t 𝚵^T e_r + λ_N
//
// Since λ_B = 0, we can rewrite those as:
//
// λ'_B - λ_B = t e_r
// λ'_N - λ_N = -t 𝚵^T e_r
//
// That is the dual ray.
const double primal_value = ComputeFormVarPrimalValue(
problem, /*num_cstrs=*/num_cstrs, /*k=*/basic_variable);
const double upper_bound = ComputeFormVarUpperBound(
problem, /*num_cstrs=*/num_cstrs, /*k=*/basic_variable);
const double lower_bound = ComputeFormVarLowerBound(
problem, /*num_cstrs=*/num_cstrs, /*k=*/basic_variable);
if (!(primal_value > upper_bound || primal_value < lower_bound)) {
return absl::InternalError(
"dual ray computation failed: GLPK identified a basic variable as the "
"source of unboundness but its primal value is within its bounds");
}
// As we have seen in the maths above, depending on which primal bound is
// violated and the optimization direction, we choose the sign of t.
//
// Here the problem is unbounded so we can pick any value for t we want.
const double t = (glp_get_obj_dir(problem) == GLP_MAX ? 1.0 : -1.0) *
(primal_value > upper_bound ? 1.0 : -1.0);
GlpkRay::SparseVector ray_non_zeros;
// As seen in the math above:
//
// λ'_B - λ_B = t e_r
//
ray_non_zeros.emplace_back(basic_variable, t);
// As we have seen above, to keep the dual feasible, we must update the
// reduced costs of the non-basic variables by the formula:
//
// λ'_N - λ_N = -t 𝚵^T e_r
//
// Here 𝚵^T e_r is the r-th row of the tableau. We thus use the GLPK function
// that returns this row.
const int num_structural_vars = glp_get_num_cols(problem);
std::vector<int> inds(num_structural_vars + 1);
std::vector<double> vals(num_structural_vars + 1);
const int non_zeros =
glp_eval_tab_row(problem, basic_variable, inds.data(), vals.data());
for (int i = 1; i <= non_zeros; ++i) {
ray_non_zeros.emplace_back(inds[i], -t * vals[i]);
}
return GlpkRay(GlpkRayType::kDual, std::move(ray_non_zeros));
}
} // namespace
GlpkRay::GlpkRay(const GlpkRayType type, SparseVector non_zero_components)
: type(type), non_zero_components(std::move(non_zero_components)) {}
absl::StatusOr<std::optional<GlpkRay>> GlpkComputeUnboundRay(
glp_prob* const problem) {
const int unbound_ray = glp_get_unbnd_ray(problem);
if (unbound_ray <= 0) {
// No ray, do nothing.
DCHECK_EQ(unbound_ray, 0);
return std::nullopt;
}
// The factorization may not exists when GLPK's trivial_lp() is used to solve
// a trivial LP. Here we force the computation of the factorization if
// necessary.
if (!glp_bf_exists(problem)) {
const int factorization_rc = glp_factorize(problem);
if (factorization_rc != 0) {
return util::InternalErrorBuilder() << "glp_factorize() failed: "
<< ReturnCodeString(factorization_rc);
}
}
// The function glp_get_unbnd_ray() returns either:
// - a non-basic tableau variable if we have primal unboundness.
// - a basic tableau variable if we have dual unboundness.
const bool is_dual_ray =
ComputeFormVarStatus(problem,
/*num_cstrs=*/glp_get_num_rows(problem),
/*k=*/unbound_ray) == GLP_BS;
ASSIGN_OR_RETURN(const GlpkRay ray,
(is_dual_ray ? ComputeDualRay(problem, unbound_ray)
: ComputePrimalRay(problem, unbound_ray)));
return ray;
}
} // namespace operations_research::math_opt

85
ortools/math_opt/solvers/glpk/rays.h Normal file
View File

@@ -0,0 +1,85 @@
// Copyright 2010-2021 Google LLC
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
// This header defines primal/dual unboundness ray computation functions for
// GLPK solver. They use the index space of the computation form of the model as
// defined in operations_research/glpk/glpk_computational_form.h.
#ifndef OR_TOOLS_MATH_OPT_SOLVERS_GLPK_RAYS_H_
#define OR_TOOLS_MATH_OPT_SOLVERS_GLPK_RAYS_H_
#include <optional>
#include <utility>
#include <vector>
#include "absl/status/statusor.h"
extern "C" {
#include <glpk.h>
}
namespace operations_research::math_opt {
// The type of the GlpkRay.
enum GlpkRayType {
// A primal ray.
//
// If x (vector of variables) is a primal feasible solution to a primal
// unbounded problem and r is the ray, then x' = x + t r is also a primal
// feasible solution for all t >= 0.
kPrimal,
// A dual ray.
//
// If λ (vector of reduced costs) is a dual feasible solution to a dual
// unbounded problem and r is the ray, then λ' = λ + t r is also a dual
// feasible solution for all t >= 0.
kDual,
};
// A primal or dual unbound ray for the model in computational form.
//
// See the top comment of operations_research/glpk/glpk_computational_form.h to
// understand that the computational form is. This structure uses the word
// "variable" to mean a variable in the joint set of structural and auxiliary
// variables.
struct GlpkRay {
using SparseVector = std::vector<std::pair<int, double>>;
GlpkRay(GlpkRayType type, SparseVector non_zero_components);
// The type of ray, primal or dual.
GlpkRayType type;
// The non zero components of the vector, in no particular order.
//
// The first member of the pair is the index of the variable (or of its
// corresponding reduced cost) and the second is the component's value.
//
// A given index can only appear once.
//
// The indices in GLPK are one-based. Here the indices are defined by:
// - if 1 <= k <= m: k is the index of the k-th auxiliary variable
// (a.k.a. row, a.k.a. constraint)
// - if m + 1 <= k <= m + n: k is the index of the (k-m)-th structural
// variable (a.k.a. column)
// Note that the value k = 0 is not used.
SparseVector non_zero_components;
};
// Returns the primal or dual ray if one is identified by
// glp_get_unbnd_ray(). Returns an error status if an internal error occurs.
absl::StatusOr<std::optional<GlpkRay>> GlpkComputeUnboundRay(glp_prob* problem);
} // namespace operations_research::math_opt
#endif // OR_TOOLS_MATH_OPT_SOLVERS_GLPK_RAYS_H_

1748
ortools/math_opt/solvers/glpk_solver.cc Normal file
View File

File diff suppressed because it is too large Load Diff

221
ortools/math_opt/solvers/glpk_solver.h Normal file
View File

@@ -0,0 +1,221 @@
// Copyright 2010-2021 Google LLC
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
#ifndef OR_TOOLS_MATH_OPT_SOLVERS_GLPK_SOLVER_H_
#define OR_TOOLS_MATH_OPT_SOLVERS_GLPK_SOLVER_H_
#include <cstdint>
#include <memory>
#include <optional>
#include <vector>
#include "absl/container/flat_hash_map.h"
#include "absl/status/status.h"
#include "absl/status/statusor.h"
#include "ortools/base/logging.h"
#include "ortools/math_opt/callback.pb.h"
#include "ortools/math_opt/core/solve_interrupter.h"
#include "ortools/math_opt/core/solver_interface.h"
#include "ortools/math_opt/model.pb.h"
#include "ortools/math_opt/model_parameters.pb.h"
#include "ortools/math_opt/model_update.pb.h"
#include "ortools/math_opt/parameters.pb.h"
#include "ortools/math_opt/result.pb.h"
extern "C" {
#include <glpk.h>
}
namespace operations_research {
namespace math_opt {
class GlpkSolver : public SolverInterface {
public:
static absl::StatusOr<std::unique_ptr<SolverInterface>> New(
const ModelProto& model, const InitArgs& init_args);
~GlpkSolver() override;
absl::StatusOr<SolveResultProto> Solve(
const SolveParametersProto& parameters,
const ModelSolveParametersProto& model_parameters,
MessageCallback message_cb,
const CallbackRegistrationProto& callback_registration, Callback cb,
SolveInterrupter* interrupter) override;
absl::Status Update(const ModelUpdateProto& model_update) override;
bool CanUpdate(const ModelUpdateProto& model_update) override;
private:
// The columns of the GPLK problem.
//
// This structures is intentionally similar to LinearConstraints so that some
// template function can accept either of those to share code between rows and
// columns. For this purpose it also defines some aliases to some GLPK
// functions and the IsInteger() (which is trivially implemented for
// LinearConstraints).
struct Variables {
static constexpr auto kSetBounds = glp_set_col_bnds;
static constexpr auto kGetLb = glp_get_col_lb;
static constexpr auto kGetUb = glp_get_col_ub;
static constexpr auto kGetType = glp_get_col_type;
static constexpr auto kDelElts = glp_del_cols;
// Returns true if the given one-based column is an integer variable.
static inline bool IsInteger(glp_prob* const problem, const int j);
// The MathOpt variable id of each column in GLPK. This is zero-based, the
// first column corresponds to the 0 and the ids.size() matches the number
// of columns.
//
// The id_to_index map can be used to get the GLPK column index of a given
// MathOpt variable id but the return value will be one-based (the
// convention used in GLPK). Thus this invariant holds:
//
// for all i in [0, num_cols), id_to_index.at(ids[i]) == i + 1
std::vector<int64_t> ids;
// Map each MathOpt variable id to the column one-based index in GLPK (thus
// values are in [1, num_cols]). See the ids vector for the counter part.
absl::flat_hash_map<int64_t, int> id_to_index;
// The unrounded lower bound value of each column.
//
// We keep this value since GLPK's glp_intopt() expects integer bounds for
// integer variables. We need the unrounded value when the type of a
// variable is changed to continuous though by an update.
std::vector<double> unrounded_lower_bounds;
// The unrounded upper bound value of each column.
//
// See unrounded_lower_bounds documentation for details..
std::vector<double> unrounded_upper_bounds;
};
// The columns of the GPLK problem.
//
// See the comment on Variables for details.
struct LinearConstraints {
static constexpr auto kSetBounds = glp_set_row_bnds;
static constexpr auto kGetLb = glp_get_row_lb;
static constexpr auto kGetUb = glp_get_row_ub;
static constexpr auto kGetType = glp_get_row_type;
static constexpr auto kDelElts = glp_del_rows;
// Returns false. This function mirrors Variables::IsInteger() and enables
// sharing code between variables and constraints.
static bool IsInteger(glp_prob*, int) { return false; }
// The MathOpt linear constraint id of each row in GLPK. This is zero-based,
// the first row corresponds to the 0 and ids.size() matches the number of
// rows.
//
// The id_to_index map can be used to get the GLPK row index of a given
// MathOpt variable id but the return value will be one-based (the
// convention used in GLPK). Thus this invariant holds:
//
// for all i in [0, num_rows), id_to_index.at(ids[i]) == i + 1
std::vector<int64_t> ids;
// Map each MathOpt linear constraint id to the row one-based index in GLPK
// (thus values are in [1, num_rows]). See the ids vector for the counter
// part.
absl::flat_hash_map<int64_t, int> id_to_index;
};
explicit GlpkSolver(const ModelProto& model);
// Appends the variables to GLPK cols.
void AddVariables(const VariablesProto& new_variables);
// Appends the variables to GLPK rows.
void AddLinearConstraints(
const LinearConstraintsProto& new_linear_constraints);
// Updates the objective coefficients with the new values in
// coefficients_proto.
void UpdateObjectiveCoefficients(
const SparseDoubleVectorProto& coefficients_proto);
// Updates the constraints matrix with the new values in matrix_updates.
//
// The first_new_(var|cstr)_id are the smallest ids of the new
// variables/constraints (in MathOpt the same id is never reused thus all
// variables with ids greater or equal to these values are new). A nullopt
// value means that there are not new variables/constraints.
void UpdateLinearConstraintMatrix(
const SparseDoubleMatrixProto& matrix_updates,
std::optional<int64_t> first_new_var_id,
std::optional<int64_t> first_new_cstr_id);
// Adds the primal solution (if it exists) to the result using the provided
// functions to get the status of the solution (GLP_FEAS, ...), its
// objective value and the structural variables values.
//
// Here `col_val` is a functions that takes a column index (i.e. the index of
// a structural variable) and returns its primal value in the solution.
void AddPrimalSolution(int (*get_prim_stat)(glp_prob*),
double (*obj_val)(glp_prob*),
double (*col_val)(glp_prob*, int),
const ModelSolveParametersProto& model_parameters,
SolutionProto& solution_proto);
// Adds the dual solution (if it exists) to the result. This function must
// only be called after having solved an LP, with the provided methods
// depending on the type of LP solved.
//
// Here `col_dual` is a function that takes a column index (i.e. the index of
// a structural variable) and returns its dual value in the solution. The
// `row_dual` does the same for a row index (i.e. the index of an auxiliary
// variable associated to a constraint).
void AddDualSolution(int (*get_dual_stat)(glp_prob*),
double (*obj_val)(glp_prob*),
double (*row_dual)(glp_prob*, int),
double (*col_dual)(glp_prob*, int),
const ModelSolveParametersProto& model_parameters,
SolutionProto& solution_proto);
// Adds a primal or dual ray to the result depending on the value returned by
// glp_get_unbnd_ray().
absl::Status AddPrimalOrDualRay(
const ModelSolveParametersProto& model_parameters,
SolveResultProto& result);
glp_prob* const problem_;
Variables variables_;
LinearConstraints linear_constraints_;
};
////////////////////////////////////////////////////////////////////////////////
// Inline functions implementation.
////////////////////////////////////////////////////////////////////////////////
bool GlpkSolver::Variables::IsInteger(glp_prob* const problem, const int j) {
const int kind = glp_get_col_kind(problem, j);
switch (kind) {
case GLP_IV:
case GLP_BV:
// GLP_BV is returned when the GLPK internal kind is GLP_IV and the
// bounds are [0,1].
return true;
case GLP_CV:
return false;
default:
LOG(FATAL) << "Unexpected column kind: " << kind;
}
}
} // namespace math_opt
} // namespace operations_research
#endif // OR_TOOLS_MATH_OPT_SOLVERS_GLPK_SOLVER_H_

58
ortools/math_opt/solvers/gscip_solver.cc
View File

@@ -122,8 +122,8 @@ absl::flat_hash_map<int64_t, double> SparseDoubleVectorAsMap(
// order, does a linear scan from index scan_start to find the index of the
// first entry with row >= row_id. Returns the size the tuple list if there is
// no such entry.
inline int FindRowStart(const SparseDoubleMatrixProto& matrix, const int row_id,
const int scan_start) {
inline int FindRowStart(const SparseDoubleMatrixProto& matrix,
const int64_t row_id, const int scan_start) {
int result = scan_start;
while (result < matrix.row_ids_size() && matrix.row_ids(result) < row_id) {
++result;
@@ -178,11 +178,11 @@ class LinearConstraintIterator {
result.name = SafeName(*linear_constraints_, current_con_);
result.linear_constraint_id = SafeId(*linear_constraints_, current_con_);
const auto vars_begin = linear_constraint_matrix_->column_ids().begin();
const auto vars_begin = linear_constraint_matrix_->column_ids().data();
result.variable_ids = absl::MakeConstSpan(vars_begin + matrix_start_,
vars_begin + matrix_end_);
const auto coefficients_begins =
linear_constraint_matrix_->coefficients().begin();
linear_constraint_matrix_->coefficients().data();
result.coefficients = absl::MakeConstSpan(
coefficients_begins + matrix_start_, coefficients_begins + matrix_end_);
return result;
@@ -433,14 +433,18 @@ absl::StatusOr<GScipParameters> GScipSolver::MergeParameters(
GScipSetMaxNumThreads(solve_parameters.threads(), &result);
}
if (solve_parameters.has_relative_gap_limit()) {
if (solve_parameters.has_relative_gap_tolerance()) {
(*result.mutable_real_params())["limits/gap"] =
solve_parameters.relative_gap_limit();
solve_parameters.relative_gap_tolerance();
}
if (solve_parameters.has_absolute_gap_limit()) {
if (solve_parameters.has_absolute_gap_tolerance()) {
(*result.mutable_real_params())["limits/absgap"] =
solve_parameters.absolute_gap_limit();
solve_parameters.absolute_gap_tolerance();
}
if (solve_parameters.has_node_limit()) {
(*result.mutable_long_params())["limits/totalnodes"] =
solve_parameters.node_limit();
}
if (solve_parameters.has_objective_limit()) {
@@ -688,13 +692,6 @@ absl::StatusOr<SolveResultProto> GScipSolver::CreateSolveResultProto(
GScipResult gscip_result, const ModelSolveParametersProto& model_parameters,
const std::optional<double> cutoff) {
SolveResultProto solve_result;
ASSIGN_OR_RETURN(
*solve_result.mutable_termination(),
ConvertTerminationReason(
gscip_result.gscip_output.status(),
gscip_result.gscip_output.status_detail(),
/*has_feasible_solution=*/!gscip_result.solutions.empty(),
/*had_cutoff=*/cutoff.has_value()));
const bool is_maximize = gscip_->ObjectiveIsMaximize();
// When an objective limit is set, SCIP returns the solutions worse than the
// limit, we need to filter these out manually.
@@ -740,6 +737,12 @@ absl::StatusOr<SolveResultProto> GScipSolver::CreateSolveResultProto(
model_parameters.variable_values_filter());
}
const bool has_feasible_solution = solve_result.solutions_size() > 0;
ASSIGN_OR_RETURN(
*solve_result.mutable_termination(),
ConvertTerminationReason(gscip_result.gscip_output.status(),
gscip_result.gscip_output.status_detail(),
/*has_feasible_solution=*/has_feasible_solution,
/*had_cutoff=*/cutoff.has_value()));
*solve_result.mutable_solve_stats()->mutable_problem_status() =
GetProblemStatusProto(
gscip_result.gscip_output.status(),
@@ -838,6 +841,9 @@ absl::StatusOr<SolveResultProto> GScipSolver::Solve(
const auto interrupter_cleanup = absl::MakeCleanup(
[&]() { interrupt_event_handler_.interrupter = nullptr; });
// SCIP returns "infeasible" when the model contain invalid bounds.
RETURN_IF_ERROR(ListInvertedBounds().ToStatus());
ASSIGN_OR_RETURN(GScipResult gscip_result,
gscip_->Solve(gscip_parameters,
/*legacy_params=*/"",
@@ -873,6 +879,28 @@ absl::flat_hash_set<SCIP_VAR*> GScipSolver::LookupAllVariables(
return result;
}
// Returns the ids of variables and linear constraints with inverted bounds.
InvertedBounds GScipSolver::ListInvertedBounds() const {
// Get the SCIP variables/constraints with inverted bounds.
InvertedBounds inverted_bounds;
for (const auto& [id, var] : variables_) {
if (gscip_->Lb(var) > gscip_->Ub(var)) {
inverted_bounds.variables.push_back(id);
}
}
for (const auto& [id, cstr] : linear_constraints_) {
if (gscip_->LinearConstraintLb(cstr) > gscip_->LinearConstraintUb(cstr)) {
inverted_bounds.linear_constraints.push_back(id);
}
}
// Above code have inserted ids in non-stable order.
std::sort(inverted_bounds.variables.begin(), inverted_bounds.variables.end());
std::sort(inverted_bounds.linear_constraints.begin(),
inverted_bounds.linear_constraints.end());
return inverted_bounds;
}
bool GScipSolver::CanUpdate(const ModelUpdateProto& model_update) {
return gscip_
->CanSafeBulkDelete(

4
ortools/math_opt/solvers/gscip_solver.h
View File

@@ -27,6 +27,7 @@
#include "ortools/gscip/gscip.pb.h"
#include "ortools/gscip/gscip_event_handler.h"
#include "ortools/math_opt/callback.pb.h"
#include "ortools/math_opt/core/inverted_bounds.h"
#include "ortools/math_opt/core/solve_interrupter.h"
#include "ortools/math_opt/core/solver_interface.h"
#include "ortools/math_opt/model.pb.h"
@@ -128,6 +129,9 @@ class GScipSolver : public SolverInterface {
const ModelSolveParametersProto& model_parameters,
std::optional<double> cutoff);
// Returns the ids of variables and linear constraints with inverted bounds.
InvertedBounds ListInvertedBounds() const;
const std::unique_ptr<GScip> gscip_;
InterruptEventHandler interrupt_event_handler_;
absl::flat_hash_map<int64_t, SCIP_VAR*> variables_;

2
ortools/math_opt/solvers/gurobi.proto
View File

@@ -21,7 +21,7 @@ message GurobiInitializerProto {
message ISVKey {
string name = 1;
string application_name = 2;
int64 expiration = 3;
int32 expiration = 3;
string key = 4;
}

2
ortools/math_opt/solvers/gurobi/g_gurobi.h
View File

@@ -52,7 +52,7 @@ namespace operations_research::math_opt {
struct GurobiIsvKey {
std::string name;
std::string application_name;
int64_t expiration = 0;
int32_t expiration = 0;
std::string key;
};

1
ortools/math_opt/solvers/gurobi_callback.cc
View File

@@ -50,7 +50,6 @@ namespace {
// at, see the table here:
// https://www.gurobi.com/documentation/9.1/refman/cb_codes.html
constexpr int kNumGurobiEvents = 9;
constexpr int kGrbOk = 0;
constexpr double kInf = std::numeric_limits<double>::infinity();
template <int where>

71
ortools/math_opt/solvers/gurobi_solver.cc
View File

@@ -68,8 +68,6 @@ namespace operations_research {
namespace math_opt {
namespace {
constexpr int kGrbOk = 0;
absl::StatusOr<std::unique_ptr<Gurobi>> GurobiFromInitArgs(
const SolverInterface::InitArgs& init_args) {
// We don't test or return an error for incorrect non streamable arguments
@@ -146,6 +144,13 @@ GurobiParametersProto MergeParameters(
parameter->set_value(absl::StrCat(time_limit));
}
if (solve_parameters.has_node_limit()) {
GurobiParametersProto::Parameter* const parameter =
merged_parameters.add_parameters();
parameter->set_name(GRB_DBL_PAR_NODELIMIT);
parameter->set_value(absl::StrCat(solve_parameters.node_limit()));
}
if (solve_parameters.has_threads()) {
const int threads = solve_parameters.threads();
GurobiParametersProto::Parameter* const parameter =
@@ -154,20 +159,22 @@ GurobiParametersProto MergeParameters(
parameter->set_value(absl::StrCat(threads));
}
if (solve_parameters.has_absolute_gap_limit()) {
const double absolute_gap_limit = solve_parameters.absolute_gap_limit();
if (solve_parameters.has_absolute_gap_tolerance()) {
const double absolute_gap_tolerance =
solve_parameters.absolute_gap_tolerance();
GurobiParametersProto::Parameter* const parameter =
merged_parameters.add_parameters();
parameter->set_name(GRB_DBL_PAR_MIPGAPABS);
parameter->set_value(absl::StrCat(absolute_gap_limit));
parameter->set_value(absl::StrCat(absolute_gap_tolerance));
}
if (solve_parameters.has_relative_gap_limit()) {
const double relative_gap_limit = solve_parameters.relative_gap_limit();
if (solve_parameters.has_relative_gap_tolerance()) {
const double relative_gap_tolerance =
solve_parameters.relative_gap_tolerance();
GurobiParametersProto::Parameter* const parameter =
merged_parameters.add_parameters();
parameter->set_name(GRB_DBL_PAR_MIPGAP);
parameter->set_value(absl::StrCat(relative_gap_limit));
parameter->set_value(absl::StrCat(relative_gap_tolerance));
}
if (solve_parameters.has_cutoff_limit()) {
@@ -1373,7 +1380,7 @@ absl::Status GurobiSolver::UpdateDoubleListAttribute(
}
std::vector<int> index;
index.reserve(update.ids_size());
for (const int id : update.ids()) {
for (const int64_t id : update.ids()) {
index.push_back(id_hash_map.at(id));
}
return gurobi_->SetDoubleAttrList(attribute_name, index, update.values());
@@ -1387,7 +1394,7 @@ absl::Status GurobiSolver::UpdateInt32ListAttribute(
}
std::vector<int> index;
index.reserve(update.ids_size());
for (const int id : update.ids()) {
for (const int64_t id : update.ids()) {
index.push_back(id_hash_map.at(id));
}
return gurobi_->SetIntAttrList(attribute_name, index, update.values());
@@ -1599,11 +1606,13 @@ absl::Status GurobiSolver::UpdateLinearConstraints(
}
}
// If the constraint had a slack, and now is marked for deletion, we reset
// the stored slack_index in linear_constraints_map_[id] and save the index
// in the list of variables to be deleted later on.
// the stored slack_index in linear_constraints_map_[id], save the index
// in the list of variables to be deleted later on and remove the constraint
// from slack_map_.
if (delete_slack && update_data.source.slack_index != kUnspecifiedIndex) {
deleted_variables_index.emplace_back(update_data.source.slack_index);
update_data.source.slack_index = kUnspecifiedIndex;
slack_map_.erase(update_data.constraint_id);
}
}
@@ -1717,7 +1726,7 @@ absl::Status GurobiSolver::Update(const ModelUpdateProto& model_update) {
model_update.variable_updates().integers();
std::vector<GurobiVariableIndex> index;
index.reserve(update.ids_size());
for (int id : update.ids()) {
for (const int64_t id : update.ids()) {
index.push_back(variables_map_.at(id));
}
std::vector<char> value;
@@ -1860,6 +1869,34 @@ GurobiSolver::RegisterCallback(const CallbackRegistrationProto& registration,
local_interrupter);
}
absl::StatusOr<InvertedBounds> GurobiSolver::ListInvertedBounds() const {
InvertedBounds inverted_bounds;
{
ASSIGN_OR_RETURN(
const std::vector<double> var_lbs,
gurobi_->GetDoubleAttrArray(GRB_DBL_ATTR_LB, num_gurobi_variables_));
ASSIGN_OR_RETURN(
const std::vector<double> var_ubs,
gurobi_->GetDoubleAttrArray(GRB_DBL_ATTR_UB, num_gurobi_variables_));
for (const auto& [id, index] : variables_map_) {
if (var_lbs[index] > var_ubs[index]) {
inverted_bounds.variables.push_back(id);
}
}
}
for (const auto& [id, cstr_data] : linear_constraints_map_) {
if (cstr_data.lower_bound > cstr_data.upper_bound) {
inverted_bounds.linear_constraints.push_back(id);
}
}
// Above code have inserted ids in non-stable order.
std::sort(inverted_bounds.variables.begin(), inverted_bounds.variables.end());
std::sort(inverted_bounds.linear_constraints.begin(),
inverted_bounds.linear_constraints.end());
return inverted_bounds;
}
absl::StatusOr<SolveResultProto> GurobiSolver::Solve(
const SolveParametersProto& parameters,
const ModelSolveParametersProto& model_parameters,
@@ -1934,6 +1971,14 @@ absl::StatusOr<SolveResultProto> GurobiSolver::Solve(
gurobi_cb_data->local_interrupter);
};
}
// Gurobi returns "infeasible" when bounds are inverted.
{
ASSIGN_OR_RETURN(const InvertedBounds inverted_bounds,
ListInvertedBounds());
RETURN_IF_ERROR(inverted_bounds.ToStatus());
}
RETURN_IF_ERROR(gurobi_->Optimize(grb_cb));
// We flush message callbacks before testing for Gurobi error in case where

4
ortools/math_opt/solvers/gurobi_solver.h
View File

@@ -31,6 +31,7 @@
#include "ortools/base/logging.h"
#include "ortools/gurobi/environment.h"
#include "ortools/math_opt/callback.pb.h"
#include "ortools/math_opt/core/inverted_bounds.h"
#include "ortools/math_opt/core/solve_interrupter.h"
#include "ortools/math_opt/core/solver_interface.h"
#include "ortools/math_opt/model.pb.h"
@@ -227,6 +228,9 @@ class GurobiSolver : public SolverInterface {
const MessageCallback message_cb, absl::Time start,
SolveInterrupter* interrupter);
// Returns the ids of variables and linear constraints with inverted bounds.
absl::StatusOr<InvertedBounds> ListInvertedBounds() const;
const std::unique_ptr<Gurobi> gurobi_;
// Note that we use linked_hash_map because the index of the gurobi_model_

42
ortools/math_opt/solvers/pdlp_bridge.cc
View File

@@ -24,6 +24,7 @@
#include "absl/status/status.h"
#include "absl/status/statusor.h"
#include "absl/strings/str_cat.h"
#include "ortools/math_opt/core/inverted_bounds.h"
#include "ortools/math_opt/core/math_opt_proto_utils.h"
#include "ortools/math_opt/core/sparse_vector_view.h"
#include "ortools/math_opt/model.pb.h"
@@ -103,26 +104,29 @@ absl::StatusOr<PdlpBridge> PdlpBridge::FromProto(
}
const SparseDoubleMatrixProto& quadratic_objective =
model_proto.objective().quadratic_coefficients();
std::vector<Eigen::Triplet<double, int64_t>> obj_triplets;
const int obj_nnz = quadratic_objective.row_ids().size();
if (obj_nnz > 0) {
pdlp_lp.objective_matrix.emplace();
pdlp_lp.objective_matrix->setZero(variables.ids_size());
}
for (int i = 0; i < obj_nnz; ++i) {
const int64_t row_index =
result.var_id_to_pdlp_index_.at(quadratic_objective.row_ids(i));
const int64_t column_index =
result.var_id_to_pdlp_index_.at(quadratic_objective.column_ids(i));
const double value = obj_scale * quadratic_objective.coefficients(i);
if (row_index != column_index) {
return absl::InvalidArgumentError(
"PDLP cannot solve problems with non-diagonal objective matrices");
}
// MathOpt represents quadratic objectives in "terms" form, i.e. as a sum
// of double * Variable * Variable terms. They are stored in upper
// triangular form with row_index <= column_index. In contrast, PDLP
// represents quadratic objectives in "matrix" form as 1/2 x'Qx; it wants
// the symmetric matrix Q. To get to the right format, we simply add each
// term and its transposed entry, and defer to Eigen to sum duplicate
// entries along the diagonal.
obj_triplets.push_back({row_index, column_index, value});
obj_triplets.push_back({column_index, row_index, value});
// represents quadratic objectives in "matrix" form as 1/2 x'Qx, where Q is
// diagonal. To get to the right format, we simply double each diagonal
// entry.
pdlp_lp.objective_matrix->diagonal()[row_index] = 2 * value;
}
pdlp_lp.objective_matrix.setFromTriplets(obj_triplets.begin(),
obj_triplets.end());
pdlp_lp.objective_scaling_factor = obj_scale;
// Note: MathOpt stores the constraint data in row major order, but PDLP
// wants the data in column major order. There is probably a more efficient
@@ -145,6 +149,26 @@ absl::StatusOr<PdlpBridge> PdlpBridge::FromProto(
return result;
}
InvertedBounds PdlpBridge::ListInvertedBounds() const {
InvertedBounds inverted_bounds;
for (int64_t var_index = 0; var_index < pdlp_index_to_var_id_.size();
++var_index) {
if (pdlp_lp_.variable_lower_bounds[var_index] >
pdlp_lp_.variable_upper_bounds[var_index]) {
inverted_bounds.variables.push_back(pdlp_index_to_var_id_[var_index]);
}
}
for (int64_t lin_con_index = 0;
lin_con_index < pdlp_index_to_lin_con_id_.size(); ++lin_con_index) {
if (pdlp_lp_.constraint_lower_bounds[lin_con_index] >
pdlp_lp_.constraint_upper_bounds[lin_con_index]) {
inverted_bounds.linear_constraints.push_back(
pdlp_index_to_lin_con_id_[lin_con_index]);
}
}
return inverted_bounds;
}
absl::StatusOr<SparseDoubleVectorProto> PdlpBridge::PrimalVariablesToProto(
const Eigen::VectorXd& primal_values,
const SparseVectorFilterProto& variable_filter) const {

4
ortools/math_opt/solvers/pdlp_bridge.h
View File

@@ -20,6 +20,7 @@
#include "Eigen/Core"
#include "absl/container/flat_hash_map.h"
#include "absl/status/statusor.h"
#include "ortools/math_opt/core/inverted_bounds.h"
#include "ortools/math_opt/model.pb.h"
#include "ortools/math_opt/sparse_containers.pb.h"
#include "ortools/pdlp/quadratic_program.h"
@@ -48,6 +49,9 @@ class PdlpBridge {
const pdlp::QuadraticProgram& pdlp_lp() const { return pdlp_lp_; }
// Returns the ids of variables and linear constraints with inverted bounds.
InvertedBounds ListInvertedBounds() const;
// TODO(b/183616124): we need to support the inverse of these methods for
// warm start.
absl::StatusOr<SparseDoubleVectorProto> PrimalVariablesToProto(

27
ortools/math_opt/solvers/pdlp_solver.cc
View File

@@ -13,6 +13,7 @@
#include "ortools/math_opt/solvers/pdlp_solver.h"
#include <algorithm>
#include <cmath>
#include <cstdint>
#include <functional>
@@ -34,6 +35,7 @@
#include "ortools/base/protoutil.h"
#include "ortools/base/status_macros.h"
#include "ortools/math_opt/callback.pb.h"
#include "ortools/math_opt/core/inverted_bounds.h"
#include "ortools/math_opt/core/math_opt_proto_utils.h"
#include "ortools/math_opt/core/solve_interrupter.h"
#include "ortools/math_opt/core/solver_interface.h"
@@ -71,9 +73,7 @@ absl::StatusOr<PrimalDualHybridGradientParams> PdlpSolver::MergeParameters(
PrimalDualHybridGradientParams result;
std::vector<std::string> warnings;
if (parameters.enable_output()) {
// TODO(b/183502493): this is not a robust solution. It is not thread safe
// and will interfere with the global vlog state.
SetVLOGLevel("primal_dual_hybrid_gradient", 2);
result.set_verbosity_level(3);
}
if (parameters.has_threads()) {
result.set_num_threads(parameters.threads());
@@ -83,6 +83,9 @@ absl::StatusOr<PrimalDualHybridGradientParams> PdlpSolver::MergeParameters(
absl::ToDoubleSeconds(
util_time::DecodeGoogleApiProto(parameters.time_limit()).value()));
}
if (parameters.has_node_limit()) {
warnings.push_back("parameter node_limit not supported for PDLP");
}
if (parameters.has_cutoff_limit()) {
warnings.push_back("parameter cutoff_limit not supported for PDLP");
}
@@ -149,6 +152,8 @@ absl::StatusOr<TerminationProto> ConvertReason(
case pdlp::TERMINATION_REASON_NUMERICAL_ERROR:
return TerminateForReason(TERMINATION_REASON_NUMERICAL_ERROR,
pdlp_detail);
case pdlp::TERMINATION_REASON_INTERRUPTED_BY_USER:
return NoSolutionFoundTermination(LIMIT_INTERRUPTED, pdlp_detail);
case pdlp::TERMINATION_REASON_INVALID_PROBLEM:
// Indicates that the solver detected invalid problem data, e.g.,
// inconsistent bounds.
@@ -240,7 +245,8 @@ absl::StatusOr<SolveResultProto> PdlpSolver::MakeSolveResult(
case pdlp::TERMINATION_REASON_TIME_LIMIT:
case pdlp::TERMINATION_REASON_ITERATION_LIMIT:
case pdlp::TERMINATION_REASON_KKT_MATRIX_PASS_LIMIT:
case pdlp::TERMINATION_REASON_NUMERICAL_ERROR: {
case pdlp::TERMINATION_REASON_NUMERICAL_ERROR:
case pdlp::TERMINATION_REASON_INTERRUPTED_BY_USER: {
SolutionProto* solution_proto = result.add_solutions();
{
auto maybe_primal = pdlp_bridge_.PrimalVariablesToProto(
@@ -334,8 +340,6 @@ absl::StatusOr<SolveResultProto> PdlpSolver::Solve(
const MessageCallback message_cb,
const CallbackRegistrationProto& callback_registration, const Callback cb,
SolveInterrupter* const interrupter) {
// TODO(b/192274409): Use interrupter if PDLP supports interruption.
// TODO(b/183502493): Implement message callback when PDLP supports that.
if (message_cb != nullptr) {
return absl::InvalidArgumentError(internal::kMessageCallbackNotSupported);
@@ -345,8 +349,17 @@ absl::StatusOr<SolveResultProto> PdlpSolver::Solve(
/*supported_events=*/{}));
ASSIGN_OR_RETURN(auto pdlp_params, MergeParameters(parameters));
// PDLP returns `(TERMINATION_REASON_INVALID_PROBLEM): The input problem has
// inconsistent bounds.` but we want a more detailed error.
RETURN_IF_ERROR(pdlp_bridge_.ListInvertedBounds().ToStatus());
std::atomic<bool> interrupt = false;
const ScopedSolveInterrupterCallback set_interrupt(
interrupter, [&]() { interrupt = true; });
const SolverResult pdlp_result =
PrimalDualHybridGradient(pdlp_bridge_.pdlp_lp(), pdlp_params);
PrimalDualHybridGradient(pdlp_bridge_.pdlp_lp(), pdlp_params, &interrupt);
return MakeSolveResult(pdlp_result, model_parameters);
}
absl::Status PdlpSolver::Update(const ModelUpdateProto& model_update) {

24
ortools/math_opt/tools/BUILD.bazel Normal file
View File

@@ -0,0 +1,24 @@
package(default_visibility = ["//visibility:private"])
cc_binary(
name = "mathopt_solve",
srcs = ["mathopt_solve_main.cc"],
deps = [
"//ortools/base",
"//ortools/base:status_macros",
"//ortools/math_opt:parameters_cc_proto",
"//ortools/math_opt/core:solver_interface",
"//ortools/math_opt/cpp:math_opt",
"//ortools/math_opt/io:mps_converter",
"//ortools/math_opt/solvers:cp_sat_solver",
"//ortools/math_opt/solvers:glop_solver",
"//ortools/math_opt/solvers:glpk_solver",
"//ortools/math_opt/solvers:gscip_solver",
"//ortools/math_opt/solvers:gurobi_solver",
"@com_google_absl//absl/flags:flag",
"@com_google_absl//absl/status",
"@com_google_absl//absl/status:statusor",
"@com_google_absl//absl/strings",
"@com_google_absl//absl/time",
],
)

266
ortools/math_opt/tools/mathopt_solve_main.cc Normal file
View File

@@ -0,0 +1,266 @@
// Copyright 2010-2021 Google LLC
// Licensed under the Apache License, Version 2.0 (the "License");
// you may not use this file except in compliance with the License.
// You may obtain a copy of the License at
//
// http://www.apache.org/licenses/LICENSE-2.0
//
// Unless required by applicable law or agreed to in writing, software
// distributed under the License is distributed on an "AS IS" BASIS,
// WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
// See the License for the specific language governing permissions and
// limitations under the License.
// Tool to run MathOpt on the given problems.
//
// Examples:
//
// mathopt_solve --input_file model.pb
//
// mathopt_solve --input_file model.mps.gz --solver_type=glop
//
// mathopt_solve --input_file model --solver_logs --format=mathopt
//
#include <iostream>
#include <memory>
#include <optional>
#include <string>
#include <utility>
#include <vector>
#include "absl/flags/flag.h"
#include "absl/status/status.h"
#include "absl/status/statusor.h"
#include "absl/strings/match.h"
#include "absl/strings/str_cat.h"
#include "absl/strings/str_join.h"
#include "absl/strings/string_view.h"
#include "absl/time/time.h"
#include "ortools/base/file.h"
#include "ortools/base/init_google.h"
#include "ortools/base/logging.h"
#include "ortools/base/status_macros.h"
#include "ortools/math_opt/core/solver_interface.h"
#include "ortools/math_opt/cpp/math_opt.h"
#include "ortools/math_opt/io/mps_converter.h"
#include "ortools/math_opt/parameters.pb.h"
#include "ortools/util/status_macros.h"
inline constexpr absl::string_view kMathOptBinaryFormat = "mathopt";
inline constexpr absl::string_view kMathOptTextFormat = "mathopt_txt";
inline constexpr absl::string_view kMPSFormat = "mps";
inline constexpr absl::string_view kAutoFormat = "auto";
inline constexpr absl::string_view kPbExt = ".pb";
inline constexpr absl::string_view kProtoExt = ".proto";
inline constexpr absl::string_view kPbTxtExt = ".pb.txt";
inline constexpr absl::string_view kTextProtoExt = ".textproto";
inline constexpr absl::string_view kMPSExt = ".mps";
inline constexpr absl::string_view kMPSGzipExt = ".mps.gz";
namespace {
struct SolverTypeProtoFormatter {
void operator()(
std::string* const out,
const operations_research::math_opt::SolverTypeProto solver_type) {
out->append(EnumToString(EnumFromProto(solver_type).value()));
}
};
} // namespace
ABSL_FLAG(std::string, input_file, "",
"the file containing the model to solve; use --format to specify the "
"file format");
ABSL_FLAG(std::string, format, "auto",
absl::StrCat(
"the format of the --input_file; possible values:\n", "* ",
kMathOptBinaryFormat, ": for a MathOpt ModelProto in binary\n",
"* ", kMathOptTextFormat, ": when the proto is in text\n", "* ",
kMPSFormat, ": for MPS file (which can be GZiped)\n", "* ",
kAutoFormat, ": to guess the format from the file extension:\n",
" - '", kPbExt, "', '", kProtoExt, "': ", kMathOptBinaryFormat,
"\n", " - '", kPbTxtExt, "', '", kTextProtoExt,
"': ", kMathOptTextFormat, "\n", " - '", kMPSExt, "', '",
kMPSGzipExt, "': ", kMPSFormat));
ABSL_FLAG(
std::vector<std::string>, update_files, {},
absl::StrCat(
"the file containing ModelUpdateProto to apply to the --input_file; "
"when this flag is used, the --format must be either ",
kMathOptBinaryFormat, " or ", kMathOptTextFormat));
ABSL_FLAG(operations_research::math_opt::SolverType, solver_type,
operations_research::math_opt::SolverType::kGscip,
absl::StrCat(
"the solver to use, possible values: ",
absl::StrJoin(
operations_research::math_opt::AllSolversRegistry::Instance()
->RegisteredSolvers(),
", ", SolverTypeProtoFormatter())));
ABSL_FLAG(bool, solver_logs, false,
"use a message callback to print the solver convergence logs");
ABSL_FLAG(absl::Duration, time_limit, absl::InfiniteDuration(),
"the time limit to use for the solve");
namespace operations_research {
namespace math_opt {
namespace {
// Returned the guessed format (one of the kXxxFormat constant) from the file
// extension; or nullopt.
std::optional<absl::string_view> FormatFromFilePath(
const absl::string_view file_path) {
const std::vector<std::pair<absl::string_view, absl::string_view>>
extension_to_format = {
{kPbExt, kMathOptBinaryFormat}, {kProtoExt, kMathOptBinaryFormat},
{kPbTxtExt, kMathOptTextFormat}, {kTextProtoExt, kMathOptTextFormat},
{kMPSExt, kMPSFormat}, {kMPSGzipExt, kMPSFormat},
};
for (const auto& [ext, format] : extension_to_format) {
if (absl::EndsWith(file_path, ext)) {
return format;
}
}
return std::nullopt;
}
// Returns the ModelProto read from the given file. The format must not be
// kAutoFormat; other invalid values will be reported as QFATAL log mentioning
// the --format flag.
absl::StatusOr<ModelProto> ReadModel(const absl::string_view file_path,
const absl::string_view format) {
if (format == kMathOptBinaryFormat) {
return file::GetBinaryProto<ModelProto>(file_path, file::Defaults());
}
if (format == kMathOptTextFormat) {
return file::GetTextProto<ModelProto>(file_path, file::Defaults());
}
if (format == kMPSFormat) {
return ReadMpsFile(file_path);
}
LOG(QFATAL) << "Unsupported value of --format: " << format;
}
// Returns the ModelUpdateProto read from the given file. The format must be
// kMathOptBinaryFormat or kMathOptTextFormat; other values will generate an
// error.
absl::StatusOr<ModelUpdateProto> ReadModelUpdate(
const absl::string_view file_path, const absl::string_view format) {
if (format == kMathOptBinaryFormat) {
return file::GetBinaryProto<ModelUpdateProto>(file_path, file::Defaults());
}
if (format == kMathOptTextFormat) {
return file::GetTextProto<ModelUpdateProto>(file_path, file::Defaults());
}
return absl::InternalError(
absl::StrCat("invalid format in ReadModelUpdate(): ", format));
}
// Prints the summary of the solve result.
absl::Status PrintSummary(const SolveResult& result) {
std::cout << "Solve finished:\n"
<< " termination: " << result.termination << "\n"
<< " solve time: " << result.solve_stats.solve_time
<< "\n best primal bound: " << result.solve_stats.best_primal_bound
<< "\n best dual bound: " << result.solve_stats.best_dual_bound
<< std::endl;
if (result.solutions.empty()) {
std::cout << " no solution" << std::endl;
}
for (int i = 0; i < result.solutions.size(); ++i) {
const Solution& solution = result.solutions[i];
std::cout << " solution #" << (i + 1) << " objective: ";
if (solution.primal_solution.has_value()) {
std::cout << solution.primal_solution->objective_value;
} else {
std::cout << "n/a";
}
std::cout << std::endl;
}
return absl::OkStatus();
}
absl::Status RunSolver() {
const std::string input_file_path = absl::GetFlag(FLAGS_input_file);
if (input_file_path.empty()) {
LOG(QFATAL) << "The flag --input_file is mandatory.";
}
// Parses --format.
std::string format = absl::GetFlag(FLAGS_format);
if (format == kAutoFormat) {
const std::optional<absl::string_view> guessed_format =
FormatFromFilePath(input_file_path);
if (!guessed_format) {
LOG(QFATAL) << "Can't guess the format from the file extension, please "
"use --format to specify the file format explicitly.";
}
format = *guessed_format;
}
// We deal with input validation in the ReadModel() function.
// Read the model and the optional updates.
const std::vector<std::string> update_file_paths =
absl::GetFlag(FLAGS_update_files);
if (!update_file_paths.empty() && format != kMathOptBinaryFormat &&
format != kMathOptTextFormat) {
LOG(QFATAL) << "Can't use --update_files with a input of format " << format
<< ".";
}
OR_ASSIGN_OR_RETURN3(const ModelProto model_proto,
ReadModel(input_file_path, format),
_ << "failed to read " << input_file_path);
std::vector<ModelUpdateProto> model_updates;
for (const std::string& update_file_path : update_file_paths) {
ASSIGN_OR_RETURN(ModelUpdateProto update,
ReadModelUpdate(update_file_path, format));
model_updates.emplace_back(std::move(update));
}
// Solve the problem.
ASSIGN_OR_RETURN(const std::unique_ptr<Model> model,
Model::FromModelProto(model_proto));
for (int u = 0; u < model_updates.size(); ++u) {
const ModelUpdateProto& update = model_updates[u];
RETURN_IF_ERROR(model->ApplyUpdateProto(update))
<< "failed to apply the update file: " << update_file_paths[u];
}
SolveArguments solve_args = {
.parameters = {.time_limit = absl::GetFlag(FLAGS_time_limit)},
};
if (absl::GetFlag(FLAGS_solver_logs)) {
solve_args.message_callback = PrinterMessageCallback(std::cout, "logs| ");
}
OR_ASSIGN_OR_RETURN3(
const SolveResult result,
Solve(*model, absl::GetFlag(FLAGS_solver_type), solve_args),
_ << "the solver failed");
RETURN_IF_ERROR(PrintSummary(result));
return absl::OkStatus();
}
} // namespace
} // namespace math_opt
} // namespace operations_research
int main(int argc, char* argv[]) {
InitGoogle(argv[0], &argc, &argv, /*remove_flags=*/true);
const absl::Status status = operations_research::math_opt::RunSolver();
// We don't use QCHECK_OK() here since the logged message contains more than
// the failing status.
if (!status.ok()) {
LOG(QFATAL) << status;
}
return 0;
}

7
ortools/math_opt/validators/BUILD.bazel
View File

@@ -101,8 +101,8 @@ cc_library(
deps = [
":scalar_validator",
"//ortools/base",
"//ortools/base:status_macros",
"//ortools/base:protoutil",
"//ortools/base:status_macros",
"//ortools/math_opt:result_cc_proto",
"//ortools/port:proto_utils",
"@com_google_absl//absl/status",
@@ -161,6 +161,7 @@ cc_library(
"//ortools/base:protoutil",
"//ortools/base:status_macros",
"//ortools/math_opt:parameters_cc_proto",
"//ortools/util:status_macros",
"@com_google_absl//absl/memory",
"@com_google_absl//absl/status",
"@com_google_absl//absl/status:statusor",
@@ -180,9 +181,9 @@ cc_library(
"//ortools/base",
"//ortools/base:status_macros",
"//ortools/math_opt:callback_cc_proto",
"//ortools/math_opt/core:model_summary",
"//ortools/math_opt:solution_cc_proto",
"//ortools/math_opt:sparse_containers_cc_proto",
"//ortools/math_opt/core:model_summary",
"//ortools/math_opt/core:sparse_vector_view",
"//ortools/port:proto_utils",
"@com_google_absl//absl/status",
@@ -200,8 +201,8 @@ cc_library(
"//ortools/base",
"//ortools/base:status_macros",
"//ortools/math_opt:model_parameters_cc_proto",
"//ortools/math_opt/core:model_summary",
"//ortools/math_opt:sparse_containers_cc_proto",
"//ortools/math_opt/core:model_summary",
"//ortools/math_opt/validators:solution_validator",
"@com_google_absl//absl/status",
"@com_google_absl//absl/strings",

2
ortools/math_opt/validators/ids_validator.cc
View File

@@ -33,8 +33,6 @@
namespace operations_research {
namespace math_opt {
constexpr double kInf = std::numeric_limits<double>::infinity();
namespace {
absl::Status CheckSortedIdsSubsetWithIndexOffset(
const absl::Span<const int64_t> ids,

2
ortools/math_opt/validators/model_validator.cc
View File

@@ -38,8 +38,6 @@ namespace operations_research {
namespace math_opt {
namespace {
constexpr double kInf = std::numeric_limits<double>::infinity();
////////////////////////////////////////////////////////////////////////////////
// Submessages
////////////////////////////////////////////////////////////////////////////////

12
ortools/math_opt/validators/result_validator.cc
View File

@@ -31,18 +31,6 @@ namespace operations_research {
namespace math_opt {
namespace {
constexpr double kInf = std::numeric_limits<double>::infinity();
absl::Status ValidateSolutionStatus(const SolutionStatusProto& status) {
if (!SolutionStatusProto_IsValid(status)) {
return absl::InvalidArgumentError(absl::StrCat("status = ", status));
}
if (status == SOLUTION_STATUS_UNSPECIFIED) {
return absl::InvalidArgumentError("status = SOLUTION_STATUS_UNSPECIFIED");
}
return absl::OkStatus();
}
absl::Status ValidateTermination(const TerminationProto& termination) {
if (termination.reason() == TERMINATION_REASON_UNSPECIFIED) {
return absl::InvalidArgumentError("termination reason must be specified");

2
ortools/math_opt/validators/solution_validator.cc
View File

@@ -34,8 +34,6 @@ namespace operations_research {
namespace math_opt {
namespace {
constexpr double kInf = std::numeric_limits<double>::infinity();
absl::Status ValidateSolutionStatus(const SolutionStatusProto& status) {
if (!SolutionStatusProto_IsValid(status)) {
return absl::InvalidArgumentError(absl::StrCat("status = ", status));

41
ortools/math_opt/validators/solve_parameters_validator.cc
View File

@@ -13,19 +13,33 @@
#include "ortools/math_opt/validators/solve_parameters_validator.h"
#include <cmath>
#include <string>
#include <type_traits>
#include "absl/status/status.h"
#include "absl/status/statusor.h"
#include "absl/strings/str_cat.h"
#include "absl/time/time.h"
#include "ortools/base/protoutil.h"
#include "ortools/base/status_macros.h"
#include "ortools/math_opt/parameters.pb.h"
#include "ortools/util/status_macros.h"
namespace operations_research {
namespace math_opt {
absl::Status ValidateSolveParameters(const SolveParametersProto& parameters) {
RETURN_IF_ERROR(
util_time::DecodeGoogleApiProto(parameters.time_limit()).status())
<< "invalid SolveParameters.time_limit";
{
OR_ASSIGN_OR_RETURN3(
const absl::Duration time_limit,
util_time::DecodeGoogleApiProto(parameters.time_limit()),
_ << "invalid SolveParameters.time_limit");
if (time_limit < absl::ZeroDuration()) {
return util::InvalidArgumentErrorBuilder()
<< "SolveParameters.time_limit = " << time_limit << " < 0";
}
}
if (parameters.has_threads()) {
if (parameters.threads() <= 0) {
@@ -34,21 +48,26 @@ absl::Status ValidateSolveParameters(const SolveParametersProto& parameters) {
}
}
if (parameters.has_relative_gap_limit()) {
if (parameters.relative_gap_limit() < 0) {
if (parameters.has_relative_gap_tolerance()) {
if (parameters.relative_gap_tolerance() < 0) {
return absl::InvalidArgumentError(
absl::StrCat("SolveParameters.relative_gap_limit = ",
parameters.relative_gap_limit(), " < 0"));
absl::StrCat("SolveParameters.relative_gap_tolerance = ",
parameters.relative_gap_tolerance(), " < 0"));
}
}
if (parameters.has_absolute_gap_limit()) {
if (parameters.absolute_gap_limit() < 0) {
if (parameters.has_absolute_gap_tolerance()) {
if (parameters.absolute_gap_tolerance() < 0) {
return absl::InvalidArgumentError(
absl::StrCat("SolveParameters.absolute_gap_limit = ",
parameters.absolute_gap_limit(), " < 0"));
absl::StrCat("SolveParameters.absolute_gap_tolerance = ",
parameters.absolute_gap_tolerance(), " < 0"));
}
}
if (parameters.has_node_limit() && parameters.node_limit() < 0) {
return util::InvalidArgumentErrorBuilder()
<< "SolveParameters.node_limit = " << parameters.node_limit()
<< " should be nonnegative.";
}
if (parameters.has_solution_limit() && parameters.solution_limit() <= 0) {
return util::InvalidArgumentErrorBuilder()

2
ortools/math_opt/validators/solve_stats_validator.cc
View File

@@ -31,8 +31,6 @@ namespace operations_research {
namespace math_opt {
namespace {
constexpr double kInf = std::numeric_limits<double>::infinity();
absl::Status ValidateFeasibilityStatus(const FeasibilityStatusProto& status) {
if (!FeasibilityStatusProto_IsValid(status)) {
return absl::InvalidArgumentError(absl::StrCat("invalid status ", status));