diff --git a/ortools/linear_solver/linear_solver.cc b/ortools/linear_solver/linear_solver.cc index 45f08accea..d735b70667 100644 --- a/ortools/linear_solver/linear_solver.cc +++ b/ortools/linear_solver/linear_solver.cc @@ -1395,6 +1395,8 @@ void MPSolver::SetStartingLpBasis( interface_->SetStartingLpBasis(variable_statuses, constraint_statuses); } +double MPSolver::solver_infinity() { return interface_->infinity(); } + MPVariable* MPSolver::MakeVar(double lb, double ub, bool integer, const std::string& name) { const int var_index = NumVariables(); diff --git a/ortools/linear_solver/linear_solver.h b/ortools/linear_solver/linear_solver.h index 522d91b506..6497c9da3a 100644 --- a/ortools/linear_solver/linear_solver.h +++ b/ortools/linear_solver/linear_solver.h @@ -717,6 +717,7 @@ class MPSolver { * You can use -MPSolver::infinity() for negative infinity. */ static double infinity() { return std::numeric_limits::infinity(); } + double solver_infinity(); /** * Controls (or queries) the amount of output produced by the underlying @@ -1764,6 +1765,8 @@ class MPSolverInterface { LOG(FATAL) << "Not supported by this solver."; } + virtual double infinity() { return std::numeric_limits::infinity(); } + virtual bool InterruptSolve() { return false; } // See MPSolver::NextSolution() for contract. diff --git a/ortools/linear_solver/proto_solver/scip_proto_solver.cc b/ortools/linear_solver/proto_solver/scip_proto_solver.cc index def833aa71..60b35311cd 100644 --- a/ortools/linear_solver/proto_solver/scip_proto_solver.cc +++ b/ortools/linear_solver/proto_solver/scip_proto_solver.cc @@ -37,6 +37,7 @@ #include "absl/time/time.h" #include "ortools/base/cleanup.h" #include "ortools/base/commandlineflags.h" +#include "ortools/base/logging.h" #include "ortools/base/status_macros.h" #include "ortools/base/timer.h" #include "ortools/gscip/legacy_scip_params.h" @@ -49,6 +50,7 @@ #include "scip/cons_quadratic.h" #include "scip/pub_var.h" #include "scip/scip.h" +#include "scip/scip_numerics.h" #include "scip/scip_param.h" #include "scip/scip_prob.h" #include "scip/scip_var.h" @@ -65,6 +67,22 @@ ABSL_FLAG(std::string, scip_proto_solver_output_cip_file, "", namespace operations_research { namespace { +// Clamp to SCIP's finite range, setting infinities to the SCIP infinities. +absl::StatusOr ScipInfClamp(SCIP* scip, const double d) { + if (std::isnan(d)) { + return absl::InvalidArgumentError("Cannot clamp a NaN."); + } + const double kScipInf = SCIPinfinity(scip); + if (d > -kScipInf && d < kScipInf) { + return d; + } + const double clamped = std::clamp(d, -kScipInf, kScipInf); + LOG_EVERY_N_SEC(WARNING, 5) << "A bound was clamped to SCIP's 'infinite' " + "value for safety; this may affect results. " + << "Was: " << d << ", is now: " << clamped; + return clamped; +} + // This function will create a new constraint if the indicator constraint has // both a lower bound and an upper bound. absl::Status AddIndicatorConstraint(const MPGeneralConstraintProto& gen_cst, @@ -101,10 +119,11 @@ absl::Status AddIndicatorConstraint(const MPGeneralConstraintProto& gen_cst, } if (ind.constraint().upper_bound() < kInfinity) { + ASSIGN_OR_RETURN(const double upper_bound, + ScipInfClamp(scip, ind.constraint().upper_bound())); RETURN_IF_SCIP_ERROR(SCIPcreateConsIndicator( scip, scip_cst, gen_cst.name().c_str(), ind_var, size, - tmp_variables->data(), tmp_coefficients->data(), - ind.constraint().upper_bound(), + tmp_variables->data(), tmp_coefficients->data(), upper_bound, /*initial=*/!ind.constraint().is_lazy(), /*separate=*/true, /*enforce=*/true, @@ -119,13 +138,14 @@ absl::Status AddIndicatorConstraint(const MPGeneralConstraintProto& gen_cst, scip_cst = &scip_constraints->back(); } if (ind.constraint().lower_bound() > -kInfinity) { + ASSIGN_OR_RETURN(const double lower_bound, + ScipInfClamp(scip, ind.constraint().lower_bound())); for (int i = 0; i < size; ++i) { (*tmp_coefficients)[i] *= -1; } RETURN_IF_SCIP_ERROR(SCIPcreateConsIndicator( scip, scip_cst, gen_cst.name().c_str(), ind_var, size, - tmp_variables->data(), tmp_coefficients->data(), - -ind.constraint().lower_bound(), + tmp_variables->data(), tmp_coefficients->data(), -lower_bound, /*initial=*/!ind.constraint().is_lazy(), /*separate=*/true, /*enforce=*/true, @@ -243,6 +263,10 @@ absl::Status AddQuadraticConstraint( (*tmp_qcoefficients)[i] = quad_cst.qcoefficient(i); } + ASSIGN_OR_RETURN(const double lower_bound, + ScipInfClamp(scip, quad_cst.lower_bound())); + ASSIGN_OR_RETURN(const double upper_bound, + ScipInfClamp(scip, quad_cst.upper_bound())); RETURN_IF_SCIP_ERROR( SCIPcreateConsBasicQuadratic(scip, /*cons=*/scip_cst, @@ -254,8 +278,8 @@ absl::Status AddQuadraticConstraint( /*quadvars1=*/tmp_qvariables1->data(), /*quadvars2=*/tmp_qvariables2->data(), /*quadcoefs=*/tmp_qcoefficients->data(), - /*lhs=*/quad_cst.lower_bound(), - /*rhs=*/quad_cst.upper_bound())); + /*lhs=*/lower_bound, + /*rhs=*/upper_bound)); RETURN_IF_SCIP_ERROR(SCIPaddCons(scip, *scip_cst)); return absl::OkStatus(); } @@ -377,6 +401,8 @@ absl::Status AddMinMaxConstraint(const MPGeneralConstraintProto& gen_cst, CHECK(scip_cst != nullptr); CHECK(tmp_variables != nullptr); CHECK(gen_cst.has_min_constraint() || gen_cst.has_max_constraint()); + constexpr double kInfinity = std::numeric_limits::infinity(); + const auto& minmax = gen_cst.has_min_constraint() ? gen_cst.min_constraint() : gen_cst.max_constraint(); const absl::btree_set unique_var_indices(minmax.var_index().begin(), @@ -393,10 +419,15 @@ absl::Status AddMinMaxConstraint(const MPGeneralConstraintProto& gen_cst, CHECK(vars.size() == vals.size()); const std::string name = gen_cst.has_name() ? absl::StrCat(gen_cst.name(), name_prefix) : ""; + ASSIGN_OR_RETURN(const double scip_lower_bound, + ScipInfClamp(scip, lower_bound)); + ASSIGN_OR_RETURN(const double scip_upper_bound, + ScipInfClamp(scip, upper_bound)); RETURN_IF_SCIP_ERROR(SCIPcreateConsBasicLinear( scip, /*cons=*/&scip_cons, /*name=*/name.c_str(), /*nvars=*/vars.size(), /*vars=*/vars.data(), - /*vals=*/vals.data(), /*lhs=*/lower_bound, /*rhs=*/upper_bound)); + /*vals=*/vals.data(), /*lhs=*/scip_lower_bound, + /*rhs=*/scip_upper_bound)); // Note that the constraints are, by design, not added into the model using // SCIPaddCons. cons.push_back(scip_cons); @@ -427,7 +458,6 @@ absl::Status AddMinMaxConstraint(const MPGeneralConstraintProto& gen_cst, RETURN_IF_SCIP_ERROR(SCIPaddCons(scip, *scip_cst)); // Add all of the inequality constraints. - constexpr double kInfinity = std::numeric_limits::infinity(); cons.clear(); for (const int var_index : unique_var_indices) { vars = {scip_resultant_var, scip_variables[var_index]}; @@ -733,9 +763,13 @@ absl::StatusOr ScipSolveProto( for (int v = 0; v < model.variable_size(); ++v) { const MPVariableProto& variable = model.variable(v); + ASSIGN_OR_RETURN(const double lower_bound, + ScipInfClamp(scip, variable.lower_bound())); + ASSIGN_OR_RETURN(const double upper_bound, + ScipInfClamp(scip, variable.upper_bound())); RETURN_IF_SCIP_ERROR(SCIPcreateVarBasic( scip, /*var=*/&scip_variables[v], /*name=*/variable.name().c_str(), - /*lb=*/variable.lower_bound(), /*ub=*/variable.upper_bound(), + /*lb=*/lower_bound, /*ub=*/upper_bound, /*obj=*/variable.objective_coefficient(), /*vartype=*/variable.is_integer() ? SCIP_VARTYPE_INTEGER : SCIP_VARTYPE_CONTINUOUS)); @@ -754,12 +788,16 @@ absl::StatusOr ScipSolveProto( ct_variables[i] = scip_variables[constraint.var_index(i)]; ct_coefficients[i] = constraint.coefficient(i); } + ASSIGN_OR_RETURN(const double lower_bound, + ScipInfClamp(scip, constraint.lower_bound())); + ASSIGN_OR_RETURN(const double upper_bound, + ScipInfClamp(scip, constraint.upper_bound())); RETURN_IF_SCIP_ERROR(SCIPcreateConsLinear( scip, /*cons=*/&scip_constraints[c], /*name=*/constraint.name().c_str(), /*nvars=*/constraint.var_index_size(), /*vars=*/ct_variables.data(), /*vals=*/ct_coefficients.data(), - /*lhs=*/constraint.lower_bound(), /*rhs=*/constraint.upper_bound(), + /*lhs=*/lower_bound, /*rhs=*/upper_bound, /*initial=*/!constraint.is_lazy(), /*separate=*/true, /*enforce=*/true, diff --git a/ortools/linear_solver/python/model_builder.py b/ortools/linear_solver/python/model_builder.py index e99e36682b..0e5cb1836a 100644 --- a/ortools/linear_solver/python/model_builder.py +++ b/ortools/linear_solver/python/model_builder.py @@ -1427,6 +1427,28 @@ class ModelSolver: values=variables, ) + def reduced_costs(self, variables: _IndexOrSeries) -> pd.Series: + """Returns the reduced cost of the input variables. + + If `variables` is a `pd.Index`, then the output will be indexed by the + variables. If `variables` is a `pd.Series` indexed by the underlying + dimensions, then the output will be indexed by the same underlying + dimensions. + + Args: + variables (Union[pd.Index, pd.Series]): The set of variables from which to + get the values. + + Returns: + pd.Series: The reduced cost of all variables in the set. + """ + if not self.__solve_helper.has_solution(): + return _attribute_series(func=lambda v: pd.NA, values=variables) + return _attribute_series( + func=lambda v: self.__solve_helper.reduced_cost(v.index), + values=variables, + ) + def reduced_cost(self, var: Variable) -> np.double: """Returns the reduced cost of a linear expression after solve.""" if not self.__solve_helper.has_solution(): diff --git a/ortools/linear_solver/python/model_builder_test.py b/ortools/linear_solver/python/model_builder_test.py index be4d989eac..2c6c50fdf6 100644 --- a/ortools/linear_solver/python/model_builder_test.py +++ b/ortools/linear_solver/python/model_builder_test.py @@ -1825,6 +1825,40 @@ class SolverTest(parameterized.TestCase): else: self.assertAlmostEqual(model_solver.objective_value, 0) + def test_simple_problem(self): + # max 5x1 + 4x2 + 3x3 + # s.t 2x1 + 3x2 + x3 <= 5 + # 4x1 + x2 + 2x3 <= 11 + # 3x1 + 4x2 + 2x3 <= 8 + # x1, x2, x3 >= 0 + # Values = (2,0,1) + # Reduced Costs = (0,-3,0) + model = mb.ModelBuilder() + x = model.new_var_series( + "x", pd.Index(range(3)), lower_bounds=0, is_integral=True + ) + self.assertLen(model.get_variables(), 3) + model.maximize(x.dot([5, 4, 3])) + model.add(x.dot([2, 3, 1]) <= 5) + model.add(x.dot([4, 1, 2]) <= 11) + model.add(x.dot([3, 4, 2]) <= 8) + self.assertLen(model.get_linear_constraints(), 3) + solver = mb.ModelSolver("glop") + test_red_cost = solver.reduced_costs(model.get_variables()) + self.assertLen(test_red_cost, 3) + for reduced_cost in test_red_cost: + self.assertTrue(pd.isna(reduced_cost)) + run = solver.solve(model) + self.assertEqual(run, mb.SolveStatus.OPTIMAL) + i = solver.values(model.get_variables()) + self.assertSequenceAlmostEqual(i, [2, 0, 1]) + red_cost = solver.reduced_costs(model.get_variables()) + self.assertSequenceAlmostEqual(red_cost, [0, -3, 0]) + self.assertAlmostEqual(2, solver.value(x[0])) + self.assertAlmostEqual(0, solver.reduced_cost((x[0]))) + self.assertAlmostEqual(-3, solver.reduced_cost((x[1]))) + self.assertAlmostEqual(0, solver.reduced_cost((x[2]))) + if __name__ == "__main__": absltest.main() diff --git a/ortools/linear_solver/scip_interface.cc b/ortools/linear_solver/scip_interface.cc index 958cbcdcbb..0f9e53ae31 100644 --- a/ortools/linear_solver/scip_interface.cc +++ b/ortools/linear_solver/scip_interface.cc @@ -46,6 +46,7 @@ #include "scip/cons_indicator.h" #include "scip/scip.h" #include "scip/scip_copy.h" +#include "scip/scip_numerics.h" #include "scip/scip_param.h" #include "scip/scip_prob.h" #include "scip/scipdefplugins.h" @@ -73,6 +74,8 @@ class SCIPInterface : public MPSolverInterface { const MPModelRequest& request, std::atomic* interrupt) override; void Reset() override; + double infinity() override; + void SetVariableBounds(int var_index, double lb, double ub) override; void SetVariableInteger(int var_index, bool integer) override; void SetConstraintBounds(int row_index, double lb, double ub) override; @@ -308,6 +311,8 @@ absl::Status SCIPInterface::CreateSCIP() { return absl::OkStatus(); } +double SCIPInterface::infinity() { return SCIPinfinity(scip_); } + SCIP* SCIPInterface::DeleteSCIP(bool return_scip) { // NOTE(user): DeleteSCIP() shouldn't "give up" mid-stage if it fails, since // it might be the user's chance to reset the solver to start fresh without diff --git a/ortools/linear_solver/solve.cc b/ortools/linear_solver/solve.cc index 63ef76f5bf..12d949e880 100644 --- a/ortools/linear_solver/solve.cc +++ b/ortools/linear_solver/solve.cc @@ -49,11 +49,10 @@ // --dump_response=/tmp/foo.response \ // 2>/tmp/foo.err -#include #include +#include #include #include -#include #include "absl/flags/flag.h" #include "absl/status/status.h" @@ -61,11 +60,9 @@ #include "absl/strings/match.h" #include "absl/strings/str_format.h" #include "absl/time/time.h" -#include "ortools/base/commandlineflags.h" #include "ortools/base/file.h" #include "ortools/base/helpers.h" #include "ortools/base/init_google.h" -#include "ortools/base/integral_types.h" #include "ortools/base/logging.h" #include "ortools/base/options.h" #include "ortools/linear_solver/linear_solver.h"