diff --git a/ortools/math_opt/core/BUILD.bazel b/ortools/math_opt/core/BUILD.bazel index 4f6b8b0c9f..4ce3175cda 100644 --- a/ortools/math_opt/core/BUILD.bazel +++ b/ortools/math_opt/core/BUILD.bazel @@ -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", + ], +) diff --git a/ortools/math_opt/core/inverted_bounds.cc b/ortools/math_opt/core/inverted_bounds.cc new file mode 100644 index 0000000000..3a5cbdf254 --- /dev/null +++ b/ortools/math_opt/core/inverted_bounds.cc @@ -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 +#include +#include +#include + +#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& 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 diff --git a/ortools/math_opt/core/inverted_bounds.h b/ortools/math_opt/core/inverted_bounds.h new file mode 100644 index 0000000000..ef8eeec504 --- /dev/null +++ b/ortools/math_opt/core/inverted_bounds.h @@ -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 +#include + +#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 variables; + + // Ids of the linear constraints with inverted bounds. + std::vector linear_constraints; +}; + +} // namespace operations_research::math_opt + +#endif // OR_TOOLS_MATH_OPT_CORE_INVERTED_BOUNDS_H_ diff --git a/ortools/math_opt/core/math_opt_proto_utils.h b/ortools/math_opt/core/math_opt_proto_utils.h index 332cde8532..d4bc118673 100644 --- a/ortools/math_opt/core/math_opt_proto_utils.h +++ b/ortools/math_opt/core/math_opt_proto_utils.h @@ -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 }; diff --git a/ortools/math_opt/core/model_storage.cc b/ortools/math_opt/core/model_storage.cc index 0c6c95e375..1b2d534ae8 100644 --- a/ortools/math_opt/core/model_storage.cc +++ b/ortools/math_opt/core/model_storage.cc @@ -748,7 +748,7 @@ std::optional ModelStorage::ExportModelUpdate( // * Use n-way merge here if the performances justify it. const auto merge = std::make_shared(); 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 ModelStorage::ExportModelUpdate( const std::optional 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 diff --git a/ortools/math_opt/core/model_storage.h b/ortools/math_opt/core/model_storage.h index 20b46fc84b..271f6c62d9 100644 --- a/ortools/math_opt/core/model_storage.h +++ b/ortools/math_opt/core/model_storage.h @@ -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: diff --git a/ortools/math_opt/core/model_update_merge.h b/ortools/math_opt/core/model_update_merge.h index 564e615539..c5385d079d 100644 --- a/ortools/math_opt/core/model_update_merge.h +++ b/ortools/math_opt/core/model_update_merge.h @@ -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; diff --git a/ortools/math_opt/core/non_streamable_solver_init_arguments.h b/ortools/math_opt/core/non_streamable_solver_init_arguments.h index 58f1344877..2b2afd69ce 100644 --- a/ortools/math_opt/core/non_streamable_solver_init_arguments.h +++ b/ortools/math_opt/core/non_streamable_solver_init_arguments.h @@ -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 diff --git a/ortools/math_opt/cpp/BUILD.bazel b/ortools/math_opt/cpp/BUILD.bazel index 74ea804bfc..c2af36300c 100644 --- a/ortools/math_opt/cpp/BUILD.bazel +++ b/ortools/math_opt/cpp/BUILD.bazel @@ -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", ], ) diff --git a/ortools/math_opt/cpp/id_map.h b/ortools/math_opt/cpp/id_map.h index 34ef15c35d..9205293c12 100644 --- a/ortools/math_opt/cpp/id_map.h +++ b/ortools/math_opt/cpp/id_map.h @@ -163,7 +163,7 @@ class IdMap { inline std::pair 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::iterator, bool> IdMap::try_emplace( } template -int IdMap::erase(const K& k) { +typename IdMap::size_type IdMap::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; } diff --git a/ortools/math_opt/cpp/id_set.h b/ortools/math_opt/cpp/id_set.h index 96f83e2e63..5480cf180c 100644 --- a/ortools/math_opt/cpp/id_set.h +++ b/ortools/math_opt/cpp/id_set.h @@ -125,7 +125,7 @@ class IdSet { inline std::pair 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::const_iterator, bool> IdSet::emplace( } template -int IdSet::erase(const K& k) { +typename IdSet::size_type IdSet::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; } diff --git a/ortools/math_opt/cpp/matchers.cc b/ortools/math_opt/cpp/matchers.cc index dbf20a9fab..1d7f34e0fc 100644 --- a/ortools/math_opt/cpp/matchers.cc +++ b/ortools/math_opt/cpp/matchers.cc @@ -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; diff --git a/ortools/math_opt/cpp/model.cc b/ortools/math_opt/cpp/model.cc index 1afca1ff19..5990e5db5d 100644 --- a/ortools/math_opt/cpp/model.cc +++ b/ortools/math_opt/cpp/model.cc @@ -209,7 +209,7 @@ QuadraticExpression Model::ObjectiveAsQuadraticExpression() const { ModelProto Model::ExportModel() const { return storage()->ExportModel(); } -std::unique_ptr Model::NewUpdateTracker() const { +std::unique_ptr Model::NewUpdateTracker() { return std::make_unique(storage_); } diff --git a/ortools/math_opt/cpp/model.h b/ortools/math_opt/cpp/model.h index 84a6fc68c3..c0be707646 100644 --- a/ortools/math_opt/cpp/model.h +++ b/ortools/math_opt/cpp/model.h @@ -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::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 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 NewUpdateTracker() const; + std::unique_ptr 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 { diff --git a/ortools/math_opt/cpp/parameters.cc b/ortools/math_opt/cpp/parameters.cc index 32453b3afc..2c4492b61b 100644 --- a/ortools/math_opt/cpp/parameters.cc +++ b/ortools/math_opt/cpp/parameters.cc @@ -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::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::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()); diff --git a/ortools/math_opt/cpp/parameters.h b/ortools/math_opt/cpp/parameters.h index c84de098f6..ca3f3ef2a8 100644 --- a/ortools/math_opt/cpp/parameters.h +++ b/ortools/math_opt/cpp/parameters.h @@ -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 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 relative_gap_limit; - std::optional 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 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 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 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 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 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 relative_gap_tolerance; + // The algorithm for solving a linear program. If nullopt, use the solver // default algorithm. // diff --git a/ortools/math_opt/cpp/solve.cc b/ortools/math_opt/cpp/solve.cc index 60a5e8347b..0e87f60c2c 100644 --- a/ortools/math_opt/cpp/solve.cc +++ b/ortools/math_opt/cpp/solve.cc @@ -102,14 +102,18 @@ absl::StatusOr Solve(const Model& model, } absl::StatusOr> IncrementalSolver::New( - Model& model, const SolverType solver_type, SolverInitArguments arguments) { - std::unique_ptr update_tracker = model.NewUpdateTracker(); - ASSIGN_OR_RETURN( - std::unique_ptr 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 update_tracker = model->NewUpdateTracker(); + ASSIGN_OR_RETURN(const ModelProto model_proto, update_tracker->ExportModel()); + ASSIGN_OR_RETURN(std::unique_ptr solver, + Solver::New(EnumToProto(solver_type), model_proto, + ToSolverInitArgs(arguments))); return absl::WrapUnique( - 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 IncrementalSolver::Solve( } absl::StatusOr IncrementalSolver::Update() { - std::optional model_update = - update_tracker_->ExportModelUpdate(); + ASSIGN_OR_RETURN(std::optional 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)); diff --git a/ortools/math_opt/cpp/solve.h b/ortools/math_opt/cpp/solve.h index 83daf557a4..b5dbea1270 100644 --- a/ortools/math_opt/cpp/solve.h +++ b/ortools/math_opt/cpp/solve.h @@ -75,15 +75,15 @@ absl::StatusOr 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 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 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> 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. // diff --git a/ortools/math_opt/cpp/streamable_solver_init_arguments.cc b/ortools/math_opt/cpp/streamable_solver_init_arguments.cc index 0e03bc59ca..92166b6d10 100644 --- a/ortools/math_opt/cpp/streamable_solver_init_arguments.cc +++ b/ortools/math_opt/cpp/streamable_solver_init_arguments.cc @@ -16,6 +16,7 @@ #include #include +#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::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 diff --git a/ortools/math_opt/cpp/streamable_solver_init_arguments.h b/ortools/math_opt/cpp/streamable_solver_init_arguments.h index e6f4510b10..4b78e5d2c4 100644 --- a/ortools/math_opt/cpp/streamable_solver_init_arguments.h +++ b/ortools/math_opt/cpp/streamable_solver_init_arguments.h @@ -24,7 +24,7 @@ #include #include -#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 FromProto( + const SolverInitializerProto& args_proto); }; } // namespace math_opt diff --git a/ortools/math_opt/cpp/update_tracker.cc b/ortools/math_opt/cpp/update_tracker.cc index 1a60022ec0..ec7f6275b2 100644 --- a/ortools/math_opt/cpp/update_tracker.cc +++ b/ortools/math_opt/cpp/update_tracker.cc @@ -16,6 +16,8 @@ #include #include +#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& storage) : storage_(ABSL_DIE_IF_NULL(storage)), update_tracker_(storage->NewUpdateTracker()) {} -std::optional UpdateTracker::ExportModelUpdate() { +absl::StatusOr> +UpdateTracker::ExportModelUpdate() { const std::shared_ptr 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 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 UpdateTracker::ExportModel() const { const std::shared_ptr storage = storage_.lock(); - CHECK(storage != nullptr) << internal::kModelIsDestroyed; + if (storage == nullptr) { + return absl::InvalidArgumentError(internal::kModelIsDestroyed); + } return storage->ExportModel(); } diff --git a/ortools/math_opt/cpp/update_tracker.h b/ortools/math_opt/cpp/update_tracker.h index 59bab1fff0..90dcace59c 100644 --- a/ortools/math_opt/cpp/update_tracker.h +++ b/ortools/math_opt/cpp/update_tracker.h @@ -17,6 +17,8 @@ #include #include +#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 update_proto = -// update_tracker.ExportModelUpdate(); -// update_tracker.Checkpoint(); +// ASSIGN_OR_RETURN(const std::optional 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 ExportModelUpdate(); + // + // If fails if the Model has been destroyed. + absl::StatusOr> 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 ExportModel() const; private: const std::weak_ptr 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 diff --git a/ortools/math_opt/io/BUILD.bazel b/ortools/math_opt/io/BUILD.bazel index aeef34897d..2d18725b74 100644 --- a/ortools/math_opt/io/BUILD.bazel +++ b/ortools/math_opt/io/BUILD.bazel @@ -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", + ], +) diff --git a/ortools/math_opt/io/mps_converter.cc b/ortools/math_opt/io/mps_converter.cc new file mode 100644 index 0000000000..d28ef83abc --- /dev/null +++ b/ortools/math_opt/io/mps_converter.cc @@ -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 + +#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 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 diff --git a/ortools/math_opt/io/mps_converter.h b/ortools/math_opt/io/mps_converter.h new file mode 100644 index 0000000000..ccdca4d04c --- /dev/null +++ b/ortools/math_opt/io/mps_converter.h @@ -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 + +#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 ReadMpsFile(absl::string_view filename); + +} // namespace operations_research::math_opt + +#endif // OR_TOOLS_MATH_OPT_IO_MPS_CONVERTER_H_ diff --git a/ortools/math_opt/parameters.proto b/ortools/math_opt/parameters.proto index f16616b7ff..fa72e37fd8 100644 --- a/ortools/math_opt/parameters.proto +++ b/ortools/math_opt/parameters.proto @@ -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. // diff --git a/ortools/math_opt/result.proto b/ortools/math_opt/result.proto index b553301b8a..016346f87c 100644 --- a/ortools/math_opt/result.proto +++ b/ortools/math_opt/result.proto @@ -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 diff --git a/ortools/math_opt/samples/basic_example.cc b/ortools/math_opt/samples/basic_example.cc index 7e2a010a96..34b0585873 100644 --- a/ortools/math_opt/samples/basic_example.cc +++ b/ortools/math_opt/samples/basic_example.cc @@ -16,82 +16,56 @@ #include #include +#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::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; } diff --git a/ortools/math_opt/samples/cocktail_hour.cc b/ortools/math_opt/samples/cocktail_hour.cc index eb64d05064..ce524647e9 100644 --- a/ortools/math_opt/samples/cocktail_hour.cc +++ b/ortools/math_opt/samples/cocktail_hour.cc @@ -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 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& 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({"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 = + 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; } diff --git a/ortools/math_opt/samples/cutting_stock.cc b/ortools/math_opt/samples/cutting_stock.cc index ac30760d4c..824f303d48 100644 --- a/ortools/math_opt/samples/cutting_stock.cc +++ b/ortools/math_opt/samples/cutting_stock.cc @@ -137,8 +137,8 @@ absl::StatusOr> 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 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 prices; for (const math_opt::LinearConstraint d : demand_met) { prices.push_back(solve_result.dual_values().at(d)); @@ -219,11 +227,14 @@ absl::StatusOr 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) { diff --git a/ortools/math_opt/samples/facility_lp_benders.cc b/ortools/math_opt/samples/facility_lp_benders.cc index cd13dcd268..49b63eb4ad 100644 --- a/ortools/math_opt/samples/facility_lp_benders.cc +++ b/ortools/math_opt/samples/facility_lp_benders.cc @@ -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 #include @@ -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::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; @@ -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::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 z; + std::vector 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 x; - for (const auto& edge : network.edges()) { - const Variable x_edge = model.AddContinuousVariable(0.0, kInf); + absl::flat_hash_map 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 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 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 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 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> New( + FacilityLocationInstance instance, math_opt::SolverType solver_type); + + absl::StatusOr> Solve( + const std::vector& z_values, double w_value, + double fist_stage_objective); + + private: + SecondStageSolver(FacilityLocationInstance instance, + math_opt::SolveParameters second_stage_params); + + absl::StatusOr OptimalityCut( + const math_opt::SolveResult& second_stage_result); + absl::StatusOr 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 x_; + std::vector supply_constraints_; + std::vector demand_constraints_; + std::unique_ptr solver_; +}; + +absl::StatusOr 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 demand_constraints; + return parameters; +} +absl::StatusOr> 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 second_stage_solver = + absl::WrapUnique( + new SecondStageSolver(std::move(instance), parameters)); + ASSIGN_OR_RETURN(std::unique_ptr 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 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> SecondStageSolver::Solve( + const std::vector& 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 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 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 first_stage_solver, + math_opt::IncrementalSolver::New(&first_stage.model, solver_type)); + // Setup second stage solver. + ASSIGN_OR_RETURN(std::unique_ptr second_stage_solver, + SecondStageSolver::New(instance, solver_type)); // Start Benders int iteration = 0; double best_upper_bound = kInf; - const std::unique_ptr first_stage_solver = - IncrementalSolver::New(first_stage_model, SolverType::kGurobi).value(); - const std::unique_ptr second_stage_solver = - IncrementalSolver::New(second_stage_model, SolverType::kGurobi).value(); + std::vector 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; } diff --git a/ortools/math_opt/samples/integer_programming.cc b/ortools/math_opt/samples/integer_programming.cc index e5b4cba694..2ced30b269 100644 --- a/ortools/math_opt/samples/integer_programming.cc +++ b/ortools/math_opt/samples/integer_programming.cc @@ -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::infinity(); @@ -39,12 +37,12 @@ constexpr double kInf = std::numeric_limits::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; } diff --git a/ortools/math_opt/samples/lagrangian_relaxation.cc b/ortools/math_opt/samples/lagrangian_relaxation.cc index 41291dbcdf..bb42e6b8b1 100644 --- a/ortools/math_opt/samples/lagrangian_relaxation.cc +++ b/ortools/math_opt/samples/lagrangian_relaxation.cc @@ -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("LagrangianProblem")) {} - std::unique_ptr model; - LinearExpression cost; - LinearExpression resource_1; - LinearExpression resource_2; - std::vector flow_vars; + FlowModel() : model(std::make_unique("LagrangianProblem")) {} + std::unique_ptr model; + math_opt::LinearExpression cost; + math_opt::LinearExpression resource_1; + math_opt::LinearExpression resource_2; + std::vector 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 out_flow(graph.num_nodes); - std::vector in_flow(graph.num_nodes); + std::vector out_flow(graph.num_nodes); + std::vector 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 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 out_flow(graph.num_nodes); - std::vector in_flow(graph.num_nodes); + std::vector out_flow(graph.num_nodes); + std::vector 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& 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& 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 mu; - std::vector grad_mu; + std::vector 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& 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& 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; } diff --git a/ortools/math_opt/samples/linear_programming.cc b/ortools/math_opt/samples/linear_programming.cc index 8ef6a691a6..3104db9c77 100644 --- a/ortools/math_opt/samples/linear_programming.cc +++ b/ortools/math_opt/samples/linear_programming.cc @@ -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::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::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 x; + std::vector x; for (int j = 0; j < 3; j++) { x.push_back(model.AddContinuousVariable(0.0, kInf, absl::StrCat("x", j))); } // Constraints - std::vector constraints; + std::vector 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; } diff --git a/ortools/math_opt/solvers/BUILD.bazel b/ortools/math_opt/solvers/BUILD.bazel index a64559523e..05f9fe1b67 100644 --- a/ortools/math_opt/solvers/BUILD.bazel +++ b/ortools/math_opt/solvers/BUILD.bazel @@ -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"], diff --git a/ortools/math_opt/solvers/cp_sat_solver.cc b/ortools/math_opt/solvers/cp_sat_solver.cc index 7af8a8c79b..2ca20e6491 100644 --- a/ortools/math_opt/solvers/cp_sat_solver.cc +++ b/ortools/math_opt/solvers/cp_sat_solver.cc @@ -100,6 +100,9 @@ std::vector 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 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> CpSatSolver::New( @@ -323,14 +327,18 @@ absl::StatusOr> 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 CpSatSolver::Solve( @@ -435,6 +443,10 @@ absl::StatusOr 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 variable_ids) + std::vector variable_ids, + std::vector 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 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 diff --git a/ortools/math_opt/solvers/cp_sat_solver.h b/ortools/math_opt/solvers/cp_sat_solver.h index 78aeb1d089..437b02b074 100644 --- a/ortools/math_opt/solvers/cp_sat_solver.h +++ b/ortools/math_opt/solvers/cp_sat_solver.h @@ -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 variable_ids); + CpSatSolver(MPModelProto cp_sat_model, std::vector variable_ids, + std::vector linear_constraint_ids); // Extract the solution from CP-SAT's response. SparseDoubleVectorProto ExtractSolution( absl::Span 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 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 linear_constraint_ids_; }; } // namespace math_opt diff --git a/ortools/math_opt/solvers/glop_solver.cc b/ortools/math_opt/solvers/glop_solver.cc index 69764a3f45..cbe7fb1e8d 100644 --- a/ortools/math_opt/solvers/glop_solver.cc +++ b/ortools/math_opt/solvers/glop_solver.cc @@ -67,6 +67,8 @@ namespace math_opt { namespace { +constexpr double kInf = std::numeric_limits::infinity(); + absl::string_view SafeName(const VariablesProto& variables, int index) { if (variables.names().empty()) { return {}; @@ -299,6 +301,9 @@ absl::StatusOr 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 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 +glop::StrictITIVector IndexToId( + const absl::flat_hash_map& id_map) { + // Guard value used to identify not-yet-set elements of index_to_id. + constexpr int64_t kEmptyId = -1; + glop::StrictITIVector 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 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 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 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 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::infinity(); // Set default status and bounds. solve_stats.mutable_problem_status()->set_primal_status( @@ -760,6 +830,11 @@ absl::StatusOr 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; diff --git a/ortools/math_opt/solvers/glop_solver.h b/ortools/math_opt/solvers/glop_solver.h index 978681fae8..cfccc9ad66 100644 --- a/ortools/math_opt/solvers/glop_solver.h +++ b/ortools/math_opt/solvers/glop_solver.h @@ -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); diff --git a/ortools/math_opt/solvers/glpk/BUILD.bazel b/ortools/math_opt/solvers/glpk/BUILD.bazel new file mode 100644 index 0000000000..0c5c0173e2 --- /dev/null +++ b/ortools/math_opt/solvers/glpk/BUILD.bazel @@ -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", + ], +) diff --git a/ortools/math_opt/solvers/glpk/rays.cc b/ortools/math_opt/solvers/glpk/rays.cc new file mode 100644 index 0000000000..d406c00938 --- /dev/null +++ b/ortools/math_opt/solvers/glpk/rays.cc @@ -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 +#include +#include +#include + +#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 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 inds(num_cstrs + 1); + std::vector 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 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 inds(num_structural_vars + 1); + std::vector 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> 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 diff --git a/ortools/math_opt/solvers/glpk/rays.h b/ortools/math_opt/solvers/glpk/rays.h new file mode 100644 index 0000000000..911305bc3a --- /dev/null +++ b/ortools/math_opt/solvers/glpk/rays.h @@ -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 +#include +#include + +#include "absl/status/statusor.h" + +extern "C" { +#include +} + +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>; + + 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> GlpkComputeUnboundRay(glp_prob* problem); + +} // namespace operations_research::math_opt + +#endif // OR_TOOLS_MATH_OPT_SOLVERS_GLPK_RAYS_H_ diff --git a/ortools/math_opt/solvers/glpk_solver.cc b/ortools/math_opt/solvers/glpk_solver.cc new file mode 100644 index 0000000000..3cff3a26dc --- /dev/null +++ b/ortools/math_opt/solvers/glpk_solver.cc @@ -0,0 +1,1748 @@ +// 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_solver.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "absl/base/thread_annotations.h" +#include "absl/container/flat_hash_map.h" +#include "absl/memory/memory.h" +#include "absl/status/status.h" +#include "absl/status/statusor.h" +#include "absl/strings/str_cat.h" +#include "absl/strings/str_join.h" +#include "absl/strings/string_view.h" +#include "absl/synchronization/mutex.h" +#include "absl/time/clock.h" +#include "absl/time/time.h" +#include "ortools/base/cleanup.h" +#include "ortools/base/logging.h" +#include "ortools/base/protoutil.h" +#include "ortools/base/status_macros.h" +#include "ortools/glpk/glpk_env_deleter.h" +#include "ortools/glpk/glpk_formatters.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" +#include "ortools/math_opt/core/sparse_submatrix.h" +#include "ortools/math_opt/core/sparse_vector_view.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" +#include "ortools/math_opt/solution.pb.h" +#include "ortools/math_opt/solvers/glpk/rays.h" +#include "ortools/math_opt/solvers/message_callback_data.h" +#include "ortools/math_opt/sparse_containers.pb.h" +#include "ortools/math_opt/validators/callback_validator.h" +#include "ortools/port/proto_utils.h" + +namespace operations_research { +namespace math_opt { + +namespace { + +constexpr double kInf = std::numeric_limits::infinity(); + +// Bounds of rows or columns. +struct Bounds { + double lower = -kInf; + double upper = kInf; +}; + +// Sets either a row or a column bounds. The index k is the one-based index of +// the row or the column. +// +// The Dimension type should be either GlpkSolver::Variable or +// GlpkSolver::LinearConstraints. +// +// When Dimension::IsInteger() returns true, the bounds are rounded before being +// applied which is mandatory for integer variables (solvers fail if a model +// contains non-integer bounds for integer variables). Thus the integrality of +// variables must be set/updated before calling this function. +template +void SetBounds(glp_prob* const problem, const int k, const Bounds& bounds) { + // GLPK wants integer bounds for integer variables. + const bool is_integer = Dimension::IsInteger(problem, k); + const double lb = is_integer ? std::ceil(bounds.lower) : bounds.lower; + const double ub = is_integer ? std::floor(bounds.upper) : bounds.upper; + int type = GLP_FR; + if (std::isinf(lb) && std::isinf(ub)) { + type = GLP_FR; + } else if (std::isinf(lb)) { + type = GLP_UP; + } else if (std::isinf(ub)) { + type = GLP_LO; + } else if (lb == ub) { + type = GLP_FX; + } else { // Bounds not inf and not equal. + type = GLP_DB; + } + Dimension::kSetBounds(problem, k, type, lb, ub); +} + +// Gets either a row or a column bounds. The index k is the one-based index of +// the row or the column. +// +// The Dimension type should be either GlpkSolver::Variable or +// GlpkSolver::LinearConstraints. +template +Bounds GetBounds(glp_prob* const problem, const int k) { + const int type = Dimension::kGetType(problem, k); + switch (type) { + case GLP_FR: + return {}; + case GLP_LO: + return {.lower = Dimension::kGetLb(problem, k)}; + case GLP_UP: + return {.upper = Dimension::kGetUb(problem, k)}; + case GLP_DB: + case GLP_FX: + return {.lower = Dimension::kGetLb(problem, k), + .upper = Dimension::kGetUb(problem, k)}; + default: + LOG(FATAL) << type; + } +} + +// Updates the bounds of either rows or columns. +// +// The Dimension type should be either GlpkSolver::Variable or +// GlpkSolver::LinearConstraints. +// +// When Dimension::IsInteger() returns true, the bounds are rounded before being +// applied which is mandatory for integer variables (solvers fail if a model +// contains non-integer bounds for integer variables). Thus the integrality of +// variables must be updated before calling this function. +template +void UpdateBounds(glp_prob* const problem, const Dimension& dimension, + const SparseDoubleVectorProto& lower_bounds_proto, + const SparseDoubleVectorProto& upper_bounds_proto) { + const auto lower_bounds = MakeView(lower_bounds_proto); + const auto upper_bounds = MakeView(upper_bounds_proto); + + auto current_lower_bound = lower_bounds.begin(); + auto current_upper_bound = upper_bounds.begin(); + for (;;) { + // Get the smallest unvisited id from either sparse container. + std::optional next_id; + if (current_lower_bound != lower_bounds.end()) { + if (!next_id.has_value() || current_lower_bound->first < *next_id) { + next_id = current_lower_bound->first; + } + } + if (current_upper_bound != upper_bounds.end()) { + if (!next_id.has_value() || current_upper_bound->first < *next_id) { + next_id = current_upper_bound->first; + } + } + + if (!next_id.has_value()) { + // We exhausted all collections. + break; + } + + // Find the corresponding row or column. + const int row_or_col_index = dimension.id_to_index.at(*next_id); + CHECK_EQ(dimension.ids[row_or_col_index - 1], *next_id); + + // Get the updated values for bounds and move the iterator for consumed + // updates. + Bounds bounds = GetBounds(problem, + /*k=*/row_or_col_index); + if (current_lower_bound != lower_bounds.end() && + current_lower_bound->first == *next_id) { + bounds.lower = current_lower_bound->second; + ++current_lower_bound; + } + if (current_upper_bound != upper_bounds.end() && + current_upper_bound->first == *next_id) { + bounds.upper = current_upper_bound->second; + ++current_upper_bound; + } + SetBounds(problem, /*k=*/row_or_col_index, + /*bounds=*/bounds); + } + + CHECK(current_lower_bound == lower_bounds.end()); + CHECK(current_upper_bound == upper_bounds.end()); +} + +// Deletes in-place the data corresponding to the indices of rows/cols. +// +// The vector of one-based indices sorted_deleted_rows_or_cols is expected to be +// sorted and its first element of index 0 is ignored (this is the GLPK +// convention). +template +void DeleteRowOrColData(std::vector& data, + const std::vector& sorted_deleted_rows_or_cols) { + if (sorted_deleted_rows_or_cols.empty()) { + // Avoid looping when not necessary. + return; + } + + std::size_t next_insertion_point = 0; + std::size_t current_row_or_col = 0; + for (std::size_t i = 1; i < sorted_deleted_rows_or_cols.size(); ++i) { + const int deleted_row_or_col = sorted_deleted_rows_or_cols[i]; + for (; current_row_or_col + 1 < deleted_row_or_col; + ++current_row_or_col, ++next_insertion_point) { + DCHECK_LT(current_row_or_col, data.size()); + data[next_insertion_point] = data[current_row_or_col]; + } + // Skip the deleted row/col. + ++current_row_or_col; + } + for (; current_row_or_col < data.size(); + ++current_row_or_col, ++next_insertion_point) { + data[next_insertion_point] = data[current_row_or_col]; + } + data.resize(next_insertion_point); +} + +// Deletes the row or cols of the GLPK problem and returns their indices. As a +// side effect it updates dimension.ids and dimension.id_to_index. +// +// The Dimension type should be either GlpkSolver::Variable or +// GlpkSolver::LinearConstraints. +// +// The returned vector is sorted and the first element (index 0) must be ignored +// (this is the GLPK convention). It can be used with DeleteRowOrColData(). +template +std::vector DeleteRowsOrCols( + glp_prob* const problem, Dimension& dimension, + const google::protobuf::RepeatedField& deleted_ids) { + if (deleted_ids.empty()) { + // This is not only an optimization. Functions glp_del_rows() and + // glp_del_cols() fails if the number of deletion is 0. + return {}; + } + + // Delete GLPK rows or columns. + std::vector deleted_rows_or_cols; + // Functions glp_del_rows() and glp_del_cols() only use values in ranges + // [1,n]. The first element is not used. + deleted_rows_or_cols.reserve(deleted_ids.size() + 1); + deleted_rows_or_cols.push_back(-1); + for (const int64_t deleted_id : deleted_ids) { + deleted_rows_or_cols.push_back(dimension.id_to_index.at(deleted_id)); + } + Dimension::kDelElts(problem, deleted_rows_or_cols.size() - 1, + deleted_rows_or_cols.data()); + + // Since deleted_ids are in strictly increasing order and we allocate + // rows/cols in orders of MathOpt ids; deleted_rows_or_cols should also be + // sorted. + CHECK( + std::is_sorted(deleted_rows_or_cols.begin(), deleted_rows_or_cols.end())); + + // Update the ids vector. + DeleteRowOrColData(dimension.ids, deleted_rows_or_cols); + + // Update the id_to_index map. + for (const int64_t deleted_id : deleted_ids) { + CHECK(dimension.id_to_index.erase(deleted_id)); + } + for (int i = 0; i < dimension.ids.size(); ++i) { + dimension.id_to_index.at(dimension.ids[i]) = i + 1; + } + + return deleted_rows_or_cols; +} + +// Translates the input MathOpt indices in row/column GLPK indices to use with +// glp_load_matrix(). The returned vector first element is always 0 and unused +// as it is required by GLPK (which uses one-based indices for arrays as well). +// +// The id_to_index is supposed to contain GLPK's one-based indices for rows and +// columns. +std::vector MatrixIds( + const google::protobuf::RepeatedField& proto_ids, + const absl::flat_hash_map& id_to_index) { + std::vector ids; + ids.reserve(proto_ids.size() + 1); + // First item (index 0) is not used by GLPK. + ids.push_back(0); + for (const int64_t proto_id : proto_ids) { + ids.push_back(id_to_index.at(proto_id)); + } + return ids; +} + +// Returns a vector of coefficients starting at index 1 (as used by GLPK) to use +// with glp_load_matrix(). The returned vector first element is always 0 and it +// is ignored by GLPK. +std::vector MatrixCoefficients( + const google::protobuf::RepeatedField& proto_coeffs) { + std::vector coeffs(proto_coeffs.size() + 1); + // First item (index 0) is not used by GLPK. + coeffs[0] = 0.0; + std::copy(proto_coeffs.begin(), proto_coeffs.end(), coeffs.begin() + 1); + return coeffs; +} + +// Returns true if the input GLPK problem contains integer variables. +bool IsMip(glp_prob* const problem) { + const int num_vars = glp_get_num_cols(problem); + for (int v = 1; v <= num_vars; ++v) { + if (glp_get_col_kind(problem, v) != GLP_CV) { + return true; + } + } + return false; +} + +// Returns true if the input GLPK problem has no rows and no cols. +bool IsEmpty(glp_prob* const problem) { + return glp_get_num_cols(problem) == 0 && glp_get_num_rows(problem) == 0; +} + +// Returns a sparse vector with the values returned by the getter for the input +// ids and taking into account the provided filter. +SparseDoubleVectorProto FilteredVector(glp_prob* const problem, + const SparseVectorFilterProto& filter, + const std::vector& ids, + double (*const getter)(glp_prob*, int)) { + SparseDoubleVectorProto vec; + vec.mutable_ids()->Reserve(ids.size()); + vec.mutable_values()->Reserve(ids.size()); + + SparseVectorFilterPredicate predicate(filter); + for (int i = 0; i < ids.size(); ++i) { + const double value = getter(problem, i + 1); + if (predicate.AcceptsAndUpdate(ids[i], value)) { + vec.add_ids(ids[i]); + vec.add_values(value); + } + } + return vec; +} + +// Returns the ray data the corresponds to element id having the given value and +// all other elements of ids having 0. +SparseDoubleVectorProto FilteredRay(const SparseVectorFilterProto& filter, + const std::vector& ids, + const std::vector& values) { + CHECK_EQ(ids.size(), values.size()); + SparseDoubleVectorProto vec; + SparseVectorFilterPredicate predicate(filter); + for (int i = 0; i < ids.size(); ++i) { + if (predicate.AcceptsAndUpdate(ids[i], values[i])) { + vec.add_ids(ids[i]); + vec.add_values(values[i]); + } + } + return vec; +} + +// Sets the parameters shared between MIP and LP and returns warnings for bad +// parameters. +// +// The input Parameters type should be glp_smcp (for LP), glp_iptcp (for LP with +// interior point) or glp_iocp (for MIP). +template +absl::Status SetSharedParameters(const SolveParametersProto& parameters, + const bool has_message_callback, + Parameters& glpk_parameters) { + std::vector warnings; + if (parameters.has_threads() && parameters.threads() > 1) { + warnings.push_back( + absl::StrCat("GLPK only supports parameters.threads = 1; value ", + parameters.threads(), " is not supported")); + } + if (parameters.enable_output() || has_message_callback) { + glpk_parameters.msg_lev = GLP_MSG_ALL; + } else { + glpk_parameters.msg_lev = GLP_MSG_OFF; + } + if (parameters.has_node_limit()) { + warnings.push_back("Parameter node_limit not supported by GLPK"); + } + if (parameters.has_objective_limit()) { + warnings.push_back("Parameter objective_limit not supported by GLPK"); + } + if (parameters.has_best_bound_limit()) { + warnings.push_back("Parameter best_bound_limit not supported by GLPK"); + } + if (parameters.has_cutoff_limit()) { + warnings.push_back("Parameter cutoff_limit not supported by GLPK"); + } + if (parameters.has_solution_limit()) { + warnings.push_back("Parameter solution_limit not supported by GLPK"); + } + if (!warnings.empty()) { + return absl::InvalidArgumentError(absl::StrJoin(warnings, "; ")); + } + return absl::OkStatus(); +} + +// Sets the time limit parameter which is only supported by some LP algorithm +// and MIP, but not by interior point. +// +// The input Parameters type should be glp_smcp (for LP), or glp_iocp (for MIP). +template +void SetTimeLimitParameter(const SolveParametersProto& parameters, + Parameters& glpk_parameters) { + if (parameters.has_time_limit()) { + const int64_t time_limit_ms = absl::ToInt64Milliseconds( + util_time::DecodeGoogleApiProto(parameters.time_limit()).value()); + glpk_parameters.tm_lim = static_cast(std::min( + static_cast(std::numeric_limits::max()), time_limit_ms)); + } +} + +// Sets the LP specific parameters and returns an InvalidArgumentError for +// invalid parameters or parameter values. +absl::Status SetLPParameters(const SolveParametersProto& parameters, + glp_smcp& glpk_parameters) { + std::vector warnings; + switch (parameters.presolve()) { + case EMPHASIS_UNSPECIFIED: + // Keep the default. + // + // TODO(b/187027049): default is off, which may be surprising for users. + break; + case EMPHASIS_OFF: + glpk_parameters.presolve = GLP_OFF; + break; + default: + glpk_parameters.presolve = GLP_ON; + break; + } + switch (parameters.lp_algorithm()) { + case LP_ALGORITHM_UNSPECIFIED: + break; + case LP_ALGORITHM_PRIMAL_SIMPLEX: + glpk_parameters.meth = GLP_PRIMAL; + break; + case LP_ALGORITHM_DUAL_SIMPLEX: + // Use GLP_DUALP to switch back to primal simplex if the dual simplex + // fails. + // + // TODO(b/187027049): GLPK also supports GLP_DUAL to only try dual + // simplex. We should have an option to support it. + glpk_parameters.meth = GLP_DUALP; + break; + default: + warnings.push_back(absl::StrCat( + "GLPK does not support ", + operations_research::ProtoEnumToString(parameters.lp_algorithm()), + " for parameters.lp_algorithm")); + break; + } + if (!warnings.empty()) { + return absl::InvalidArgumentError(absl::StrJoin(warnings, "; ")); + } + return absl::OkStatus(); +} + +class MipCallbackData { + public: + explicit MipCallbackData(SolveInterrupter* const interrupter) + : interrupter_(interrupter) {} + + void Callback(glp_tree* const tree) { + // We only update the best bound on some specific events since it makes a + // traversal of all active nodes. + switch (glp_ios_reason(tree)) { + case GLP_ISELECT: + // The ISELECT call is the first one that happens after a node has been + // split on two sub-nodes (IBRANCH) with updated `bound`s based on the + // integer value of the branched variable. + case GLP_IBINGO: + // We found a new integer solution, the `bound` has been updated. + case GLP_IROWGEN: + // The IROWGEN call is the first one that happens on a current node + // after the relaxed problem has been solved and the `bound` field + // updated. + // + // Note that the model/cut pool changes done in IROWGEN and ICUTGEN have + // no influence on the `bound` and IROWGEN is the first call to happen. + if (const int best_node = glp_ios_best_node(tree); best_node != 0) { + best_bound_ = glp_ios_node_bound(tree, best_node); + } + break; + default: + // We can ignore: + // - IPREPRO: since the `bound` of the current node has not been + // computed yet. + // - IHEUR: since we have IBINGO if the integer solution is better. + // - ICUTGEN: since the `bound` is not updated with the rows added at + // IROWGEN so we would get the same best bound. + // - IBRANCH: since the sub-nodes will be created after that and their + // `bound`s taken into account at ISELECT. + break; + } + if (interrupter_ != nullptr && interrupter_->IsInterrupted()) { + glp_ios_terminate(tree); + interrupted_by_interrupter_ = true; + return; + } + } + + bool HasBeenInterruptedByInterrupter() const { + return interrupted_by_interrupter_.load(); + } + + std::optional best_bound() const { return best_bound_; } + + private: + // Optional interrupter. + SolveInterrupter* const interrupter_; + + // Set to true if glp_ios_terminate() has been called due to the interrupter. + std::atomic interrupted_by_interrupter_ = false; + + // Set on each callback that may update the best bound. + std::optional best_bound_; +}; + +void MipCallback(glp_tree* const tree, void* const info) { + static_cast(info)->Callback(tree); +} + +// Returns the MathOpt ids of the rows/columns with lower_bound > upper_bound. +InvertedBounds ListInvertedBounds( + glp_prob* const problem, const std::vector& variable_ids, + const std::vector& linear_constraint_ids) { + InvertedBounds inverted_bounds; + + const int num_cols = glp_get_num_cols(problem); + for (int c = 1; c <= num_cols; ++c) { + if (glp_get_col_lb(problem, c) > glp_get_col_ub(problem, c)) { + inverted_bounds.variables.push_back(variable_ids[c - 1]); + } + } + + const int num_rows = glp_get_num_rows(problem); + for (int r = 1; r <= num_rows; ++r) { + if (glp_get_row_lb(problem, r) > glp_get_row_ub(problem, r)) { + inverted_bounds.linear_constraints.push_back( + linear_constraint_ids[r - 1]); + } + } + + return inverted_bounds; +} + +// Returns the termination reason based on the current MIP data of the problem +// assuming that the last call to glp_intopt() returned 0 and that the model has +// not been modified since. +absl::StatusOr MipTerminationOnSuccess( + glp_prob* const problem) { + const int status = glp_mip_status(problem); + switch (status) { + case GLP_OPT: + return TerminateForReason(TERMINATION_REASON_OPTIMAL); + case GLP_FEAS: + return FeasibleTermination(LIMIT_UNDETERMINED, + "glp_mip_status() returned GLP_FEAS"); + case GLP_NOFEAS: + return TerminateForReason(TERMINATION_REASON_INFEASIBLE); + default: + return absl::InternalError( + absl::StrCat("glp_intopt() returned 0 but glp_mip_status()" + "returned the unexpected value ", + SolutionStatusString(status))); + } +} + +// Returns the termination reason based on the current interior point data of +// the problem assuming that the last call to glp_interior() returned 0 and that +// the model has not been modified since. +absl::StatusOr InteriorTerminationOnSuccess( + glp_prob* const problem) { + const int status = glp_ipt_status(problem); + switch (status) { + case GLP_OPT: + return TerminateForReason(TERMINATION_REASON_OPTIMAL); + case GLP_INFEAS: + return NoSolutionFoundTermination(LIMIT_UNDETERMINED, + "glp_ipt_status() returned GLP_INFEAS"); + case GLP_NOFEAS: + // Documentation in glpapi08.c for glp_ipt_status says this status means + // "no feasible solution exists", but the Reference Manual for GLPK + // Version 5.0 clarifies that it means "no feasible primal-dual solution + // exists." (See also the comment in glpipm.c when ipm_solve returns 1). + // Hence, GLP_NOFEAS corresponds to the solver claiming that either the + // primal problem, the dual problem (or both) are infeasible. Under this + // condition if the primal is feasible, then the dual must be infeasible + // and hence the primal is unbounded. + return TerminateForReason(TERMINATION_REASON_INFEASIBLE_OR_UNBOUNDED); + default: + return absl::InternalError( + absl::StrCat("glp_interior() returned 0 but glp_ipt_status()" + "returned the unexpected value ", + SolutionStatusString(status))); + } +} + +// Returns the termination reason based on the current interior point data of +// the problem assuming that the last call to glp_simplex() returned 0 and that +// the model has not been modified since. +absl::StatusOr SimplexTerminationOnSuccess( + glp_prob* const problem) { + // Here we don't use glp_get_status() since it is biased towards the primal + // simplex algorithm. For example if the dual simplex returns GLP_NOFEAS for + // the dual and GLP_INFEAS for the primal then glp_get_status() returns + // GLP_INFEAS. This is misleading since the dual successfully determined that + // the problem was dual infeasible. So here we use the two statuses of the + // primal and the dual to get a better result (the glp_get_status() only + // combines them anyway, it does not have any other benefit). + const int prim_status = glp_get_prim_stat(problem); + const int dual_status = glp_get_dual_stat(problem); + + // Returns the undetermined limit for cases where we can't draw a conclusion. + const auto undetermined_limit = [&]() { + const std::string detail = absl::StrCat( + "glp_get_prim_stat() returned ", SolutionStatusString(prim_status), + " and glp_get_dual_stat() returned ", + SolutionStatusString(dual_status)); + if (prim_status == GLP_FEAS) { + return FeasibleTermination(LIMIT_UNDETERMINED, detail); + } + return NoSolutionFoundTermination(LIMIT_UNDETERMINED, detail); + }; + + // Returns a status error indicating that glp_get_dual_stat() returned an + // unexpected value. + const auto unexpected_dual_stat = [&]() { + return absl::InternalError( + absl::StrCat("glp_simplex() returned 0 but glp_get_dual_stat() " + "returned the unexpected value ", + SolutionStatusString(dual_status))); + }; + + switch (prim_status) { + case GLP_FEAS: + switch (dual_status) { + case GLP_FEAS: + // Dual feasibility here means that the solution is dual feasible + // (correct signs of the residual costs) and that the complementary + // slackness condition are respected. Hence the solution is optimal. + return TerminateForReason(TERMINATION_REASON_OPTIMAL); + case GLP_INFEAS: + return undetermined_limit(); + case GLP_NOFEAS: + return TerminateForReason(TERMINATION_REASON_UNBOUNDED); + default: + return unexpected_dual_stat(); + } + case GLP_INFEAS: + switch (dual_status) { + case GLP_FEAS: + case GLP_INFEAS: + return undetermined_limit(); + case GLP_NOFEAS: + return TerminateForReason(TERMINATION_REASON_INFEASIBLE_OR_UNBOUNDED); + default: + return unexpected_dual_stat(); + } + case GLP_NOFEAS: + switch (dual_status) { + case GLP_FEAS: + case GLP_INFEAS: + case GLP_NOFEAS: + // Dual being feasible (GLP_FEAS) here would lead to dual unbounded; + // but this does not exist as a reason. + // + // If both the primal and dual are proven infeasible (GLP_NOFEAS), the + // primal wins. Maybe GLPK does never return that though since it + // implements either primal or dual simplex algorithm but does not + // combine both of them. + return TerminateForReason(TERMINATION_REASON_INFEASIBLE); + default: + return unexpected_dual_stat(); + } + default: + return absl::InternalError( + absl::StrCat("glp_simplex() returned 0 but glp_get_prim_stat() " + "returned the unexpected value ", + SolutionStatusString(prim_status))); + } +} + +// Returns the termination reason based on the return code rc of calling fn_name +// function which is glp_simplex(), glp_interior() or glp_intopt(). +// +// For return code 0 which means successful solve, the function +// termination_on_success is called to build the termination. Other return +// values (errors) are dealt with specifically. +// +// For glp_intopt(), the optional mip_cb_data is used to test for interruption +// and the LIMIT_INTERRUPTED is set if the interrupter has been triggered (even +// if the return code is 0). +// +// The parameters `(variable|linear_constraint)_ids` are the +// `GlpkSolver::(LinearConstraints|Variables)::ids`. +absl::StatusOr BuildTermination( + glp_prob* const problem, const std::string_view fn_name, const int rc, + const std::function(glp_prob*)> + termination_on_success, + MipCallbackData* const mip_cb_data, const bool has_feasible_solution, + const std::vector& variable_ids, + const std::vector& linear_constraint_ids) { + if (mip_cb_data != nullptr && + mip_cb_data->HasBeenInterruptedByInterrupter()) { + return TerminateForLimit(LIMIT_INTERRUPTED, + /*feasible=*/has_feasible_solution); + } + + // TODO(b/187027049): see if GLP_EOBJLL and GLP_EOBJUL should be handled with + // dual simplex. + switch (rc) { + case 0: + return termination_on_success(problem); + case GLP_EBOUND: { + // GLP_EBOUND is returned when a variable or a constraint has the GLP_DB + // bounds type and lower_bound >= upper_bound. The code in this file makes + // sure we don't use GLP_DB but GLP_FX when lower_bound == upper_bound + // thus we expect GLP_EBOUND only when lower_bound > upper_bound. + RETURN_IF_ERROR( + ListInvertedBounds(problem, + /*variable_ids=*/variable_ids, + /*linear_constraint_ids=*/linear_constraint_ids) + .ToStatus()); + return util::InternalErrorBuilder() + << fn_name << "() returned `" << ReturnCodeString(rc) + << "` but the model does not contain variables with inverted " + "bounds"; + } + case GLP_EITLIM: + return TerminateForLimit(LIMIT_ITERATION, + /*feasible=*/has_feasible_solution); + case GLP_ETMLIM: + return TerminateForLimit(LIMIT_TIME, /*feasible=*/has_feasible_solution); + case GLP_EMIPGAP: + return TerminateForReason( + TERMINATION_REASON_OPTIMAL, + // absl::StrCat() does not compile with std::string_view on WASM. + // + absl::StrCat(std::string(fn_name), "() returned ", + ReturnCodeString(rc))); + case GLP_ESTOP: + return TerminateForLimit(LIMIT_INTERRUPTED, + /*feasible=*/has_feasible_solution); + case GLP_ENOPFS: + // With presolve on, this error is returned if the LP has no feasible + // solution. + return TerminateForReason(TERMINATION_REASON_INFEASIBLE); + case GLP_ENODFS: + // With presolve on, this error is returned if the LP has no dual + // feasible solution. + return TerminateForReason(TERMINATION_REASON_INFEASIBLE_OR_UNBOUNDED); + case GLP_ENOCVG: + // Very slow convergence/divergence (for glp_interior). + return TerminateForLimit(LIMIT_SLOW_PROGRESS, + /*feasible=*/has_feasible_solution); + case GLP_EINSTAB: + // Numeric stability solving Newtonian system (for glp_interior). + return TerminateForReason( + TERMINATION_REASON_NUMERICAL_ERROR, + // absl::StrCat() does not compile with std::string_view on WASM. + // + absl::StrCat(std::string(fn_name), "() returned ", + ReturnCodeString(rc), + " which means that there is a numeric stability issue " + "solving Newtonian system")); + default: + return util::InternalErrorBuilder() + << fn_name + << "() returned unexpected value: " << ReturnCodeString(rc); + } +} + +class TermHookData { + public: + explicit TermHookData(SolverInterface::MessageCallback callback) + : callback_(std::move(callback)) {} + + void Parse(const std::string_view message) { + // Here we keep the lock while calling the callback. This should not be an + // issue since we don't expect code in a message callback to trigger a new + // message. On top of that, for proper interleaving it may be better to use + // the lock anyway. + const absl::MutexLock lock(&mutex_); + std::vector new_lines = buffer_.Parse(message); + if (!new_lines.empty()) { + callback_(new_lines); + } + } + + // Flushes the buffer and calls the callback if the result is not empty. + void Flush() { + // See comment in Parse() about holding the lock while calling the callback. + const absl::MutexLock lock(&mutex_); + std::vector new_lines = buffer_.Flush(); + if (!new_lines.empty()) { + callback_(new_lines); + } + } + + private: + absl::Mutex mutex_; + MessageCallbackData buffer_ ABSL_GUARDED_BY(mutex_); + const SolverInterface::MessageCallback callback_; +}; + +// Callback for glp_term_hook(). +// +// It expects `info` to be a pointer on a TermHookData. +int TermHook(void* const info, const char* const message) { + static_cast(info)->Parse(message); + + // Returns non-zero to remove any terminal output. + return 1; +} + +// Returns the objective offset. This is used as a placeholder for function +// returning the objective value for solve method not supporting solving empty +// models (glp_exact() and glp_interior()). +double OffsetOnlyObjVal(glp_prob* const problem) { + return glp_get_obj_coef(problem, 0); +} + +// Returns GLP_OPT. This is used as a placeholder for function returning the +// status for solve method not supporting solving empty models (glp_exact() and +// glp_interior()). +int OptStatus(glp_prob*) { return GLP_OPT; } + +// Returns the error when a model or an update contains a quadratic objective. +absl::Status QuadraticObjectiveError() { + return absl::InvalidArgumentError( + "GLPK does not support quadratic objectives"); +} + +} // namespace + +absl::StatusOr> GlpkSolver::New( + const ModelProto& model, const InitArgs& init_args) { + if (!model.objective().quadratic_coefficients().row_ids().empty()) { + return QuadraticObjectiveError(); + } + return absl::WrapUnique(new GlpkSolver(model)); +} + +GlpkSolver::GlpkSolver(const ModelProto& model) : problem_(glp_create_prob()) { + // Make sure glp_free_env() is called at the exit of the current thread. + SetupGlpkEnvAutomaticDeletion(); + + glp_set_prob_name(problem_, TruncateAndQuoteGLPKName(model.name()).c_str()); + + AddVariables(model.variables()); + + AddLinearConstraints(model.linear_constraints()); + + glp_set_obj_dir(problem_, model.objective().maximize() ? GLP_MAX : GLP_MIN); + // Glpk uses index 0 for the "shift" of the objective. + glp_set_obj_coef(problem_, 0, model.objective().offset()); + for (const auto [v, coeff] : + MakeView(model.objective().linear_coefficients())) { + const int col_index = variables_.id_to_index.at(v); + CHECK_EQ(variables_.ids[col_index - 1], v); + glp_set_obj_coef(problem_, col_index, coeff); + } + + const SparseDoubleMatrixProto& proto_matrix = + model.linear_constraint_matrix(); + glp_load_matrix( + problem_, proto_matrix.row_ids_size(), + MatrixIds(proto_matrix.row_ids(), linear_constraints_.id_to_index).data(), + MatrixIds(proto_matrix.column_ids(), variables_.id_to_index).data(), + MatrixCoefficients(proto_matrix.coefficients()).data()); +} + +GlpkSolver::~GlpkSolver() { glp_delete_prob(problem_); } + +namespace { + +ProblemStatusProto GetMipProblemStatusProto(const int rc, const int mip_status, + const bool has_finite_dual_bound) { + ProblemStatusProto problem_status; + problem_status.set_primal_status(FEASIBILITY_STATUS_UNDETERMINED); + problem_status.set_dual_status(FEASIBILITY_STATUS_UNDETERMINED); + + switch (rc) { + case GLP_ENOPFS: + problem_status.set_primal_status(FEASIBILITY_STATUS_INFEASIBLE); + return problem_status; + case GLP_ENODFS: + problem_status.set_dual_status(FEASIBILITY_STATUS_INFEASIBLE); + return problem_status; + } + + switch (mip_status) { + case GLP_OPT: + problem_status.set_primal_status(FEASIBILITY_STATUS_FEASIBLE); + problem_status.set_dual_status(FEASIBILITY_STATUS_FEASIBLE); + return problem_status; + case GLP_FEAS: + problem_status.set_primal_status(FEASIBILITY_STATUS_FEASIBLE); + break; + case GLP_NOFEAS: + problem_status.set_primal_status(FEASIBILITY_STATUS_INFEASIBLE); + break; + } + + if (has_finite_dual_bound) { + problem_status.set_dual_status(FEASIBILITY_STATUS_FEASIBLE); + } + return problem_status; +} + +absl::StatusOr TranslateProblemStatus( + const int glpk_status, const absl::string_view fn_name) { + switch (glpk_status) { + case GLP_FEAS: + return FEASIBILITY_STATUS_FEASIBLE; + case GLP_NOFEAS: + return FEASIBILITY_STATUS_INFEASIBLE; + case GLP_INFEAS: + case GLP_UNDEF: + return FEASIBILITY_STATUS_UNDETERMINED; + default: + return absl::InternalError( + absl::StrCat(fn_name, " returned the unexpected value ", + SolutionStatusString(glpk_status))); + } +} + +// Builds problem status from: +// * glp_simplex_rc: code returned by glp_simplex. +// * glpk_primal_status: primal status returned by glp_get_prim_stat. +// * glpk_dual_status: dual status returned by glp_get_dual_stat. +absl::StatusOr GetSimplexProblemStatusProto( + const int glp_simplex_rc, const int glpk_primal_status, + const int glpk_dual_status) { + ProblemStatusProto problem_status; + problem_status.set_primal_status(FEASIBILITY_STATUS_UNDETERMINED); + problem_status.set_dual_status(FEASIBILITY_STATUS_UNDETERMINED); + + switch (glp_simplex_rc) { + case GLP_ENOPFS: + // LP presolver concluded primal infeasibility. + problem_status.set_primal_status(FEASIBILITY_STATUS_INFEASIBLE); + return problem_status; + case GLP_ENODFS: + // LP presolver concluded dual infeasibility. + problem_status.set_dual_status(FEASIBILITY_STATUS_INFEASIBLE); + return problem_status; + default: { + // Get primal status from basic solution. + ASSIGN_OR_RETURN( + const FeasibilityStatusProto primal_status, + TranslateProblemStatus(glpk_primal_status, "glp_get_prim_stat")); + problem_status.set_primal_status(primal_status); + + // Get dual status from basic solution. + ASSIGN_OR_RETURN( + const FeasibilityStatusProto dual_status, + TranslateProblemStatus(glpk_dual_status, "glp_get_dual_stat")); + problem_status.set_dual_status(dual_status); + return problem_status; + } + } +} + +absl::StatusOr GetBarrierProblemStatusProto( + const int glp_interior_rc, const int ipt_status) { + ProblemStatusProto problem_status; + problem_status.set_primal_status(FEASIBILITY_STATUS_UNDETERMINED); + problem_status.set_dual_status(FEASIBILITY_STATUS_UNDETERMINED); + + switch (glp_interior_rc) { + case 0: + // We only use the glp_ipt_status() result when glp_interior() returned 0. + switch (ipt_status) { + case GLP_OPT: + problem_status.set_primal_status(FEASIBILITY_STATUS_FEASIBLE); + problem_status.set_dual_status(FEASIBILITY_STATUS_FEASIBLE); + return problem_status; + case GLP_INFEAS: + return problem_status; + case GLP_NOFEAS: + problem_status.set_primal_or_dual_infeasible(true); + return problem_status; + case GLP_UNDEF: + return problem_status; + default: + return absl::InternalError( + absl::StrCat("glp_ipt_status returned the unexpected value ", + SolutionStatusString(ipt_status))); + } + default: + return problem_status; + } +} + +} // namespace + +absl::StatusOr GlpkSolver::Solve( + const SolveParametersProto& parameters, + const ModelSolveParametersProto& model_parameters, + MessageCallback message_cb, + const CallbackRegistrationProto& callback_registration, const Callback cb, + SolveInterrupter* const interrupter) { + // Make sure glp_free_env() is called at the exit of the current thread. The + // environment gets created automatically for messages for example. + SetupGlpkEnvAutomaticDeletion(); + + const absl::Time start = absl::Now(); + + RETURN_IF_ERROR(CheckRegisteredCallbackEvents(callback_registration, + /*supported_events=*/{})); + + std::unique_ptr term_hook_data; + if (message_cb != nullptr) { + term_hook_data = std::make_unique(std::move(message_cb)); + + // Note that glp_term_hook() uses get_env_ptr() that relies on thread local + // storage to have a different environment per thread. Thus using + // glp_term_hook() is thread-safe. + // + glp_term_hook(TermHook, term_hook_data.get()); + } + + // We must reset the term hook when before exiting or before flushing the last + // unfinished line. + auto message_cb_cleanup = absl::MakeCleanup([&]() { + if (term_hook_data != nullptr) { + glp_term_hook(/*func=*/nullptr, /*info=*/nullptr); + } + }); + + SolveResultProto result; + + const bool is_mip = IsMip(problem_); + + // We need to use different functions depending on the solve function we used + // (or placeholders if no solve function was called in case of empty models). + int (*get_prim_stat)(glp_prob*) = nullptr; + double (*obj_val)(glp_prob*) = nullptr; + double (*col_val)(glp_prob*, int) = nullptr; + + int (*get_dual_stat)(glp_prob*) = nullptr; + double (*row_dual)(glp_prob*, int) = nullptr; + double (*col_dual)(glp_prob*, int) = nullptr; + + const bool maximize = glp_get_obj_dir(problem_) == GLP_MAX; + double best_dual_bound = maximize ? kInf : -kInf; + + // Here we use different solve algorithms depending on the type of problem: + // * For MIPs: glp_intopt() + // * For LPs: + // * glp_interior() when using BARRIER LP algorithm + // * glp_simplex() for other LP algorithms. + // + // These solve algorithms have dedicated data segments in glp_prob which use + // different access functions to get the solution; hence each branch will set + // the corresponding function pointers accordingly. They also use a custom + // struct for parameters that will be initialized and passed to the algorithm. + if (is_mip) { + get_prim_stat = glp_mip_status; + obj_val = glp_mip_obj_val; + col_val = glp_mip_col_val; + + glp_iocp glpk_parameters; + glp_init_iocp(&glpk_parameters); + RETURN_IF_ERROR(SetSharedParameters( + parameters, + /*has_message_callback=*/term_hook_data != nullptr, glpk_parameters)); + SetTimeLimitParameter(parameters, glpk_parameters); + // TODO(b/187027049): glp_intopt with presolve off requires an optional + // solution of the relaxed problem. Here we simply always enable pre-solve + // but we should support disabling the presolve and call glp_simplex() in + // that case. + glpk_parameters.presolve = GLP_ON; + MipCallbackData mip_cb_data(interrupter); + glpk_parameters.cb_func = MipCallback; + glpk_parameters.cb_info = &mip_cb_data; + const int rc = glp_intopt(problem_, &glpk_parameters); + const int mip_status = glp_mip_status(problem_); + const bool has_feasible_solution = + mip_status == GLP_OPT || mip_status == GLP_FEAS; + ASSIGN_OR_RETURN( + *result.mutable_termination(), + BuildTermination(problem_, "glp_intopt", rc, MipTerminationOnSuccess, + &mip_cb_data, + /*has_feasible_solution=*/has_feasible_solution, + /*variable_ids=*/variables_.ids, + /*linear_constraint_ids=*/linear_constraints_.ids)); + if (mip_cb_data.best_bound().has_value()) { + best_dual_bound = *mip_cb_data.best_bound(); + } + *result.mutable_solve_stats()->mutable_problem_status() = + GetMipProblemStatusProto(rc, mip_status, + std::isfinite(best_dual_bound)); + } else { + if (parameters.lp_algorithm() == LP_ALGORITHM_BARRIER) { + get_prim_stat = glp_ipt_status; + obj_val = glp_ipt_obj_val; + col_val = glp_ipt_col_prim; + + get_dual_stat = glp_ipt_status; + row_dual = glp_ipt_row_dual; + col_dual = glp_ipt_col_dual; + + glp_iptcp glpk_parameters; + glp_init_iptcp(&glpk_parameters); + if (parameters.has_time_limit()) { + return absl::InvalidArgumentError( + "Parameter time_limit not supported by GLPK for interior point " + "algorithm."); + } + RETURN_IF_ERROR(SetSharedParameters( + parameters, + /*has_message_callback=*/term_hook_data != nullptr, glpk_parameters)); + + // glp_interior() does not support being called with an empty model and + // returns GLP_EFAIL. Thus we use placeholders in that case. + if (IsEmpty(problem_)) { + get_prim_stat = OptStatus; + get_dual_stat = OptStatus; + obj_val = OffsetOnlyObjVal; + *result.mutable_termination() = TerminateForReason( + TERMINATION_REASON_OPTIMAL, + "glp_interior() not called since the model is empty"); + result.mutable_solve_stats() + ->mutable_problem_status() + ->set_primal_status(FEASIBILITY_STATUS_FEASIBLE); + result.mutable_solve_stats()->mutable_problem_status()->set_dual_status( + FEASIBILITY_STATUS_FEASIBLE); + } else { + // TODO(b/187027049): add solver specific parameters for + // glp_iptcp.ord_alg. + const int glp_interior_rc = glp_interior(problem_, &glpk_parameters); + const int ipt_status = glp_ipt_status(problem_); + const bool has_feasible_solution = ipt_status == GLP_OPT; + ASSIGN_OR_RETURN( + *result.mutable_termination(), + BuildTermination( + problem_, "glp_interior", glp_interior_rc, + InteriorTerminationOnSuccess, + /*mip_cb_data=*/nullptr, + /*has_feasible_solution=*/has_feasible_solution, + /*variable_ids=*/variables_.ids, + /*linear_constraint_ids=*/linear_constraints_.ids)); + ASSIGN_OR_RETURN( + *result.mutable_solve_stats()->mutable_problem_status(), + GetBarrierProblemStatusProto(/*glp_interior_rc=*/glp_interior_rc, + /*ipt_status=*/ipt_status)); + } + } else { + get_prim_stat = glp_get_prim_stat; + obj_val = glp_get_obj_val; + col_val = glp_get_col_prim; + + get_dual_stat = glp_get_dual_stat; + row_dual = glp_get_row_dual; + col_dual = glp_get_col_dual; + + glp_smcp glpk_parameters; + glp_init_smcp(&glpk_parameters); + RETURN_IF_ERROR(SetSharedParameters( + parameters, + /*has_message_callback=*/term_hook_data != nullptr, glpk_parameters)); + SetTimeLimitParameter(parameters, glpk_parameters); + RETURN_IF_ERROR(SetLPParameters(parameters, glpk_parameters)); + + // TODO(b/187027049): add option to use glp_exact(). + const int glp_simplex_rc = glp_simplex(problem_, &glpk_parameters); + const int prim_stat = glp_get_prim_stat(problem_); + const bool has_feasible_solution = prim_stat == GLP_FEAS; + ASSIGN_OR_RETURN( + *result.mutable_termination(), + BuildTermination(problem_, "glp_simplex", glp_simplex_rc, + SimplexTerminationOnSuccess, + /*mip_cb_data=*/nullptr, + /*has_feasible_solution=*/has_feasible_solution, + /*variable_ids=*/variables_.ids, + /*linear_constraint_ids=*/linear_constraints_.ids)); + + ASSIGN_OR_RETURN(*result.mutable_solve_stats()->mutable_problem_status(), + GetSimplexProblemStatusProto( + /*glp_simplex_rc=*/glp_simplex_rc, + /*glpk_primal_status=*/prim_stat, + /*glpk_dual_status=*/glp_get_dual_stat(problem_))); + VLOG(1) << "glp_get_status: " + << SolutionStatusString(glp_get_status(problem_)) + << " glp_get_prim_stat: " << SolutionStatusString(prim_stat) + << " glp_get_dual_stat: " + << SolutionStatusString(glp_get_dual_stat(problem_)); + } + } + + // Flushes the potential last unfinished line. + if (term_hook_data != nullptr) { + // Make sure no calls happen to the message callback before we flush. + std::move(message_cb_cleanup).Invoke(); + term_hook_data->Flush(); + term_hook_data.reset(); + } + + double best_primal_bound = maximize ? -kInf : kInf; + switch (get_prim_stat(problem_)) { + case GLP_OPT: // OPT is returned by glp_ipt_status & glp_mip_status. + case GLP_FEAS: // FEAS is returned by glp_mip_status & glp_get_prim_stat. + best_primal_bound = obj_val(problem_); + break; + } + result.mutable_solve_stats()->set_best_primal_bound(best_primal_bound); + // TODO(b/187027049): compute the dual value when the dual is feasible (or + // problem optimal for interior point) based on the bounds and the dual values + // for LPs. + result.mutable_solve_stats()->set_best_dual_bound(best_dual_bound); + SolutionProto solution; + AddPrimalSolution(get_prim_stat, obj_val, col_val, model_parameters, + solution); + if (!is_mip) { + AddDualSolution(get_dual_stat, obj_val, row_dual, col_dual, + model_parameters, solution); + } + if (solution.has_primal_solution() || solution.has_dual_solution() || + solution.has_basis()) { + *result.add_solutions() = std::move(solution); + } + // TODO(b/200695800): add a parameter to enable the computation of the + // rays. This involves matrices inversion so this is not free to compute and + // should thus be only done when the user wants it. + RETURN_IF_ERROR(AddPrimalOrDualRay(model_parameters, result)); + + CHECK_OK(util_time::EncodeGoogleApiProto( + absl::Now() - start, result.mutable_solve_stats()->mutable_solve_time())); + return result; +} + +void GlpkSolver::AddVariables(const VariablesProto& new_variables) { + if (new_variables.ids().empty()) { + return; + } + + // Indices in GLPK are one-based. + const int first_new_var_index = variables_.ids.size() + 1; + + variables_.ids.insert(variables_.ids.end(), new_variables.ids().begin(), + new_variables.ids().end()); + for (int v = 0; v < new_variables.ids_size(); ++v) { + CHECK(variables_.id_to_index + .try_emplace(new_variables.ids(v), first_new_var_index + v) + .second); + } + glp_add_cols(problem_, new_variables.ids_size()); + if (!new_variables.names().empty()) { + for (int v = 0; v < new_variables.names_size(); ++v) { + glp_set_col_name( + problem_, v + first_new_var_index, + TruncateAndQuoteGLPKName(new_variables.names(v)).c_str()); + } + } + CHECK_EQ(new_variables.lower_bounds_size(), + new_variables.upper_bounds_size()); + CHECK_EQ(new_variables.lower_bounds_size(), new_variables.ids_size()); + variables_.unrounded_lower_bounds.insert( + variables_.unrounded_lower_bounds.end(), + new_variables.lower_bounds().begin(), new_variables.lower_bounds().end()); + variables_.unrounded_upper_bounds.insert( + variables_.unrounded_upper_bounds.end(), + new_variables.upper_bounds().begin(), new_variables.upper_bounds().end()); + for (int i = 0; i < new_variables.lower_bounds_size(); ++i) { + // Here we don't use the boolean "kind" GLP_BV since it does not exist. It + // is an artifact of glp_(get|set)_col_kind() functions. When + // glp_set_col_kind() is called with GLP_BV, in addition to setting the kind + // to GLP_IV (integer) it also sets the bounds to [0,1]. Symmetrically + // glp_get_col_kind() returns GLP_BV when the kind is GLP_IV and the bounds + // are [0,1]. + glp_set_col_kind(problem_, i + first_new_var_index, + new_variables.integers(i) ? GLP_IV : GLP_CV); + SetBounds(problem_, /*k=*/i + first_new_var_index, + {.lower = new_variables.lower_bounds(i), + .upper = new_variables.upper_bounds(i)}); + } +} + +void GlpkSolver::AddLinearConstraints( + const LinearConstraintsProto& new_linear_constraints) { + if (new_linear_constraints.ids().empty()) { + return; + } + + // Indices in GLPK are one-based. + const int first_new_cstr_index = linear_constraints_.ids.size() + 1; + + linear_constraints_.ids.insert(linear_constraints_.ids.end(), + new_linear_constraints.ids().begin(), + new_linear_constraints.ids().end()); + for (int c = 0; c < new_linear_constraints.ids_size(); ++c) { + CHECK(linear_constraints_.id_to_index + .try_emplace(new_linear_constraints.ids(c), + first_new_cstr_index + c) + .second); + } + glp_add_rows(problem_, new_linear_constraints.ids_size()); + if (!new_linear_constraints.names().empty()) { + for (int c = 0; c < new_linear_constraints.names_size(); ++c) { + glp_set_row_name( + problem_, c + first_new_cstr_index, + TruncateAndQuoteGLPKName(new_linear_constraints.names(c)).c_str()); + } + } + CHECK_EQ(new_linear_constraints.lower_bounds_size(), + new_linear_constraints.upper_bounds_size()); + for (int i = 0; i < new_linear_constraints.lower_bounds_size(); ++i) { + SetBounds( + problem_, /*k=*/i + first_new_cstr_index, + {.lower = new_linear_constraints.lower_bounds(i), + .upper = new_linear_constraints.upper_bounds(i)}); + } +} + +void GlpkSolver::UpdateObjectiveCoefficients( + const SparseDoubleVectorProto& coefficients_proto) { + for (const auto [id, coeff] : MakeView(coefficients_proto)) { + const int col_index = variables_.id_to_index.at(id); + CHECK_EQ(variables_.ids[col_index - 1], id); + glp_set_obj_coef(problem_, col_index, coeff); + } +} + +void GlpkSolver::UpdateLinearConstraintMatrix( + const SparseDoubleMatrixProto& matrix_updates, + const std::optional first_new_var_id, + const std::optional first_new_cstr_id) { + // GLPK's does not have an API to set matrix elements one by one. Instead it + // can either update an entire row or update an entire column or load the + // entire matrix. On top of that there is no API to get the entire matrix at + // once. + // + // Hence to update existing coefficients we have to read rows (or columns) + // coefficients, update existing non-zero that have been changed and add new + // values and write back the result. For new rows and columns we can be more + // efficient since we don't have to read the existing values back. + // + // The strategy used below is to split the matrix in three regions: + // + // existing new + // columns columns + // / | \ + // existing | 1 | 2 | + // rows | | | + // |---------+---------| + // new | | + // rows | 3 | + // \ / + // + // We start by updating the region 1 of existing rows and columns to limit the + // number of reads of existing coefficients. Then we update region 2 with all + // new columns but we only existing rows. Finally we update region 3 with all + // new rows and include new columns. Doing updates this way remove the need to + // read existing coefficients for the updates 2 & 3 since by construction + // those values are 0. + + // Updating existing rows (constraints), ignoring the new columns. + { + // We reuse the same vectors for all calls to GLPK's API to limit + // reallocations of these temporary buffers. + // + // We use constant size vectors to remove inefficiencies of + // std::vector::resize() that has linear cost in the size change (due to + // reset to 0 of values in the existing capacity on grow on resize). + // + // The value at index 0 is never used by GLPK's API (which are one-based). + std::vector data_indices(variables_.ids.size() + 1); + std::vector data_values(variables_.ids.size() + 1); + + // This shared vector (to prevent reallocation) will be used to store for + // each GLPK column the corresponding index in data_xxx vectors. The value + // kColNotPresent being used as a guard value to indicate that the column is + // not present. + // + // This uses the GLPK convention of ignoring the element 0 and using + // one-based indices. + constexpr int kColNotPresent = std::numeric_limits::max(); + std::vector col_to_data_element(variables_.ids.size() + 1, + kColNotPresent); + for (const auto& [row_id, row_coefficients] : + SparseSubmatrixByRows(matrix_updates, + /*start_row_id=*/0, + /*end_row_id=*/first_new_cstr_id, + /*start_col_id=*/0, + /*end_col_id=*/first_new_var_id)) { + // Find the index of the row in GLPK corresponding to the MathOpt's row + // id. + const int row_index = linear_constraints_.id_to_index.at(row_id); + CHECK_EQ(linear_constraints_.ids[row_index - 1], row_id); + + // Read the current row coefficients. + const int initial_non_zeros = glp_get_mat_row( + problem_, row_index, data_indices.data(), data_values.data()); + CHECK_LE(initial_non_zeros + 1, data_indices.size()); + CHECK_LE(initial_non_zeros + 1, data_values.size()); + + // Update the col to data_xxx elements map. + for (int i = 1; i <= initial_non_zeros; ++i) { + col_to_data_element[data_indices[i]] = i; + } + + // Update the row data. + int non_zeros = initial_non_zeros; + for (const auto [col_id, coefficient] : row_coefficients) { + const int col_index = variables_.id_to_index.at(col_id); + CHECK_EQ(variables_.ids[col_index - 1], col_id); + + // Here there are two cases: either a coefficient already exists for the + // given column, and we simply replace its value (potentially with a + // zero, which will remove the coefficient in the GLPK matrix), or we + // need to append it. + if (const int i = col_to_data_element[col_index]; i == kColNotPresent) { + ++non_zeros; + data_indices[non_zeros] = col_index; + data_values[non_zeros] = coefficient; + } else { + CHECK_EQ(data_indices[i], col_index); + data_values[i] = coefficient; + } + } + CHECK_LE(non_zeros + 1, data_indices.size()); + CHECK_LE(non_zeros + 1, data_values.size()); + + // Change the row values. + glp_set_mat_row(problem_, row_index, non_zeros, data_indices.data(), + data_values.data()); + + // Cleanup the used data_xxx items to make sure we don't reuse those + // values by mistake later. + for (int i = 1; i <= non_zeros; ++i) { + data_indices[i] = 0; + data_values[i] = 0; + } + + // Resets the elements of the map we modified. Here we ignore the new + // column indices that we have added in data_indices since those have + // not been put in the map. + for (int i = 1; i <= initial_non_zeros; ++i) { + col_to_data_element[data_indices[i]] = kColNotPresent; + } + } + } + + // Add new columns's coefficients of existing rows. The coefficients of new + // columns in new rows will be added when adding new rows below. + if (first_new_var_id.has_value()) { + // See the documentation for the existing rows above of the variables with + // the same names. + std::vector data_indices(linear_constraints_.ids.size() + 1); + std::vector data_values(linear_constraints_.ids.size() + 1); + for (const auto& [col_id, col_coefficients] : TransposeSparseSubmatrix( + SparseSubmatrixByRows(matrix_updates, + /*start_row_id=*/0, + /*end_row_id=*/first_new_cstr_id, + /*start_col_id=*/*first_new_var_id, + /*end_col_id=*/std::nullopt))) { + // Find the index of the column in GLPK corresponding to the MathOpt's + // column id. + const int col_index = variables_.id_to_index.at(col_id); + CHECK_EQ(variables_.ids[col_index - 1], col_id); + + // Prepare the column data replacing MathOpt ids by GLPK one-based row + // indices. + int non_zeros = 0; + for (const auto [row_id, coefficient] : MakeView(col_coefficients)) { + const int row_index = linear_constraints_.id_to_index.at(row_id); + CHECK_EQ(linear_constraints_.ids[row_index - 1], row_id); + + ++non_zeros; + data_indices[non_zeros] = row_index; + data_values[non_zeros] = coefficient; + } + CHECK_LE(non_zeros + 1, data_indices.size()); + CHECK_LE(non_zeros + 1, data_values.size()); + + // Change the column values. + glp_set_mat_col(problem_, col_index, non_zeros, data_indices.data(), + data_values.data()); + + // Cleanup the used data_xxx items to make sure we don't reuse those + // values by mistake later. + for (int i = 1; i <= non_zeros; ++i) { + data_indices[i] = 0; + data_values[i] = 0; + } + } + } + + // Add new rows, including the new columns' coefficients. + if (first_new_cstr_id.has_value()) { + // See the documentation for the existing rows above of the variables with + // the same names. + std::vector data_indices(variables_.ids.size() + 1); + std::vector data_values(variables_.ids.size() + 1); + for (const auto& [row_id, row_coefficients] : + SparseSubmatrixByRows(matrix_updates, + /*start_row_id=*/*first_new_cstr_id, + /*end_row_id=*/std::nullopt, + /*start_col_id=*/0, + /*end_col_id=*/std::nullopt)) { + // Find the index of the row in GLPK corresponding to the MathOpt's row + // id. + const int row_index = linear_constraints_.id_to_index.at(row_id); + CHECK_EQ(linear_constraints_.ids[row_index - 1], row_id); + + // Prepare the row data replacing MathOpt ids by GLPK one-based column + // indices. + int non_zeros = 0; + for (const auto [col_id, coefficient] : row_coefficients) { + const int col_index = variables_.id_to_index.at(col_id); + CHECK_EQ(variables_.ids[col_index - 1], col_id); + + ++non_zeros; + data_indices[non_zeros] = col_index; + data_values[non_zeros] = coefficient; + } + CHECK_LE(non_zeros + 1, data_indices.size()); + CHECK_LE(non_zeros + 1, data_values.size()); + + // Change the row values. + glp_set_mat_row(problem_, row_index, non_zeros, data_indices.data(), + data_values.data()); + + // Cleanup the used data_xxx items to make sure we don't reuse those + // values by mistake later. + for (int i = 1; i <= non_zeros; ++i) { + data_indices[i] = 0; + data_values[i] = 0; + } + } + } +} + +void GlpkSolver::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) { + const int status = get_prim_stat(problem_); + if (status == GLP_OPT || status == GLP_FEAS) { + PrimalSolutionProto& primal_solution = + *solution_proto.mutable_primal_solution(); + primal_solution.set_objective_value(obj_val(problem_)); + primal_solution.set_feasibility_status(SOLUTION_STATUS_FEASIBLE); + *primal_solution.mutable_variable_values() = + FilteredVector(problem_, model_parameters.variable_values_filter(), + variables_.ids, col_val); + } +} + +void GlpkSolver::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) { + const int status = get_dual_stat(problem_); + if (status == GLP_OPT || status == GLP_FEAS) { + DualSolutionProto& dual_solution = *solution_proto.mutable_dual_solution(); + dual_solution.set_objective_value(obj_val(problem_)); + *dual_solution.mutable_dual_values() = + FilteredVector(problem_, model_parameters.dual_values_filter(), + linear_constraints_.ids, row_dual); + *dual_solution.mutable_reduced_costs() = + FilteredVector(problem_, model_parameters.reduced_costs_filter(), + variables_.ids, col_dual); + // TODO(b/197867442): Check that `status == GLP_FEAS` implies dual feasible + // solution on early termination with barrier (where both `get_dual_stat` + // and `get_prim_stat` are equal to `glp_ipt_status`). + dual_solution.set_feasibility_status(SOLUTION_STATUS_FEASIBLE); + } +} + +absl::Status GlpkSolver::AddPrimalOrDualRay( + const ModelSolveParametersProto& model_parameters, + SolveResultProto& result) { + ASSIGN_OR_RETURN(const std::optional opt_unbound_ray, + GlpkComputeUnboundRay(problem_)); + if (!opt_unbound_ray.has_value()) { + return absl::OkStatus(); + } + + const int num_cstrs = linear_constraints_.ids.size(); + switch (opt_unbound_ray->type) { + case GlpkRayType::kPrimal: { + const int num_cstrs = linear_constraints_.ids.size(); + // Note that GlpkComputeUnboundRay() returned ray considers the variables + // of the computational form. Thus it contains both structural and + // auxiliary variables. In the MathOpt's primal ray we only consider + // structural variables though. + std::vector ray_values(variables_.ids.size()); + + for (const auto [k, value] : opt_unbound_ray->non_zero_components) { + if (k <= num_cstrs) { + // Ignore auxiliary variables. + continue; + } + const int var_index = k - num_cstrs; + CHECK_GE(var_index, 1); + ray_values[var_index - 1] = value; + } + + *result.add_primal_rays()->mutable_variable_values() = + FilteredRay(model_parameters.variable_values_filter(), variables_.ids, + ray_values); + + return absl::OkStatus(); + } + case GlpkRayType::kDual: { + // Note that GlpkComputeUnboundRay() returned ray considers the variables + // of the computational form. Thus it contains reduced costs of both + // structural and auxiliary variables. In the MathOpt's dual ray we split + // the reduced costs. The ones of auxiliary variables (variables of + // constraints) are called "dual values" and the ones of structural + // variables are called "reduced costs". + std::vector ray_reduced_costs(variables_.ids.size()); + std::vector ray_dual_values(num_cstrs); + + for (const auto [k, value] : opt_unbound_ray->non_zero_components) { + if (k <= num_cstrs) { + ray_dual_values[k - 1] = value; + } else { + const int var_index = k - num_cstrs; + CHECK_GE(var_index, 1); + ray_reduced_costs[var_index - 1] = value; + } + } + + DualRayProto& dual_ray = *result.add_dual_rays(); + *dual_ray.mutable_dual_values() = + FilteredRay(model_parameters.dual_values_filter(), + linear_constraints_.ids, ray_dual_values); + *dual_ray.mutable_reduced_costs() = + FilteredRay(model_parameters.reduced_costs_filter(), variables_.ids, + ray_reduced_costs); + + return absl::OkStatus(); + } + } +} + +absl::Status GlpkSolver::Update(const ModelUpdateProto& model_update) { + if (!model_update.objective_updates() + .quadratic_coefficients() + .row_ids() + .empty()) { + return QuadraticObjectiveError(); + } + + // TODO(b/187027049): GLPK should not support modifying the model from another + // thread (the allocation depends on the per-thread environment). We should + // unit test that and see what is the actual behavior. If GLPK itself does not + // provide its own assertion we should add one here. + { + const std::vector sorted_deleted_cols = DeleteRowsOrCols( + problem_, variables_, model_update.deleted_variable_ids()); + DeleteRowOrColData(variables_.unrounded_lower_bounds, sorted_deleted_cols); + DeleteRowOrColData(variables_.unrounded_upper_bounds, sorted_deleted_cols); + CHECK_EQ(variables_.unrounded_lower_bounds.size(), + variables_.unrounded_upper_bounds.size()); + CHECK_EQ(variables_.unrounded_lower_bounds.size(), variables_.ids.size()); + } + DeleteRowsOrCols(problem_, linear_constraints_, + model_update.deleted_linear_constraint_ids()); + + for (const auto [var_id, is_integer] : + MakeView(model_update.variable_updates().integers())) { + // See comment in AddVariables() to see why we don't use GLP_BV here. + const int var_index = variables_.id_to_index.at(var_id); + glp_set_col_kind(problem_, var_index, is_integer ? GLP_IV : GLP_CV); + + // Either restore the fractional bounds if the variable was integer and is + // now integer, or rounds the existing bounds if the variable was fractional + // and is now integer. Here we use the old bounds; they will get updated + // below by the call to UpdateBounds() if they are also changed by this + // update. + SetBounds( + problem_, var_index, + {.lower = variables_.unrounded_lower_bounds[var_index - 1], + .upper = variables_.unrounded_upper_bounds[var_index - 1]}); + } + for (const auto [var_id, lower_bound] : + MakeView(model_update.variable_updates().lower_bounds())) { + variables_.unrounded_lower_bounds[variables_.id_to_index.at(var_id) - 1] = + lower_bound; + } + for (const auto [var_id, upper_bound] : + MakeView(model_update.variable_updates().upper_bounds())) { + variables_.unrounded_upper_bounds[variables_.id_to_index.at(var_id) - 1] = + upper_bound; + } + UpdateBounds( + problem_, variables_, + /*lower_bounds_proto=*/model_update.variable_updates().lower_bounds(), + /*upper_bounds_proto=*/model_update.variable_updates().upper_bounds()); + UpdateBounds(problem_, linear_constraints_, + /*lower_bounds_proto=*/ + model_update.linear_constraint_updates().lower_bounds(), + /*upper_bounds_proto=*/ + model_update.linear_constraint_updates().upper_bounds()); + + AddVariables(model_update.new_variables()); + AddLinearConstraints(model_update.new_linear_constraints()); + + if (model_update.objective_updates().has_direction_update()) { + glp_set_obj_dir(problem_, + model_update.objective_updates().direction_update() + ? GLP_MAX + : GLP_MIN); + } + if (model_update.objective_updates().has_offset_update()) { + // Glpk uses index 0 for the "shift" of the objective. + glp_set_obj_coef(problem_, 0, + model_update.objective_updates().offset_update()); + } + UpdateObjectiveCoefficients( + model_update.objective_updates().linear_coefficients()); + + UpdateLinearConstraintMatrix( + model_update.linear_constraint_matrix_updates(), + /*first_new_var_id=*/FirstVariableId(model_update.new_variables()), + /*first_new_cstr_id=*/ + FirstLinearConstraintId(model_update.new_linear_constraints())); + + return absl::OkStatus(); +} + +bool GlpkSolver::CanUpdate(const ModelUpdateProto& model_update) { + // We return true even if we have a quadratic objective so that we don't force + // the caller to create a full model to get the error that quadratic + // objectives are not supported. The caller will get the correct error + // returned by GlpkSolver::Update(). + return true; +} + +MATH_OPT_REGISTER_SOLVER(SOLVER_TYPE_GLPK, GlpkSolver::New) + +} // namespace math_opt +} // namespace operations_research diff --git a/ortools/math_opt/solvers/glpk_solver.h b/ortools/math_opt/solvers/glpk_solver.h new file mode 100644 index 0000000000..80555b35ff --- /dev/null +++ b/ortools/math_opt/solvers/glpk_solver.h @@ -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 +#include +#include +#include + +#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 +} + +namespace operations_research { +namespace math_opt { + +class GlpkSolver : public SolverInterface { + public: + static absl::StatusOr> New( + const ModelProto& model, const InitArgs& init_args); + + ~GlpkSolver() override; + + absl::StatusOr 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 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 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 unrounded_lower_bounds; + + // The unrounded upper bound value of each column. + // + // See unrounded_lower_bounds documentation for details.. + std::vector 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 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 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 first_new_var_id, + std::optional 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_ diff --git a/ortools/math_opt/solvers/gscip_solver.cc b/ortools/math_opt/solvers/gscip_solver.cc index 814ba79a11..88e246224a 100644 --- a/ortools/math_opt/solvers/gscip_solver.cc +++ b/ortools/math_opt/solvers/gscip_solver.cc @@ -122,8 +122,8 @@ absl::flat_hash_map 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 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 GScipSolver::CreateSolveResultProto( GScipResult gscip_result, const ModelSolveParametersProto& model_parameters, const std::optional 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 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 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 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( diff --git a/ortools/math_opt/solvers/gscip_solver.h b/ortools/math_opt/solvers/gscip_solver.h index f8cb3a4488..485540393c 100644 --- a/ortools/math_opt/solvers/gscip_solver.h +++ b/ortools/math_opt/solvers/gscip_solver.h @@ -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 cutoff); + // Returns the ids of variables and linear constraints with inverted bounds. + InvertedBounds ListInvertedBounds() const; + const std::unique_ptr gscip_; InterruptEventHandler interrupt_event_handler_; absl::flat_hash_map variables_; diff --git a/ortools/math_opt/solvers/gurobi.proto b/ortools/math_opt/solvers/gurobi.proto index a9b0368e43..7f2640305b 100644 --- a/ortools/math_opt/solvers/gurobi.proto +++ b/ortools/math_opt/solvers/gurobi.proto @@ -21,7 +21,7 @@ message GurobiInitializerProto { message ISVKey { string name = 1; string application_name = 2; - int64 expiration = 3; + int32 expiration = 3; string key = 4; } diff --git a/ortools/math_opt/solvers/gurobi/g_gurobi.h b/ortools/math_opt/solvers/gurobi/g_gurobi.h index 430080e5f4..8007a76be3 100644 --- a/ortools/math_opt/solvers/gurobi/g_gurobi.h +++ b/ortools/math_opt/solvers/gurobi/g_gurobi.h @@ -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; }; diff --git a/ortools/math_opt/solvers/gurobi_callback.cc b/ortools/math_opt/solvers/gurobi_callback.cc index cb756b335c..c72ad93aab 100644 --- a/ortools/math_opt/solvers/gurobi_callback.cc +++ b/ortools/math_opt/solvers/gurobi_callback.cc @@ -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::infinity(); template diff --git a/ortools/math_opt/solvers/gurobi_solver.cc b/ortools/math_opt/solvers/gurobi_solver.cc index e4bd6d93a9..bf47ab07b4 100644 --- a/ortools/math_opt/solvers/gurobi_solver.cc +++ b/ortools/math_opt/solvers/gurobi_solver.cc @@ -68,8 +68,6 @@ namespace operations_research { namespace math_opt { namespace { -constexpr int kGrbOk = 0; - absl::StatusOr> 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 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 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 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 value; @@ -1860,6 +1869,34 @@ GurobiSolver::RegisterCallback(const CallbackRegistrationProto& registration, local_interrupter); } +absl::StatusOr GurobiSolver::ListInvertedBounds() const { + InvertedBounds inverted_bounds; + { + ASSIGN_OR_RETURN( + const std::vector var_lbs, + gurobi_->GetDoubleAttrArray(GRB_DBL_ATTR_LB, num_gurobi_variables_)); + ASSIGN_OR_RETURN( + const std::vector 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 GurobiSolver::Solve( const SolveParametersProto& parameters, const ModelSolveParametersProto& model_parameters, @@ -1934,6 +1971,14 @@ absl::StatusOr 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 diff --git a/ortools/math_opt/solvers/gurobi_solver.h b/ortools/math_opt/solvers/gurobi_solver.h index 82ee086fb8..6e72773f6a 100644 --- a/ortools/math_opt/solvers/gurobi_solver.h +++ b/ortools/math_opt/solvers/gurobi_solver.h @@ -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 ListInvertedBounds() const; + const std::unique_ptr gurobi_; // Note that we use linked_hash_map because the index of the gurobi_model_ diff --git a/ortools/math_opt/solvers/pdlp_bridge.cc b/ortools/math_opt/solvers/pdlp_bridge.cc index 1c08cff07b..4b1b5d1e20 100644 --- a/ortools/math_opt/solvers/pdlp_bridge.cc +++ b/ortools/math_opt/solvers/pdlp_bridge.cc @@ -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::FromProto( } const SparseDoubleMatrixProto& quadratic_objective = model_proto.objective().quadratic_coefficients(); - std::vector> 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::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 PdlpBridge::PrimalVariablesToProto( const Eigen::VectorXd& primal_values, const SparseVectorFilterProto& variable_filter) const { diff --git a/ortools/math_opt/solvers/pdlp_bridge.h b/ortools/math_opt/solvers/pdlp_bridge.h index 7293d34903..a72f8f9d02 100644 --- a/ortools/math_opt/solvers/pdlp_bridge.h +++ b/ortools/math_opt/solvers/pdlp_bridge.h @@ -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 PrimalVariablesToProto( diff --git a/ortools/math_opt/solvers/pdlp_solver.cc b/ortools/math_opt/solvers/pdlp_solver.cc index 46144e09a7..55e7a940a1 100644 --- a/ortools/math_opt/solvers/pdlp_solver.cc +++ b/ortools/math_opt/solvers/pdlp_solver.cc @@ -13,6 +13,7 @@ #include "ortools/math_opt/solvers/pdlp_solver.h" +#include #include #include #include @@ -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 PdlpSolver::MergeParameters( PrimalDualHybridGradientParams result; std::vector 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 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 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 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 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 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 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) { diff --git a/ortools/math_opt/tools/BUILD.bazel b/ortools/math_opt/tools/BUILD.bazel new file mode 100644 index 0000000000..793973e2ba --- /dev/null +++ b/ortools/math_opt/tools/BUILD.bazel @@ -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", + ], +) diff --git a/ortools/math_opt/tools/mathopt_solve_main.cc b/ortools/math_opt/tools/mathopt_solve_main.cc new file mode 100644 index 0000000000..fbfb422f8f --- /dev/null +++ b/ortools/math_opt/tools/mathopt_solve_main.cc @@ -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 +#include +#include +#include +#include +#include + +#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, 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 FormatFromFilePath( + const absl::string_view file_path) { + const std::vector> + 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 ReadModel(const absl::string_view file_path, + const absl::string_view format) { + if (format == kMathOptBinaryFormat) { + return file::GetBinaryProto(file_path, file::Defaults()); + } + if (format == kMathOptTextFormat) { + return file::GetTextProto(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 ReadModelUpdate( + const absl::string_view file_path, const absl::string_view format) { + if (format == kMathOptBinaryFormat) { + return file::GetBinaryProto(file_path, file::Defaults()); + } + if (format == kMathOptTextFormat) { + return file::GetTextProto(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 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 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 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::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; +} diff --git a/ortools/math_opt/validators/BUILD.bazel b/ortools/math_opt/validators/BUILD.bazel index 069965c4f7..5aeca87066 100644 --- a/ortools/math_opt/validators/BUILD.bazel +++ b/ortools/math_opt/validators/BUILD.bazel @@ -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", diff --git a/ortools/math_opt/validators/ids_validator.cc b/ortools/math_opt/validators/ids_validator.cc index 2f1b910e71..92dd9c9a98 100644 --- a/ortools/math_opt/validators/ids_validator.cc +++ b/ortools/math_opt/validators/ids_validator.cc @@ -33,8 +33,6 @@ namespace operations_research { namespace math_opt { -constexpr double kInf = std::numeric_limits::infinity(); - namespace { absl::Status CheckSortedIdsSubsetWithIndexOffset( const absl::Span ids, diff --git a/ortools/math_opt/validators/model_validator.cc b/ortools/math_opt/validators/model_validator.cc index e1d878f450..8baaac71ba 100644 --- a/ortools/math_opt/validators/model_validator.cc +++ b/ortools/math_opt/validators/model_validator.cc @@ -38,8 +38,6 @@ namespace operations_research { namespace math_opt { namespace { -constexpr double kInf = std::numeric_limits::infinity(); - //////////////////////////////////////////////////////////////////////////////// // Submessages //////////////////////////////////////////////////////////////////////////////// diff --git a/ortools/math_opt/validators/result_validator.cc b/ortools/math_opt/validators/result_validator.cc index 469cd03f03..b88bf463bc 100644 --- a/ortools/math_opt/validators/result_validator.cc +++ b/ortools/math_opt/validators/result_validator.cc @@ -31,18 +31,6 @@ namespace operations_research { namespace math_opt { namespace { -constexpr double kInf = std::numeric_limits::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"); diff --git a/ortools/math_opt/validators/solution_validator.cc b/ortools/math_opt/validators/solution_validator.cc index d8153d8246..50b797c0a5 100644 --- a/ortools/math_opt/validators/solution_validator.cc +++ b/ortools/math_opt/validators/solution_validator.cc @@ -34,8 +34,6 @@ namespace operations_research { namespace math_opt { namespace { -constexpr double kInf = std::numeric_limits::infinity(); - absl::Status ValidateSolutionStatus(const SolutionStatusProto& status) { if (!SolutionStatusProto_IsValid(status)) { return absl::InvalidArgumentError(absl::StrCat("status = ", status)); diff --git a/ortools/math_opt/validators/solve_parameters_validator.cc b/ortools/math_opt/validators/solve_parameters_validator.cc index cc30fcbe09..b865597f83 100644 --- a/ortools/math_opt/validators/solve_parameters_validator.cc +++ b/ortools/math_opt/validators/solve_parameters_validator.cc @@ -13,19 +13,33 @@ #include "ortools/math_opt/validators/solve_parameters_validator.h" +#include +#include +#include + #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() diff --git a/ortools/math_opt/validators/solve_stats_validator.cc b/ortools/math_opt/validators/solve_stats_validator.cc index f833945350..ca9f6a1c76 100644 --- a/ortools/math_opt/validators/solve_stats_validator.cc +++ b/ortools/math_opt/validators/solve_stats_validator.cc @@ -31,8 +31,6 @@ namespace operations_research { namespace math_opt { namespace { -constexpr double kInf = std::numeric_limits::infinity(); - absl::Status ValidateFeasibilityStatus(const FeasibilityStatusProto& status) { if (!FeasibilityStatusProto_IsValid(status)) { return absl::InvalidArgumentError(absl::StrCat("invalid status ", status));