diff --git a/pyomo/solvers/plugins/solvers/cuopt_direct.py b/pyomo/solvers/plugins/solvers/cuopt_direct.py index 4ccdc9e631d..bc0849969a4 100644 --- a/pyomo/solvers/plugins/solvers/cuopt_direct.py +++ b/pyomo/solvers/plugins/solvers/cuopt_direct.py @@ -17,6 +17,8 @@ from pyomo.core.base import Suffix, Var, Constraint, Objective from pyomo.core.staleflag import StaleFlagManager from pyomo.repn.linear import LinearRepnVisitor +from pyomo.repn.quadratic import QuadraticRepnVisitor +from pyomo.repn.util import OrderedVarRecorder from pyomo.solvers.plugins.solvers.direct_solver import DirectSolver from pyomo.solvers.plugins.solvers.direct_or_persistent_solver import ( DirectOrPersistentSolver, @@ -35,6 +37,13 @@ def _get_cuopt_version(cuopt, avail): return CUOPTDirect._version = tuple(cuopt.__version__.split('.')) CUOPTDirect._name = "cuOpt %s.%s%s" % CUOPTDirect._version + try: + dm = cuopt.linear_programming.DataModel() + CUOPTDirect._supports_quadratic_constraint = hasattr( + dm, "add_quadratic_constraint" + ) + except Exception: + CUOPTDirect._supports_quadratic_constraint = False cuopt, cuopt_available = attempt_import("cuopt", callback=_get_cuopt_version) @@ -42,6 +51,8 @@ def _get_cuopt_version(cuopt, avail): @SolverFactory.register("cuopt", doc="Direct python interface to CUOPT") class CUOPTDirect(DirectSolver): + _supports_quadratic_constraint = False + def __init__(self, **kwds): kwds["type"] = "cuoptdirect" super().__init__(**kwds) @@ -49,12 +60,20 @@ def __init__(self, **kwds): # Note: Undefined capabilities default to None self._capabilities.linear = True self._capabilities.integer = True + self._capabilities.quadratic_objective = True + self._capabilities.quadratic_constraint = ( + CUOPTDirect._supports_quadratic_constraint + ) self.referenced_vars = ComponentSet() # remove the instance-level definition of the cuopt version: # because the version comes from an imported module, only one # version of cuopt is supported (and stored as a class attribute) del self._version + def _reset_model_flags(self): + self._has_quadratic_content = False + self._has_integer = False + def _apply_solver(self): StaleFlagManager.mark_all_as_stale() log_file = "" @@ -81,8 +100,17 @@ def _add_constraints(self, constraints): matrix_indptr = [0] matrix_indices = [] - # visitor walks expression trees and extracts linear coefficients - visitor = LinearRepnVisitor({}) + if CUOPTDirect._supports_quadratic_constraint: + visitor = QuadraticRepnVisitor( + {}, var_recorder=OrderedVarRecorder({}, {}, None) + ) + var_id_to_ndx = { + id(var): ndx for var, ndx in self._pyomo_var_to_ndx_map.items() + } + else: + visitor = LinearRepnVisitor({}) + var_id_to_ndx = None + con_idx = 0 for con in constraints: if not con.active: @@ -100,6 +128,10 @@ def _add_constraints(self, constraints): "not supported by cuOpt solver." ) + if getattr(repn, 'quadratic', None): + self._add_cuopt_quadratic_constraint(con, repn, visitor, var_id_to_ndx) + continue + # check for trivial constraints after getting repn (more efficient # than walking expression twice with is_fixed) if not repn.linear: @@ -153,6 +185,7 @@ def _add_variables(self, variables): lb, ub = v.bounds if v.is_integer(): v_type.append("I") + self._has_integer = True elif v.is_continuous(): v_type.append("C") else: @@ -169,8 +202,110 @@ def _add_variables(self, variables): self._solver_model.set_variable_types(np.array(v_type)) self._solver_model.set_variable_names(np.array(v_names)) + @staticmethod + def _linear_repn_to_coo(linear, var_id_to_ndx): + if not linear: + return None, None + indices = sorted(linear) + values = np.array([linear[i] for i in indices], dtype=np.float64) + return values, np.array([var_id_to_ndx[i] for i in indices], dtype=np.int32) + + @staticmethod + def _quadratic_repn_to_coo(quadratic, var_id_to_ndx): + nnz = len(quadratic) + vals = np.empty(nnz, dtype=np.float64) + rows = np.empty(nnz, dtype=np.int32) + cols = np.empty(nnz, dtype=np.int32) + for i, ((var_id1, var_id2), coef) in enumerate(quadratic.items()): + rows[i] = var_id_to_ndx[var_id1] + cols[i] = var_id_to_ndx[var_id2] + vals[i] = coef + return vals, rows, cols + + def _add_cuopt_quadratic_constraint(self, con, qrepn, visitor, var_id_to_ndx): + from pyomo.core.expr.numvalue import is_fixed, value + + if con.equality: + raise ValueError( + f"Equality is not supported for quadratic constraint '{con.name}' " + "in cuOpt." + ) + if con.has_lb() and con.has_ub(): + raise ValueError( + f"Ranged quadratic constraints are not supported by cuOpt " + f"('{con.name}')." + ) + if con.has_lb(): + if not is_fixed(con.lower): + raise ValueError( + f"Lower bound of constraint '{con.name}' is not constant." + ) + sense = "G" + rhs = value(con.lower) - qrepn.constant + elif con.has_ub(): + if not is_fixed(con.upper): + raise ValueError( + f"Upper bound of constraint '{con.name}' is not constant." + ) + sense = "L" + rhs = value(con.upper) - qrepn.constant + else: + raise ValueError( + f"Constraint '{con.name}' does not have a lower or an upper bound." + ) + + linear_values, linear_indices = self._linear_repn_to_coo( + qrepn.linear, var_id_to_ndx + ) + vals, rows, cols = self._quadratic_repn_to_coo(qrepn.quadratic, var_id_to_ndx) + con_name = self._symbol_map.getSymbol(con, self._labeler) + self._solver_model.add_quadratic_constraint( + constraint_row_name=con_name, + linear_values=linear_values, + linear_indices=linear_indices, + rhs_value=rhs, + vals=vals, + rows=rows, + cols=cols, + sense=sense, + ) + self._has_quadratic_content = True + for var_id, coef in qrepn.linear.items(): + if coef: + self.referenced_vars.add(visitor.var_map[var_id]) + for var_id1, var_id2 in qrepn.quadratic: + self.referenced_vars.add(visitor.var_map[var_id1]) + self.referenced_vars.add(visitor.var_map[var_id2]) + + @staticmethod + def _build_quadratic_objective_csr(quadratic, var_id_to_ndx, num_vars): + nnz = len(quadratic) + if nnz == 0: + return None, None, None + + data = np.empty(nnz, dtype=np.float64) + rows = np.empty(nnz, dtype=np.int32) + cols = np.empty(nnz, dtype=np.int32) + for i, ((var_id1, var_id2), coef) in enumerate(quadratic.items()): + rows[i] = var_id_to_ndx[var_id1] + cols[i] = var_id_to_ndx[var_id2] + data[i] = coef + + order = np.lexsort((cols, rows)) + rows = rows[order] + cols = cols[order] + data = data[order] + + indptr = np.empty(num_vars + 1, dtype=np.int32) + indptr[0] = 0 + np.cumsum(np.bincount(rows, minlength=num_vars), out=indptr[1:]) + + return data, cols, indptr + def _set_objective(self, objective): - visitor = LinearRepnVisitor({}) + visitor = QuadraticRepnVisitor( + {}, var_recorder=OrderedVarRecorder({}, {}, None) + ) repn = visitor.walk_expression(objective.expr) if repn.nonlinear is not None: raise ValueError( @@ -178,19 +313,38 @@ def _set_objective(self, objective): "not supported by cuOpt solver." ) - obj_coeffs = [0] * len(self._pyomo_var_to_ndx_map) - # repn.linear is keyed by id(var), use var_map to get actual vars + num_vars = len(self._pyomo_var_to_ndx_map) + obj_coeffs = [0] * num_vars for var_id, coef in repn.linear.items(): var = visitor.var_map[var_id] obj_coeffs[self._pyomo_var_to_ndx_map[var]] = coef self.referenced_vars.add(var) + self._solver_model.set_objective_coefficients(np.array(obj_coeffs)) + if repn.constant: + self._solver_model.set_objective_offset(repn.constant) self._solver_model.set_maximize(objective.sense == maximize) + if repn.quadratic: + var_id_to_ndx = { + id(var): ndx for var, ndx in self._pyomo_var_to_ndx_map.items() + } + q_data, q_indices, q_indptr = self._build_quadratic_objective_csr( + repn.quadratic, var_id_to_ndx, num_vars + ) + self._solver_model.set_quadratic_objective_matrix( + q_data, q_indices, q_indptr + ) + self._has_quadratic_content = True + for var_id1, var_id2 in repn.quadratic: + self.referenced_vars.add(visitor.var_map[var_id1]) + self.referenced_vars.add(visitor.var_map[var_id2]) + def _set_instance(self, model, kwds={}): DirectOrPersistentSolver._set_instance(self, model, kwds) self._pyomo_var_to_ndx_map = ComponentMap() self._ndx_count = 0 + self._reset_model_flags() try: self._solver_model = cuopt.linear_programming.DataModel() @@ -220,6 +374,11 @@ def _add_block(self, block): raise ValueError("Solver interface does not support multiple objectives.") elif objectives: self._set_objective(objectives[0]) + if self._has_quadratic_content and self._has_integer: + raise ValueError( + "cuOpt does not support mixed-integer problems with quadratic " + "or second-order cone constraints." + ) def _postsolve(self): extract_duals = False @@ -337,10 +496,21 @@ def _postsolve(self): if extract_duals: logger.warning("Cannot get duals for MIP.") else: - if extract_reduced_costs: - reduced_costs = solution.get_reduced_cost() - if extract_duals: - dual_solution = solution.get_dual_solution() + if self._has_quadratic_content: + if extract_reduced_costs: + logger.warning( + "Cannot get reduced costs for quadratic or conic " + "problems in cuOpt." + ) + if extract_duals: + logger.warning( + "Cannot get duals for quadratic or conic problems in cuOpt." + ) + else: + if extract_reduced_costs: + reduced_costs = solution.get_reduced_cost() + if extract_duals: + dual_solution = solution.get_dual_solution() if self._save_results: soln_variables = soln.variable @@ -409,6 +579,10 @@ def load_rc(self, vars_to_load=None): is_mip = self.solution.get_problem_category() if is_mip: logger.warning("Cannot get reduced costs for MIP.") + elif self._has_quadratic_content: + logger.warning( + "Cannot get reduced costs for quadratic or conic problems in cuOpt." + ) else: self._load_rc(vars_to_load) @@ -434,5 +608,7 @@ def load_duals(self, cons_to_load=None): is_mip = self.solution.get_problem_category() if is_mip: logger.warning("Cannot get duals for MIP.") + elif self._has_quadratic_content: + logger.warning("Cannot get duals for quadratic or conic problems in cuOpt.") else: self._load_duals(cons_to_load) diff --git a/pyomo/solvers/tests/checks/test_cuopt_direct.py b/pyomo/solvers/tests/checks/test_cuopt_direct.py index fa6620a603e..30bb57b3d8e 100644 --- a/pyomo/solvers/tests/checks/test_cuopt_direct.py +++ b/pyomo/solvers/tests/checks/test_cuopt_direct.py @@ -7,7 +7,6 @@ # software. This software is distributed under the 3-clause BSD License. # ____________________________________________________________________________________ -import os from pyomo.environ import ( SolverFactory, ConcreteModel, @@ -24,14 +23,15 @@ Binary, ) import pytest -from pyomo.common.dependencies import attempt_import -from pyomo.opt import check_available_solvers -from pyomo.common.tee import capture_output -from pyomo.common.tempfiles import TempfileManager +from pyomo.opt import TerminationCondition import pyomo.common.unittest as unittest from pyomo.solvers.plugins.solvers.cuopt_direct import cuopt_available, CUOPTDirect +def _cuopt_supports_qc(): + return cuopt_available and CUOPTDirect._supports_quadratic_constraint + + def _cuopt_at_least(*required): """True iff cuOpt is available and at least the given (major, minor[, patch]) version.""" if not cuopt_available: @@ -162,13 +162,90 @@ def test_nonlinear_constraint_rejected(self): m.x = Var(domain=NonNegativeReals) m.y = Var(domain=NonNegativeReals) m.obj = Objective(expr=m.x + m.y, sense=minimize) - # nonlinear constraint: x * y <= 10 - m.nonlin_con = Constraint(expr=m.x * m.y <= 10) + # nonlinear constraint (degree > 2) + m.nonlin_con = Constraint(expr=m.x**3 <= 10) opt = SolverFactory('cuopt') with pytest.raises(ValueError, match=r"contains nonlinear terms"): opt.solve(m) + @unittest.skipIf(not cuopt_available, "The CuOpt solver is not available") + def test_quadratic_objective_matrix(self): + m = ConcreteModel() + m.x = Var(domain=NonNegativeReals) + m.y = Var(domain=NonNegativeReals) + m.obj = Objective(expr=m.x**2 + 2 * m.x + 4 * m.y**2 + 3.0, sense=minimize) + m.c = Constraint(expr=m.x + m.y >= 1) + + opt = SolverFactory('cuopt') + opt._set_instance(m) + + q_values = opt._solver_model.get_quadratic_objective_values() + q_indices = opt._solver_model.get_quadratic_objective_indices() + q_offsets = opt._solver_model.get_quadratic_objective_offsets() + self.assertEqual(list(q_values), [1.0, 4.0]) + self.assertEqual(list(q_indices), [0, 1]) + self.assertEqual(list(q_offsets), [0, 1, 2]) + self.assertAlmostEqual(opt._solver_model.get_objective_offset(), 3.0) + c = opt._solver_model.get_objective_coefficients() + self.assertAlmostEqual(c[0], 2.0) + self.assertAlmostEqual(c[1], 0.0) + + @unittest.skipUnless( + _cuopt_supports_qc(), + "cuOpt quadratic constraint support requires add_quadratic_constraint", + ) + def test_quadratic_constraint_soc(self): + m = ConcreteModel() + m.x0 = Var(bounds=(None, None)) + m.x1 = Var(bounds=(None, None)) + m.x2 = Var(bounds=(None, None)) + m.y = Var(bounds=(0, 5)) + m.obj = Objective(expr=3 * m.x0 + 2 * m.x1 + m.x2, sense=minimize) + m.soc = Constraint( + expr=m.x0 * m.x0 + m.x1 * m.x1 + m.x2 * m.x2 - m.y * m.y <= 0 + ) + m.lin = Constraint(expr=m.x0 + m.x1 + 3 * m.x2 >= 1) + + opt = SolverFactory('cuopt') + opt._set_instance(m) + + qcs = opt._solver_model.get_quadratic_constraints() + self.assertEqual(len(qcs), 1) + qc = qcs[0] + self.assertEqual(qc["constraint_row_type"], "L") + self.assertAlmostEqual(qc["rhs_value"], 0.0) + self.assertEqual(list(qc["vals"]), [1.0, 1.0, 1.0, -1.0]) + self.assertEqual(list(qc["rows"]), [0, 1, 2, 3]) + self.assertEqual(list(qc["cols"]), [0, 1, 2, 3]) + self.assertTrue(opt._has_quadratic_content) + + res = opt.solve(m) + self.assertEqual(res.solver.termination_condition, TerminationCondition.optimal) + self.assertAlmostEqual(value(m.obj), -13.548638872057857, places=4) + self.assertAlmostEqual(value(m.x0), -3.874621914903146, places=4) + self.assertAlmostEqual(value(m.x1), -2.129788270186565, places=4) + self.assertAlmostEqual(value(m.x2), 2.3348034130247104, places=4) + self.assertAlmostEqual(value(m.y), 5.0, places=4) + + @unittest.skipUnless( + _cuopt_supports_qc(), + "cuOpt quadratic constraint support requires add_quadratic_constraint", + ) + def test_mip_with_quadratic_constraint_rejected(self): + m = ConcreteModel() + m.x = Var(within=Binary) + m.y = Var(bounds=(0, None)) + m.z = Var() + m.obj = Objective(expr=m.x + m.z, sense=minimize) + m.soc = Constraint(expr=m.z * m.z - m.y * m.y <= 0) + + opt = SolverFactory('cuopt') + with pytest.raises( + ValueError, match=r"does not support mixed-integer problems with quadratic" + ): + opt.solve(m) + if __name__ == "__main__": unittest.main() diff --git a/pyomo/solvers/tests/mip/test_qp.py b/pyomo/solvers/tests/mip/test_qp.py index 98c3243ac6e..3ce902d68f7 100644 --- a/pyomo/solvers/tests/mip/test_qp.py +++ b/pyomo/solvers/tests/mip/test_qp.py @@ -30,6 +30,8 @@ xpress_persistent = SolverFactory('xpress_persistent') xpress_appsi = SolverFactory('appsi_xpress') +cuopt = SolverFactory('cuopt') + class TestQuadraticModels(unittest.TestCase): def _qp_model(self): @@ -192,3 +194,25 @@ def test_qp_objective_xpress_appsi(self): xpress_appsi.set_instance(m) results = xpress_appsi.solve(m) self.assertEqual(m.obj(), results['Problem'][0]['Upper bound']) + + @unittest.skipUnless( + cuopt.available(exception_flag=False), "needs cuOpt direct interface" + ) + def test_qp_objective_cuopt_model(self): + m = self._qp_model() + cuopt._set_instance(m) + q_values = cuopt._solver_model.get_quadratic_objective_values() + q_indices = cuopt._solver_model.get_quadratic_objective_indices() + q_offsets = cuopt._solver_model.get_quadratic_objective_offsets() + self.assertEqual(list(q_values), [10000.0, 10000.0, 1000.0, 100000.0]) + self.assertEqual(list(q_indices), [0, 1, 2, 2]) + self.assertEqual(list(q_offsets), [0, 1, 3, 4]) + + @unittest.skipUnless( + cuopt.available(exception_flag=False), "needs cuOpt direct interface" + ) + def test_qp_objective_cuopt(self): + m = self._qp_model() + results = cuopt.solve(m) + # cuOpt's barrier method may return slightly inexact primals + self.assertAlmostEqual(m.obj(), results['Problem'][0]['Upper bound'], places=6) diff --git a/pyomo/solvers/tests/solvers.py b/pyomo/solvers/tests/solvers.py index 98c204d1851..beffbf648a0 100644 --- a/pyomo/solvers/tests/solvers.py +++ b/pyomo/solvers/tests/solvers.py @@ -440,7 +440,14 @@ def test_solver_cases(*args): # # CUOPT # - _cuopt_capabilities = set(['linear', 'integer']) + _cuopt_capabilities = set(['linear', 'integer', 'quadratic_objective']) + try: + from pyomo.solvers.plugins.solvers.cuopt_direct import CUOPTDirect + + if CUOPTDirect._supports_quadratic_constraint: + _cuopt_capabilities.add('quadratic_constraint') + except ImportError: + pass _test_solver_cases['cuopt', 'python'] = initialize( name='cuopt',