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

math_opt: add pdlp support

Browse Source
This commit is contained in:
Corentin Le Molgat
2023-12-06 17:33:45 +01:00
parent 47fcebdea6
commit 1b35546a32
14 changed files with 855 additions and 5 deletions
Download Patch File Download Diff File Expand all files Collapse all files

1
ortools/math_opt/BUILD.bazel
View File

@@ -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",
],

4
ortools/math_opt/core/python/BUILD.bazel
View File

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

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

@@ -73,6 +73,8 @@ std::optional<absl::string_view> Enum<SolverType>::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<absl::string_view> Enum<SolverType>::ToOptString(
absl::Span<const SolverType> Enum<SolverType>::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> 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;

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

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

9
ortools/math_opt/parameters.proto
View File

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

5
ortools/math_opt/python/parameters.py
View File

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

3
ortools/math_opt/samples/cpp/BUILD.bazel
View File

@@ -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",
],
)

2
ortools/math_opt/samples/cpp/linear_regression.cc
View File

@@ -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.");

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

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

222
ortools/math_opt/solvers/pdlp_bridge.cc Normal file
View File

@@ -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 <cstdint>
#include <optional>
#include <string>
#include <vector>
#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<SparseDoubleVectorProto> ExtractSolution(
const Eigen::VectorXd& values, const std::vector<int64_t>& 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<int64_t, int64_t>& 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> 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<Eigen::Triplet<double, int64_t>> 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<SparseDoubleVectorProto> 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<SparseDoubleVectorProto> 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<SparseDoubleVectorProto> 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

85
ortools/math_opt/solvers/pdlp_bridge.h Normal file
View File

@@ -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 <cstdint>
#include <optional>
#include <vector>
#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<PdlpBridge> 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<SparseDoubleVectorProto> PrimalVariablesToProto(
const Eigen::VectorXd& primal_values,
const SparseVectorFilterProto& variable_filter) const;
absl::StatusOr<SparseDoubleVectorProto> DualVariablesToProto(
const Eigen::VectorXd& dual_values,
const SparseVectorFilterProto& linear_constraint_filter) const;
absl::StatusOr<SparseDoubleVectorProto> 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<int64_t, int64_t> var_id_to_pdlp_index_;
// NOTE: this vector is strictly increasing
std::vector<int64_t> pdlp_index_to_var_id_;
absl::flat_hash_map<int64_t, int64_t> lin_con_id_to_pdlp_index_;
// NOTE: this vector is strictly increasing
std::vector<int64_t> pdlp_index_to_lin_con_id_;
};
} // namespace math_opt
} // namespace operations_research
#endif // OR_TOOLS_MATH_OPT_SOLVERS_PDLP_BRIDGE_H_

377
ortools/math_opt/solvers/pdlp_solver.cc Normal file
View File

@@ -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 <algorithm>
#include <atomic>
#include <cstdint>
#include <limits>
#include <memory>
#include <optional>
#include <string>
#include <utility>
#include <vector>
#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<std::unique_ptr<SolverInterface>> 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<PrimalDualHybridGradientParams> PdlpSolver::MergeParameters(
const SolveParametersProto& parameters, const bool has_message_callback) {
PrimalDualHybridGradientParams result;
std::vector<std::string> 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<int64_t>(std::numeric_limits<int32_t>::max(),
parameters.iteration_limit());
result.mutable_termination_criteria()->set_iteration_limit(
static_cast<int32_t>(limit));
}
result.MergeFrom(parameters.pdlp());
if (!warnings.empty()) {
return absl::InvalidArgumentError(absl::StrJoin(warnings, "; "));
}
return result;
}
namespace {
absl::StatusOr<TerminationProto> 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<SolveResultProto> 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<pdlp::ConvergenceInformation> 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<SolveResultProto> 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<bool> interrupt = false;
const ScopedSolveInterrupterCallback set_interrupt(
interrupter, [&]() { interrupt = true; });
std::optional<PrimalAndDualSolution> 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<bool> PdlpSolver::Update(const ModelUpdateProto&) {
return false;
}
absl::StatusOr<ComputeInfeasibleSubsystemResultProto>
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

70
ortools/math_opt/solvers/pdlp_solver.h Normal file
View File

@@ -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 <memory>
#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<std::unique_ptr<SolverInterface>> New(
const ModelProto& model, const InitArgs& init_args);
absl::StatusOr<SolveResultProto> Solve(
const SolveParametersProto& parameters,
const ModelSolveParametersProto& model_parameters,
MessageCallback message_cb,
const CallbackRegistrationProto& callback_registration, Callback cb,
SolveInterrupter* interrupter) override;
absl::StatusOr<bool> Update(const ModelUpdateProto& model_update) override;
absl::StatusOr<ComputeInfeasibleSubsystemResultProto>
ComputeInfeasibleSubsystem(const SolveParametersProto& parameters,
MessageCallback message_cb,
SolveInterrupter* interrupter) override;
// Returns the merged parameters and a list of warnings.
static absl::StatusOr<pdlp::PrimalDualHybridGradientParams> MergeParameters(
const SolveParametersProto& parameters, bool has_message_callback);
private:
PdlpSolver() = default;
absl::StatusOr<SolveResultProto> 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_

1
ortools/math_opt/tools/BUILD.bazel
View File

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