From 1b35546a32f584be54ade8a83329bb0c95a1ddcb Mon Sep 17 00:00:00 2001 From: Corentin Le Molgat Date: Wed, 6 Dec 2023 17:33:45 +0100 Subject: [PATCH] math_opt: add pdlp support --- ortools/math_opt/BUILD.bazel | 1 + ortools/math_opt/core/python/BUILD.bazel | 4 + ortools/math_opt/cpp/parameters.cc | 9 +- ortools/math_opt/cpp/parameters.h | 8 + ortools/math_opt/parameters.proto | 9 +- ortools/math_opt/python/parameters.py | 5 + ortools/math_opt/samples/cpp/BUILD.bazel | 3 + .../math_opt/samples/cpp/linear_regression.cc | 2 +- ortools/math_opt/solvers/BUILD.bazel | 64 +++ ortools/math_opt/solvers/pdlp_bridge.cc | 222 +++++++++++ ortools/math_opt/solvers/pdlp_bridge.h | 85 ++++ ortools/math_opt/solvers/pdlp_solver.cc | 377 ++++++++++++++++++ ortools/math_opt/solvers/pdlp_solver.h | 70 ++++ ortools/math_opt/tools/BUILD.bazel | 1 + 14 files changed, 855 insertions(+), 5 deletions(-) create mode 100644 ortools/math_opt/solvers/pdlp_bridge.cc create mode 100644 ortools/math_opt/solvers/pdlp_bridge.h create mode 100644 ortools/math_opt/solvers/pdlp_solver.cc create mode 100644 ortools/math_opt/solvers/pdlp_solver.h diff --git a/ortools/math_opt/BUILD.bazel b/ortools/math_opt/BUILD.bazel index 3c91e4fd14..6d3b767195 100644 --- a/ortools/math_opt/BUILD.bazel +++ b/ortools/math_opt/BUILD.bazel @@ -109,6 +109,7 @@ proto_library( "//ortools/math_opt/solvers:gurobi_proto", "//ortools/math_opt/solvers:highs_proto", "//ortools/math_opt/solvers:osqp_proto", + "//ortools/pdlp:solvers_proto", "//ortools/sat:sat_parameters_proto", "@com_google_protobuf//:duration_proto", ], diff --git a/ortools/math_opt/core/python/BUILD.bazel b/ortools/math_opt/core/python/BUILD.bazel index a74312aaef..de24937390 100644 --- a/ortools/math_opt/core/python/BUILD.bazel +++ b/ortools/math_opt/core/python/BUILD.bazel @@ -34,6 +34,10 @@ pybind_extension( select({ "//ortools/linear_solver:use_scip": ["//ortools/math_opt/solvers:gscip_solver"], "//conditions:default": [], + }) + + select({ + "//ortools/linear_solver:use_pdlp": ["//ortools/math_opt/solvers:pdlp_solver"], + "//conditions:default": [], }) + [ "//ortools/math_opt:result_cc_proto", "//ortools/math_opt/core:solve_interrupter", diff --git a/ortools/math_opt/cpp/parameters.cc b/ortools/math_opt/cpp/parameters.cc index 3a9b350f79..a134b90507 100644 --- a/ortools/math_opt/cpp/parameters.cc +++ b/ortools/math_opt/cpp/parameters.cc @@ -73,6 +73,8 @@ std::optional Enum::ToOptString( return "glop"; case SolverType::kCpSat: return "cp_sat"; + case SolverType::kPdlp: + return "pdlp"; case SolverType::kGlpk: return "glpk"; case SolverType::kEcos: @@ -90,8 +92,9 @@ std::optional Enum::ToOptString( absl::Span Enum::AllValues() { static constexpr SolverType kSolverTypeValues[] = { SolverType::kGscip, SolverType::kGurobi, SolverType::kGlop, - SolverType::kCpSat, SolverType::kGlpk, SolverType::kEcos, - SolverType::kScs, SolverType::kHighs, SolverType::kSantorini, + SolverType::kCpSat, SolverType::kPdlp, SolverType::kGlpk, + SolverType::kEcos, SolverType::kScs, SolverType::kHighs, + SolverType::kSantorini, }; return absl::MakeConstSpan(kSolverTypeValues); } @@ -258,6 +261,7 @@ SolveParametersProto SolveParameters::Proto() const { *result.mutable_gurobi() = gurobi.Proto(); *result.mutable_glop() = glop; *result.mutable_cp_sat() = cp_sat; + *result.mutable_pdlp() = pdlp; *result.mutable_glpk() = glpk.Proto(); *result.mutable_highs() = highs; return result; @@ -316,6 +320,7 @@ absl::StatusOr SolveParameters::FromProto( result.gurobi = GurobiParameters::FromProto(proto.gurobi()); result.glop = proto.glop(); result.cp_sat = proto.cp_sat(); + result.pdlp = proto.pdlp(); result.glpk = GlpkParameters::FromProto(proto.glpk()); result.highs = proto.highs(); return result; diff --git a/ortools/math_opt/cpp/parameters.h b/ortools/math_opt/cpp/parameters.h index 77142f247e..0935084102 100644 --- a/ortools/math_opt/cpp/parameters.h +++ b/ortools/math_opt/cpp/parameters.h @@ -32,6 +32,7 @@ #include "ortools/math_opt/parameters.pb.h" #include "ortools/math_opt/solvers/gurobi.pb.h" // IWYU pragma: export #include "ortools/math_opt/solvers/highs.pb.h" // IWYU pragma: export +#include "ortools/pdlp/solvers.pb.h" // IWYU pragma: export #include "ortools/sat/sat_parameters.pb.h" // IWYU pragma: export namespace operations_research { @@ -63,6 +64,12 @@ enum class SolverType { // problems with continuous variables. kCpSat = SOLVER_TYPE_CP_SAT, + // Google's PDLP solver. + // + // Supports LP and convex diagonal quadratic objectives. Uses first order + // methods rather than simplex. Can solve very large problems. + kPdlp = SOLVER_TYPE_PDLP, + // GNU Linear Programming Kit (GLPK) (third party). // // Supports MIP and LP. @@ -409,6 +416,7 @@ struct SolveParameters { GurobiParameters gurobi; glop::GlopParameters glop; sat::SatParameters cp_sat; + pdlp::PrimalDualHybridGradientParams pdlp; GlpkParameters glpk; HighsOptionsProto highs; diff --git a/ortools/math_opt/parameters.proto b/ortools/math_opt/parameters.proto index dd683bbe47..4eca316ff6 100644 --- a/ortools/math_opt/parameters.proto +++ b/ortools/math_opt/parameters.proto @@ -17,6 +17,7 @@ syntax = "proto3"; package operations_research.math_opt; import "google/protobuf/duration.proto"; +import "ortools/pdlp/solvers.proto"; import "ortools/glop/parameters.proto"; import "ortools/gscip/gscip.proto"; import "ortools/math_opt/solvers/glpk.proto"; @@ -56,7 +57,11 @@ enum SolverTypeProto { // problems with continuous variables. SOLVER_TYPE_CP_SAT = 4; - reserved 5; + // Google's PDLP solver. + // + // Supports LP and convex diagonal quadratic objectives. Uses first order + // methods rather than simplex. Can solve very large problems. + SOLVER_TYPE_PDLP = 5; // GNU Linear Programming Kit (GLPK) (third party). // @@ -346,7 +351,7 @@ message SolveParametersProto { sat.SatParameters cp_sat = 15; - reserved 16; + pdlp.PrimalDualHybridGradientParams pdlp = 16; // Users should prefer the generic MathOpt parameters over OSQP-level // parameters, when available: diff --git a/ortools/math_opt/python/parameters.py b/ortools/math_opt/python/parameters.py index 002ac77528..7bdd4ec8fd 100644 --- a/ortools/math_opt/python/parameters.py +++ b/ortools/math_opt/python/parameters.py @@ -82,6 +82,7 @@ class SolverType(enum.Enum): GUROBI = math_opt_parameters_pb2.SOLVER_TYPE_GUROBI GLOP = math_opt_parameters_pb2.SOLVER_TYPE_GLOP CP_SAT = math_opt_parameters_pb2.SOLVER_TYPE_CP_SAT + PDLP = math_opt_parameters_pb2.SOLVER_TYPE_PDLP GLPK = math_opt_parameters_pb2.SOLVER_TYPE_GLPK OSQP = math_opt_parameters_pb2.SOLVER_TYPE_OSQP ECOS = math_opt_parameters_pb2.SOLVER_TYPE_ECOS @@ -403,6 +404,9 @@ class SolveParameters: cp_sat: sat_parameters_pb2.SatParameters = dataclasses.field( default_factory=sat_parameters_pb2.SatParameters ) + pdlp: pdlp_solvers_pb2.PrimalDualHybridGradientParams = dataclasses.field( + default_factory=pdlp_solvers_pb2.PrimalDualHybridGradientParams + ) osqp: osqp_pb2.OsqpSettingsProto = dataclasses.field( default_factory=osqp_pb2.OsqpSettingsProto ) @@ -424,6 +428,7 @@ class SolveParameters: gurobi=self.gurobi.to_proto(), glop=self.glop, cp_sat=self.cp_sat, + pdlp=self.pdlp, osqp=self.osqp, glpk=self.glpk.to_proto(), highs=self.highs, diff --git a/ortools/math_opt/samples/cpp/BUILD.bazel b/ortools/math_opt/samples/cpp/BUILD.bazel index 6b990b8c33..4708a5dbd3 100644 --- a/ortools/math_opt/samples/cpp/BUILD.bazel +++ b/ortools/math_opt/samples/cpp/BUILD.bazel @@ -108,6 +108,7 @@ cc_binary( "//ortools/math_opt/solvers:glop_solver", "//ortools/math_opt/solvers:glpk_solver", "//ortools/math_opt/solvers:gurobi_solver", + "//ortools/math_opt/solvers:pdlp_solver", "//ortools/util:status_macros", "@com_google_absl//absl/container:flat_hash_map", "@com_google_absl//absl/flags:flag", @@ -150,6 +151,7 @@ cc_binary( "//ortools/base:status_macros", "//ortools/math_opt/cpp:math_opt", "//ortools/math_opt/solvers:gurobi_solver", + "//ortools/math_opt/solvers:pdlp_solver", "@com_google_absl//absl/algorithm:container", "@com_google_absl//absl/random", "@com_google_absl//absl/status", @@ -168,6 +170,7 @@ cc_binary( "//ortools/math_opt/solvers:glpk_solver", "//ortools/math_opt/solvers:gscip_solver", "//ortools/math_opt/solvers:gurobi_solver", + "//ortools/math_opt/solvers:pdlp_solver", ], ) diff --git a/ortools/math_opt/samples/cpp/linear_regression.cc b/ortools/math_opt/samples/cpp/linear_regression.cc index 5dcdd82e87..00ead161a5 100644 --- a/ortools/math_opt/samples/cpp/linear_regression.cc +++ b/ortools/math_opt/samples/cpp/linear_regression.cc @@ -56,7 +56,7 @@ #include "ortools/math_opt/cpp/math_opt.h" ABSL_FLAG(operations_research::math_opt::SolverType, solver_type, - operations_research::math_opt::SolverType::kGurobi, + operations_research::math_opt::SolverType::kPdlp, "The solver needs to support quadratic objectives, e.g. pdlp, " "gurobi, or osqp."); diff --git a/ortools/math_opt/solvers/BUILD.bazel b/ortools/math_opt/solvers/BUILD.bazel index 6203c2df1e..8e49557aba 100644 --- a/ortools/math_opt/solvers/BUILD.bazel +++ b/ortools/math_opt/solvers/BUILD.bazel @@ -281,6 +281,70 @@ cc_library( ], ) +cc_library( + name = "pdlp_bridge", + srcs = ["pdlp_bridge.cc"], + hdrs = ["pdlp_bridge.h"], + deps = [ + "//ortools/base:status_macros", + "//ortools/math_opt:model_cc_proto", + "//ortools/math_opt:model_parameters_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:sparse_vector_view", + "//ortools/pdlp:primal_dual_hybrid_gradient", + "//ortools/pdlp:quadratic_program", + "@com_google_absl//absl/container:flat_hash_map", + "@com_google_absl//absl/status", + "@com_google_absl//absl/status:statusor", + "@com_google_absl//absl/strings", + "@eigen//:eigen3", + ], +) + +cc_library( + name = "pdlp_solver", + srcs = [ + "pdlp_solver.cc", + "pdlp_solver.h", + ], + visibility = ["//visibility:public"], + deps = [ + ":pdlp_bridge", + "//ortools/base:protoutil", + "//ortools/base:status_macros", + "//ortools/math_opt:callback_cc_proto", + "//ortools/math_opt:infeasible_subsystem_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/validators:callback_validator", + "//ortools/pdlp:iteration_stats", + "//ortools/pdlp:primal_dual_hybrid_gradient", + "//ortools/pdlp:quadratic_program", + "//ortools/pdlp:solve_log_cc_proto", + "//ortools/pdlp:solvers_cc_proto", + "//ortools/port:proto_utils", + "@com_google_absl//absl/log", + "@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/time", + ], + alwayslink = 1, +) + cc_library( name = "glpk_solver", srcs = [ diff --git a/ortools/math_opt/solvers/pdlp_bridge.cc b/ortools/math_opt/solvers/pdlp_bridge.cc new file mode 100644 index 0000000000..089724e345 --- /dev/null +++ b/ortools/math_opt/solvers/pdlp_bridge.cc @@ -0,0 +1,222 @@ +// Copyright 2010-2022 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/pdlp_bridge.h" + +#include +#include +#include +#include + +#include "Eigen/Core" +#include "Eigen/SparseCore" +#include "absl/container/flat_hash_map.h" +#include "absl/status/status.h" +#include "absl/status/statusor.h" +#include "absl/strings/str_cat.h" +#include "ortools/base/status_macros.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" +#include "ortools/math_opt/solution.pb.h" +#include "ortools/math_opt/sparse_containers.pb.h" +#include "ortools/pdlp/quadratic_program.h" + +namespace operations_research { +namespace math_opt { +namespace { + +constexpr SupportedProblemStructures kPdlpSupportedStructures = { + .quadratic_objectives = SupportType::kSupported}; + +absl::StatusOr ExtractSolution( + const Eigen::VectorXd& values, const std::vector& pdlp_index_to_id, + const SparseVectorFilterProto& filter, const double scale) { + if (values.size() != pdlp_index_to_id.size()) { + return absl::InternalError( + absl::StrCat("Expected solution vector with ", pdlp_index_to_id.size(), + " elements, found: ", values.size())); + } + SparseVectorFilterPredicate predicate(filter); + SparseDoubleVectorProto result; + for (int i = 0; i < pdlp_index_to_id.size(); ++i) { + const double value = scale * values[i]; + const int64_t id = pdlp_index_to_id[i]; + if (predicate.AcceptsAndUpdate(id, value)) { + result.add_ids(id); + result.add_values(value); + } + } + return result; +} + +// We are implicitly assuming that all missing IDs have correspoding value 0. +Eigen::VectorXd EncodeSolution( + const SparseDoubleVectorProto& values, + const absl::flat_hash_map& id_to_pdlp_index, + const double scale) { + Eigen::VectorXd pdlp_vector(Eigen::VectorXd::Zero(id_to_pdlp_index.size())); + const int num_values = values.values_size(); + for (int k = 0; k < num_values; ++k) { + const int64_t index = id_to_pdlp_index.at(values.ids(k)); + pdlp_vector[index] = values.values(k) / scale; + } + return pdlp_vector; +} + +} // namespace + +absl::StatusOr PdlpBridge::FromProto( + const ModelProto& model_proto) { + RETURN_IF_ERROR( + ModelIsSupported(model_proto, kPdlpSupportedStructures, "PDLP")); + PdlpBridge result; + pdlp::QuadraticProgram& pdlp_lp = result.pdlp_lp_; + const VariablesProto& variables = model_proto.variables(); + const LinearConstraintsProto& linear_constraints = + model_proto.linear_constraints(); + pdlp_lp.ResizeAndInitialize(variables.ids_size(), + linear_constraints.ids_size()); + if (!model_proto.name().empty()) { + pdlp_lp.problem_name = model_proto.name(); + } + if (variables.names_size() > 0) { + pdlp_lp.variable_names = {variables.names().begin(), + variables.names().end()}; + } + if (linear_constraints.names_size() > 0) { + pdlp_lp.constraint_names = {linear_constraints.names().begin(), + linear_constraints.names().end()}; + } + for (int i = 0; i < variables.ids_size(); ++i) { + result.var_id_to_pdlp_index_[variables.ids(i)] = i; + result.pdlp_index_to_var_id_.push_back(variables.ids(i)); + pdlp_lp.variable_lower_bounds[i] = variables.lower_bounds(i); + pdlp_lp.variable_upper_bounds[i] = variables.upper_bounds(i); + } + for (int i = 0; i < linear_constraints.ids_size(); ++i) { + result.lin_con_id_to_pdlp_index_[linear_constraints.ids(i)] = i; + result.pdlp_index_to_lin_con_id_.push_back(linear_constraints.ids(i)); + pdlp_lp.constraint_lower_bounds[i] = linear_constraints.lower_bounds(i); + pdlp_lp.constraint_upper_bounds[i] = linear_constraints.upper_bounds(i); + } + const bool is_maximize = model_proto.objective().maximize(); + const double obj_scale = is_maximize ? -1.0 : 1.0; + pdlp_lp.objective_offset = obj_scale * model_proto.objective().offset(); + for (const auto [var_id, coef] : + MakeView(model_proto.objective().linear_coefficients())) { + pdlp_lp.objective_vector[result.var_id_to_pdlp_index_.at(var_id)] = + obj_scale * coef; + } + const SparseDoubleMatrixProto& quadratic_objective = + model_proto.objective().quadratic_coefficients(); + 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, 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_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 + // method to do this transformation. + std::vector> mat_triplets; + const int nnz = model_proto.linear_constraint_matrix().row_ids_size(); + mat_triplets.reserve(nnz); + const SparseDoubleMatrixProto& proto_mat = + model_proto.linear_constraint_matrix(); + for (int i = 0; i < nnz; ++i) { + const int64_t row_index = + result.lin_con_id_to_pdlp_index_.at(proto_mat.row_ids(i)); + const int64_t column_index = + result.var_id_to_pdlp_index_.at(proto_mat.column_ids(i)); + const double value = proto_mat.coefficients(i); + mat_triplets.emplace_back(row_index, column_index, value); + } + pdlp_lp.constraint_matrix.setFromTriplets(mat_triplets.begin(), + mat_triplets.end()); + 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 { + return ExtractSolution(primal_values, pdlp_index_to_var_id_, variable_filter, + /*scale=*/1.0); +} +absl::StatusOr PdlpBridge::DualVariablesToProto( + const Eigen::VectorXd& dual_values, + const SparseVectorFilterProto& linear_constraint_filter) const { + return ExtractSolution(dual_values, pdlp_index_to_lin_con_id_, + linear_constraint_filter, + /*scale=*/pdlp_lp_.objective_scaling_factor); +} +absl::StatusOr PdlpBridge::ReducedCostsToProto( + const Eigen::VectorXd& reduced_costs, + const SparseVectorFilterProto& variable_filter) const { + return ExtractSolution(reduced_costs, pdlp_index_to_var_id_, variable_filter, + /*scale=*/pdlp_lp_.objective_scaling_factor); +} + +pdlp::PrimalAndDualSolution PdlpBridge::SolutionHintToWarmStart( + const SolutionHintProto& solution_hint) const { + // We are implicitly assuming that all missing IDs have correspoding value 0. + pdlp::PrimalAndDualSolution result; + result.primal_solution = EncodeSolution(solution_hint.variable_values(), + var_id_to_pdlp_index_, /*scale=*/1.0); + result.dual_solution = + EncodeSolution(solution_hint.dual_values(), lin_con_id_to_pdlp_index_, + /*scale=*/pdlp_lp_.objective_scaling_factor); + return result; +} + +} // namespace math_opt +} // namespace operations_research diff --git a/ortools/math_opt/solvers/pdlp_bridge.h b/ortools/math_opt/solvers/pdlp_bridge.h new file mode 100644 index 0000000000..22b299b577 --- /dev/null +++ b/ortools/math_opt/solvers/pdlp_bridge.h @@ -0,0 +1,85 @@ +// Copyright 2010-2022 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_PDLP_BRIDGE_H_ +#define OR_TOOLS_MATH_OPT_SOLVERS_PDLP_BRIDGE_H_ + +#include +#include +#include + +#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/model_parameters.pb.h" +#include "ortools/math_opt/sparse_containers.pb.h" +#include "ortools/pdlp/primal_dual_hybrid_gradient.h" +#include "ortools/pdlp/quadratic_program.h" + +namespace operations_research { +namespace math_opt { + +// Builds a PDLP model (QuadraticProgram) from ModelProto, and provides methods +// to translate solutions back and forth. +// +// The primary difference in the models are: +// 1. PDLP maps the variable/constraint ids to consecutive indices +// [0, 1, ..., n). +// 2. PDLP does not support maximization. If the ModelProto is a maximization +// problem, the objective is negated (coefficients and offset) before +// passing to PDLP. On the way back, the objective value, and all dual +// variables/reduced costs (also for rays) must be negated. +// +// Throughout, it is assumed that the MathOpt protos have been validated, but +// no assumption is made on the PDLP output. Any Status errors resulting from +// invalid PDLP output use the status code kInternal. +class PdlpBridge { + public: + PdlpBridge() = default; + static absl::StatusOr FromProto(const ModelProto& model_proto); + + 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( + const Eigen::VectorXd& primal_values, + const SparseVectorFilterProto& variable_filter) const; + absl::StatusOr DualVariablesToProto( + const Eigen::VectorXd& dual_values, + const SparseVectorFilterProto& linear_constraint_filter) const; + absl::StatusOr ReducedCostsToProto( + const Eigen::VectorXd& reduced_costs, + const SparseVectorFilterProto& variable_filter) const; + pdlp::PrimalAndDualSolution SolutionHintToWarmStart( + const SolutionHintProto& solution_hint) const; + + private: + pdlp::QuadraticProgram pdlp_lp_; + absl::flat_hash_map var_id_to_pdlp_index_; + // NOTE: this vector is strictly increasing + std::vector pdlp_index_to_var_id_; + absl::flat_hash_map lin_con_id_to_pdlp_index_; + // NOTE: this vector is strictly increasing + std::vector pdlp_index_to_lin_con_id_; +}; + +} // namespace math_opt +} // namespace operations_research + +#endif // OR_TOOLS_MATH_OPT_SOLVERS_PDLP_BRIDGE_H_ diff --git a/ortools/math_opt/solvers/pdlp_solver.cc b/ortools/math_opt/solvers/pdlp_solver.cc new file mode 100644 index 0000000000..03d8cf0679 --- /dev/null +++ b/ortools/math_opt/solvers/pdlp_solver.cc @@ -0,0 +1,377 @@ +// Copyright 2010-2022 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/pdlp_solver.h" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#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/time/time.h" +#include "ortools/base/logging.h" +#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" +#include "ortools/math_opt/infeasible_subsystem.pb.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/pdlp_bridge.h" +#include "ortools/math_opt/sparse_containers.pb.h" +#include "ortools/math_opt/validators/callback_validator.h" +#include "ortools/pdlp/iteration_stats.h" +#include "ortools/pdlp/primal_dual_hybrid_gradient.h" +#include "ortools/pdlp/quadratic_program.h" +#include "ortools/pdlp/solve_log.pb.h" +#include "ortools/pdlp/solvers.pb.h" +#include "ortools/port/proto_utils.h" + +namespace operations_research { +namespace math_opt { + +using pdlp::PrimalAndDualSolution; +using pdlp::PrimalDualHybridGradientParams; +using pdlp::SolverResult; + +absl::StatusOr> PdlpSolver::New( + const ModelProto& model, const InitArgs&) { + auto result = absl::WrapUnique(new PdlpSolver); + ASSIGN_OR_RETURN(result->pdlp_bridge_, PdlpBridge::FromProto(model)); + return result; +} + +absl::StatusOr PdlpSolver::MergeParameters( + const SolveParametersProto& parameters, const bool has_message_callback) { + PrimalDualHybridGradientParams result; + std::vector warnings; + if (parameters.enable_output() && !has_message_callback) { + result.set_verbosity_level(3); + } + if (parameters.has_threads()) { + result.set_num_threads(parameters.threads()); + } + if (parameters.has_time_limit()) { + result.mutable_termination_criteria()->set_time_sec_limit( + 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"); + } + if (parameters.has_objective_limit()) { + warnings.push_back("parameter best_objective_limit not supported for PDLP"); + } + if (parameters.has_best_bound_limit()) { + warnings.push_back("parameter best_bound_limit not supported for PDLP"); + } + if (parameters.has_solution_limit()) { + warnings.push_back("parameter solution_limit not supported for PDLP"); + } + if (parameters.has_random_seed()) { + warnings.push_back("parameter random_seed not supported for PDLP"); + } + if (parameters.lp_algorithm() != LP_ALGORITHM_UNSPECIFIED && + parameters.lp_algorithm() != LP_ALGORITHM_FIRST_ORDER) { + warnings.push_back("parameter lp_algorithm not supported for PDLP"); + } + if (parameters.presolve() != EMPHASIS_UNSPECIFIED) { + warnings.push_back("parameter presolve not supported for PDLP"); + } + if (parameters.cuts() != EMPHASIS_UNSPECIFIED) { + warnings.push_back("parameter cuts not supported for PDLP"); + } + if (parameters.heuristics() != EMPHASIS_UNSPECIFIED) { + warnings.push_back("parameter heuristics not supported for PDLP"); + } + if (parameters.scaling() != EMPHASIS_UNSPECIFIED) { + warnings.push_back("parameter scaling not supported for PDLP"); + } + if (parameters.has_iteration_limit()) { + const int64_t limit = std::min(std::numeric_limits::max(), + parameters.iteration_limit()); + result.mutable_termination_criteria()->set_iteration_limit( + static_cast(limit)); + } + result.MergeFrom(parameters.pdlp()); + if (!warnings.empty()) { + return absl::InvalidArgumentError(absl::StrJoin(warnings, "; ")); + } + return result; +} + +namespace { + +absl::StatusOr ConvertReason( + const pdlp::TerminationReason pdlp_reason, absl::string_view pdlp_detail, + const bool is_maximize, const double primal_objective, + const double dual_objective) { + switch (pdlp_reason) { + case pdlp::TERMINATION_REASON_UNSPECIFIED: + return TerminateForReason( + is_maximize, TERMINATION_REASON_OTHER_ERROR, + absl::StrCat("PDLP unspecified reason " + "(TERMINATION_REASON_UNSPECIFIED): ", + pdlp_detail)); + return TerminateForReason(is_maximize, TERMINATION_REASON_UNSPECIFIED, + pdlp_detail); + case pdlp::TERMINATION_REASON_OPTIMAL: + return OptimalTerminationProto(primal_objective, dual_objective, + pdlp_detail); + case pdlp::TERMINATION_REASON_PRIMAL_INFEASIBLE: + return InfeasibleTerminationProto( + is_maximize, + /*dual_feasibility_status=*/FEASIBILITY_STATUS_UNDETERMINED, + pdlp_detail); + case pdlp::TERMINATION_REASON_DUAL_INFEASIBLE: + return InfeasibleOrUnboundedTerminationProto( + is_maximize, + /*dual_feasibility_status=*/FEASIBILITY_STATUS_INFEASIBLE, + pdlp_detail); + case pdlp::TERMINATION_REASON_TIME_LIMIT: + return NoSolutionFoundTerminationProto( + is_maximize, LIMIT_TIME, + /*optional_dual_objective=*/std::nullopt, pdlp_detail); + case pdlp::TERMINATION_REASON_ITERATION_LIMIT: + return NoSolutionFoundTerminationProto( + is_maximize, LIMIT_ITERATION, + /*optional_dual_objective=*/std::nullopt, pdlp_detail); + case pdlp::TERMINATION_REASON_KKT_MATRIX_PASS_LIMIT: + return NoSolutionFoundTerminationProto( + is_maximize, LIMIT_OTHER, + /*optional_dual_objective=*/std::nullopt, pdlp_detail); + case pdlp::TERMINATION_REASON_NUMERICAL_ERROR: + return TerminateForReason(is_maximize, TERMINATION_REASON_NUMERICAL_ERROR, + pdlp_detail); + case pdlp::TERMINATION_REASON_INTERRUPTED_BY_USER: + return NoSolutionFoundTerminationProto( + is_maximize, LIMIT_INTERRUPTED, + /*optional_dual_objective=*/std::nullopt, pdlp_detail); + case pdlp::TERMINATION_REASON_INVALID_PROBLEM: + // Indicates that the solver detected invalid problem data, e.g., + // inconsistent bounds. + return absl::InternalError( + absl::StrCat("Invalid problem sent to PDLP solver " + "(TERMINATION_REASON_INVALID_PROBLEM): ", + pdlp_detail)); + case pdlp::TERMINATION_REASON_INVALID_INITIAL_SOLUTION: + return absl::InvalidArgumentError( + absl::StrCat("PDLP solution hint invalid " + "(TERMINATION_REASON_INVALID_INITIAL_SOLUTION): ", + pdlp_detail)); + case pdlp::TERMINATION_REASON_INVALID_PARAMETER: + // Indicates that an invalid value for the parameters was detected. + return absl::InvalidArgumentError(absl::StrCat( + "PDLP parameters invalid (TERMINATION_REASON_INVALID_PARAMETER): ", + pdlp_detail)); + case pdlp::TERMINATION_REASON_PRIMAL_OR_DUAL_INFEASIBLE: + return InfeasibleOrUnboundedTerminationProto( + is_maximize, + /*dual_feasibility_status=*/FEASIBILITY_STATUS_UNDETERMINED, + pdlp_detail); + case pdlp::TERMINATION_REASON_OTHER: + return TerminateForReason(is_maximize, TERMINATION_REASON_OTHER_ERROR, + pdlp_detail); + } + return absl::InvalidArgumentError(absl::StrCat( + "PDLP status: ", ProtoEnumToString(pdlp_reason), " not implemented.")); +} + +} // namespace + +absl::StatusOr PdlpSolver::MakeSolveResult( + const pdlp::SolverResult& pdlp_result, + const ModelSolveParametersProto& model_params) { + // PDLP's default is a minimization problem for which the default primal and + // dual bounds are infinity and -infinity respectively. PDLP provides a + // scaling factor to flip the signs for maximization problems. + const double objective_scaling_factor = + pdlp_bridge_.pdlp_lp().objective_scaling_factor; + const bool is_maximize = objective_scaling_factor < 0.0; + SolveResultProto result; + const std::optional convergence_information = + pdlp::GetConvergenceInformation(pdlp_result.solve_log.solution_stats(), + pdlp_result.solve_log.solution_type()); + if (convergence_information.has_value()) { + *result.mutable_pdlp_output()->mutable_convergence_information() = + *convergence_information; + } + ObjectiveBoundsProto objective_bounds = MakeTrivialBounds(is_maximize); + if (convergence_information.has_value()) { + objective_bounds.set_primal_bound( + convergence_information->primal_objective()); + objective_bounds.set_dual_bound(convergence_information->dual_objective()); + } + ASSIGN_OR_RETURN(*result.mutable_termination(), + ConvertReason(pdlp_result.solve_log.termination_reason(), + pdlp_result.solve_log.termination_string(), + is_maximize, objective_bounds.primal_bound(), + objective_bounds.dual_bound())); + ASSIGN_OR_RETURN(*result.mutable_solve_stats()->mutable_solve_time(), + util_time::EncodeGoogleApiProto( + absl::Seconds(pdlp_result.solve_log.solve_time_sec()))); + result.mutable_solve_stats()->set_first_order_iterations( + pdlp_result.solve_log.iteration_count()); + + switch (pdlp_result.solve_log.termination_reason()) { + case pdlp::TERMINATION_REASON_OPTIMAL: + 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_INTERRUPTED_BY_USER: { + SolutionProto* solution_proto = result.add_solutions(); + { + auto maybe_primal = pdlp_bridge_.PrimalVariablesToProto( + pdlp_result.primal_solution, model_params.variable_values_filter()); + RETURN_IF_ERROR(maybe_primal.status()); + PrimalSolutionProto* primal_proto = + solution_proto->mutable_primal_solution(); + primal_proto->set_feasibility_status(SOLUTION_STATUS_UNDETERMINED); + *primal_proto->mutable_variable_values() = *std::move(maybe_primal); + // Note: the solution could be primal feasible for termination reasons + // other than TERMINATION_REASON_OPTIMAL, but in theory, PDLP could + // also be modified to return the best feasible solution encountered in + // an early termination run (if any). + if (pdlp_result.solve_log.termination_reason() == + pdlp::TERMINATION_REASON_OPTIMAL) { + primal_proto->set_feasibility_status(SOLUTION_STATUS_FEASIBLE); + } + if (convergence_information.has_value()) { + primal_proto->set_objective_value( + convergence_information->primal_objective()); + } + } + { + auto maybe_dual = pdlp_bridge_.DualVariablesToProto( + pdlp_result.dual_solution, model_params.dual_values_filter()); + RETURN_IF_ERROR(maybe_dual.status()); + auto maybe_reduced = pdlp_bridge_.ReducedCostsToProto( + pdlp_result.reduced_costs, model_params.reduced_costs_filter()); + RETURN_IF_ERROR(maybe_reduced.status()); + DualSolutionProto* dual_proto = solution_proto->mutable_dual_solution(); + dual_proto->set_feasibility_status(SOLUTION_STATUS_UNDETERMINED); + *dual_proto->mutable_dual_values() = *std::move(maybe_dual); + *dual_proto->mutable_reduced_costs() = *std::move(maybe_reduced); + // Note: same comment on primal solution status holds here. + if (pdlp_result.solve_log.termination_reason() == + pdlp::TERMINATION_REASON_OPTIMAL) { + dual_proto->set_feasibility_status(SOLUTION_STATUS_FEASIBLE); + } + if (convergence_information.has_value()) { + dual_proto->set_objective_value( + convergence_information->dual_objective()); + } + } + break; + } + case pdlp::TERMINATION_REASON_PRIMAL_INFEASIBLE: { + // NOTE: for primal infeasible problems, PDLP stores the infeasibility + // certificate (dual ray) in the dual variables and reduced costs. + auto maybe_dual = pdlp_bridge_.DualVariablesToProto( + pdlp_result.dual_solution, model_params.dual_values_filter()); + RETURN_IF_ERROR(maybe_dual.status()); + auto maybe_reduced = pdlp_bridge_.ReducedCostsToProto( + pdlp_result.reduced_costs, model_params.reduced_costs_filter()); + RETURN_IF_ERROR(maybe_reduced.status()); + DualRayProto* dual_ray_proto = result.add_dual_rays(); + *dual_ray_proto->mutable_dual_values() = *std::move(maybe_dual); + *dual_ray_proto->mutable_reduced_costs() = *std::move(maybe_reduced); + break; + } + case pdlp::TERMINATION_REASON_DUAL_INFEASIBLE: { + // NOTE: for dual infeasible problems, PDLP stores the infeasibility + // certificate (primal ray) in the primal variables. + auto maybe_primal = pdlp_bridge_.PrimalVariablesToProto( + pdlp_result.primal_solution, model_params.variable_values_filter()); + RETURN_IF_ERROR(maybe_primal.status()); + PrimalRayProto* primal_ray_proto = result.add_primal_rays(); + *primal_ray_proto->mutable_variable_values() = *std::move(maybe_primal); + break; + } + default: + break; + } + return result; +} + +absl::StatusOr PdlpSolver::Solve( + const SolveParametersProto& parameters, + const ModelSolveParametersProto& model_parameters, + const MessageCallback message_cb, + const CallbackRegistrationProto& callback_registration, const Callback, + SolveInterrupter* const interrupter) { + RETURN_IF_ERROR(CheckRegisteredCallbackEvents(callback_registration, + /*supported_events=*/{})); + + ASSIGN_OR_RETURN( + auto pdlp_params, + MergeParameters(parameters, + /*has_message_callback=*/message_cb != nullptr)); + + // 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; }); + + std::optional initial_solution; + if (!model_parameters.solution_hints().empty()) { + initial_solution = pdlp_bridge_.SolutionHintToWarmStart( + model_parameters.solution_hints(0)); + } + + const SolverResult pdlp_result = PrimalDualHybridGradient( + pdlp_bridge_.pdlp_lp(), pdlp_params, initial_solution, &interrupt); + return MakeSolveResult(pdlp_result, model_parameters); +} + +absl::StatusOr PdlpSolver::Update(const ModelUpdateProto&) { + return false; +} + +absl::StatusOr +PdlpSolver::ComputeInfeasibleSubsystem(const SolveParametersProto&, + MessageCallback, + SolveInterrupter* const) { + return absl::UnimplementedError( + "PDLP does not provide a method to compute an infeasible subsystem"); +} + +MATH_OPT_REGISTER_SOLVER(SOLVER_TYPE_PDLP, PdlpSolver::New); + +} // namespace math_opt +} // namespace operations_research diff --git a/ortools/math_opt/solvers/pdlp_solver.h b/ortools/math_opt/solvers/pdlp_solver.h new file mode 100644 index 0000000000..f6533cef1b --- /dev/null +++ b/ortools/math_opt/solvers/pdlp_solver.h @@ -0,0 +1,70 @@ +// Copyright 2010-2022 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_PDLP_SOLVER_H_ +#define OR_TOOLS_MATH_OPT_SOLVERS_PDLP_SOLVER_H_ + +#include + +#include "absl/status/statusor.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/infeasible_subsystem.pb.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/solvers/pdlp_bridge.h" +#include "ortools/pdlp/primal_dual_hybrid_gradient.h" +#include "ortools/pdlp/solvers.pb.h" + +namespace operations_research { +namespace math_opt { + +class PdlpSolver : public SolverInterface { + public: + static absl::StatusOr> New( + const ModelProto& model, const InitArgs& init_args); + + absl::StatusOr Solve( + const SolveParametersProto& parameters, + const ModelSolveParametersProto& model_parameters, + MessageCallback message_cb, + const CallbackRegistrationProto& callback_registration, Callback cb, + SolveInterrupter* interrupter) override; + absl::StatusOr Update(const ModelUpdateProto& model_update) override; + absl::StatusOr + ComputeInfeasibleSubsystem(const SolveParametersProto& parameters, + MessageCallback message_cb, + SolveInterrupter* interrupter) override; + + // Returns the merged parameters and a list of warnings. + static absl::StatusOr MergeParameters( + const SolveParametersProto& parameters, bool has_message_callback); + + private: + PdlpSolver() = default; + + absl::StatusOr MakeSolveResult( + const pdlp::SolverResult& pdlp_result, + const ModelSolveParametersProto& model_params); + + PdlpBridge pdlp_bridge_; +}; + +} // namespace math_opt +} // namespace operations_research + +#endif // OR_TOOLS_MATH_OPT_SOLVERS_PDLP_SOLVER_H_ diff --git a/ortools/math_opt/tools/BUILD.bazel b/ortools/math_opt/tools/BUILD.bazel index f0b5ccfe97..93d94e76b1 100644 --- a/ortools/math_opt/tools/BUILD.bazel +++ b/ortools/math_opt/tools/BUILD.bazel @@ -33,6 +33,7 @@ cc_binary( "//ortools/math_opt/solvers:glpk_solver", "//ortools/math_opt/solvers:gscip_solver", "//ortools/math_opt/solvers:gurobi_solver", + "//ortools/math_opt/solvers:pdlp_solver", "//ortools/util:status_macros", "@com_google_absl//absl/flags:flag", "@com_google_absl//absl/log:check",