diff --git a/.github/workflows/test_branches.yml b/.github/workflows/test_branches.yml index f96aef36e1f..a902a74e15a 100644 --- a/.github/workflows/test_branches.yml +++ b/.github/workflows/test_branches.yml @@ -334,6 +334,12 @@ jobs: || echo "WARNING: Xpress Community Edition is not available" python -m pip install --cache-dir cache/pip maingopy \ || echo "WARNING: MAiNGO is not available" + if [[ ${{matrix.python}} == pypy* ]]; then + echo "skipping SCIP for pypy" + else + python -m pip install --cache-dir cache/pip pyscipopt \ + || echo "WARNING: SCIP is not available" + fi if [[ ${{matrix.python}} == pypy* ]]; then echo "skipping wntr for pypy" else @@ -423,7 +429,7 @@ jobs: else XPRESS='xpress' fi - for PKG in "$CPLEX" docplex gurobi "$XPRESS" cyipopt pymumps scip; do + for PKG in "$CPLEX" docplex gurobi "$XPRESS" cyipopt pymumps scip pyscipopt; do echo "" echo "*** Install $PKG ***" echo "" diff --git a/.github/workflows/test_pr_and_main.yml b/.github/workflows/test_pr_and_main.yml index fc5bf334435..1f2f648d47c 100644 --- a/.github/workflows/test_pr_and_main.yml +++ b/.github/workflows/test_pr_and_main.yml @@ -381,6 +381,12 @@ jobs: || echo "WARNING: Xpress Community Edition is not available" python -m pip install --cache-dir cache/pip maingopy \ || echo "WARNING: MAiNGO is not available" + if [[ ${{matrix.python}} == pypy* ]]; then + echo "skipping SCIP for pypy" + else + python -m pip install --cache-dir cache/pip pyscipopt \ + || echo "WARNING: SCIP is not available" + fi if [[ ${{matrix.python}} == pypy* ]]; then echo "skipping wntr for pypy" else @@ -470,7 +476,7 @@ jobs: else XPRESS='xpress' fi - for PKG in "$CPLEX" docplex gurobi "$XPRESS" cyipopt pymumps scip; do + for PKG in "$CPLEX" docplex gurobi "$XPRESS" cyipopt pymumps scip pyscipopt; do echo "" echo "*** Install $PKG ***" echo "" diff --git a/doc/OnlineDocs/explanation/experimental/solvers.rst b/doc/OnlineDocs/explanation/experimental/solvers.rst index c04a7d870b8..5701e37d188 100644 --- a/doc/OnlineDocs/explanation/experimental/solvers.rst +++ b/doc/OnlineDocs/explanation/experimental/solvers.rst @@ -45,6 +45,9 @@ with existing interfaces). * - Ipopt - ``ipopt`` - ``ipopt_v2`` + * - GAMS + - ``gams`` + - ``gams_v2`` * - Gurobi (persistent) - ``gurobi_persistent`` - ``gurobi_persistent_v2`` @@ -57,9 +60,12 @@ with existing interfaces). * - KNITRO - ``knitro_direct`` - ``knitro_direct`` - * - GAMS - - ``gams`` - - ``gams_v2`` + * - SCIP (direct) + - ``scip_direct`` + - ``scip_direct`` + * - SCIP (persistent) + - ``scip_persistent`` + - ``scip_persistent`` Using the new interfaces through the legacy interface ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^ diff --git a/doc/OnlineDocs/getting_started/solvers.rst b/doc/OnlineDocs/getting_started/solvers.rst index ba04217b8f0..8183905c093 100644 --- a/doc/OnlineDocs/getting_started/solvers.rst +++ b/doc/OnlineDocs/getting_started/solvers.rst @@ -76,11 +76,16 @@ the license requirements for their desired solver. - ``conda install ‑c conda‑forge pymumps`` - `License `__ `Docs `__ - * - SCIP + * - SCIP (Command-line) - N/A - ``conda install ‑c conda‑forge scip`` - `License `__ `Docs `__ + * - SCIP (Python) + - ``pip install pyscipopt`` + - ``conda install ‑c conda‑forge pyscipopt`` + - `License `__ + `Docs `__ * - XPRESS - ``pip install xpress`` - ``conda install ‑c fico‑xpress xpress`` diff --git a/pyomo/contrib/solver/common/solution_loader.py b/pyomo/contrib/solver/common/solution_loader.py index 8223f33760a..faaf7c75685 100644 --- a/pyomo/contrib/solver/common/solution_loader.py +++ b/pyomo/contrib/solver/common/solution_loader.py @@ -20,6 +20,9 @@ from pyomo.core.staleflag import StaleFlagManager from pyomo.core.base.suffix import Suffix from .util import NoSolutionError +import logging + +logger = logging.getLogger(__name__) class SolutionLoader: diff --git a/pyomo/contrib/solver/plugins.py b/pyomo/contrib/solver/plugins.py index c16c369abd3..4fa18dc9694 100644 --- a/pyomo/contrib/solver/plugins.py +++ b/pyomo/contrib/solver/plugins.py @@ -14,6 +14,8 @@ from .solvers.gurobi.gurobi_persistent import GurobiPersistent from .solvers.gurobi.gurobi_direct_minlp import GurobiDirectMINLP from .solvers.highs import Highs +from .solvers.scip.scip_direct import ScipDirect +from .solvers.scip.scip_persistent import ScipPersistent from .solvers.gams import GAMS from .solvers.knitro.direct import KnitroDirectSolver @@ -22,6 +24,9 @@ def load(): SolverFactory.register( name="ipopt", legacy_name="ipopt_v2", doc="The IPOPT NLP solver" )(Ipopt, LegacyIpoptSolver) + SolverFactory.register(name='gams', legacy_name='gams_v2', doc='Interface to GAMS')( + GAMS + ) SolverFactory.register( name="gurobi_persistent", legacy_name="gurobi_persistent_v2", @@ -40,11 +45,16 @@ def load(): SolverFactory.register( name="highs", legacy_name="highs", doc="Persistent interface to HiGHS" )(Highs) - SolverFactory.register(name='gams', legacy_name='gams_v2', doc='Interface to GAMS')( - GAMS - ) SolverFactory.register( name="knitro_direct", legacy_name="knitro_direct", doc="Direct interface to KNITRO solver", )(KnitroDirectSolver) + SolverFactory.register( + name='scip_direct', legacy_name='scip_direct', doc='Direct interface pyscipopt' + )(ScipDirect) + SolverFactory.register( + name='scip_persistent', + legacy_name='scip_persistent', + doc='Persistent interface pyscipopt', + )(ScipPersistent) diff --git a/pyomo/contrib/solver/solvers/scip/__init__.py b/pyomo/contrib/solver/solvers/scip/__init__.py new file mode 100644 index 00000000000..231b44987f6 --- /dev/null +++ b/pyomo/contrib/solver/solvers/scip/__init__.py @@ -0,0 +1,8 @@ +# ____________________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2026 National Technology and Engineering Solutions of Sandia, LLC +# Under the terms of Contract DE-NA0003525 with National Technology and Engineering +# Solutions of Sandia, LLC, the U.S. Government retains certain rights in this +# software. This software is distributed under the 3-clause BSD License. +# ____________________________________________________________________________________ diff --git a/pyomo/contrib/solver/solvers/scip/base.py b/pyomo/contrib/solver/solvers/scip/base.py new file mode 100644 index 00000000000..73bf69e1528 --- /dev/null +++ b/pyomo/contrib/solver/solvers/scip/base.py @@ -0,0 +1,323 @@ +# ____________________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2026 National Technology and Engineering Solutions of Sandia, LLC +# Under the terms of Contract DE-NA0003525 with National Technology and Engineering +# Solutions of Sandia, LLC, the U.S. Government retains certain rights in this +# software. This software is distributed under the 3-clause BSD License. +# ____________________________________________________________________________________ + +from __future__ import annotations + +import logging +import math +from typing import List, Optional, Sequence, Mapping + +from pyomo.common.collections import ComponentMap +from pyomo.common.numeric_types import native_numeric_types +from pyomo.common.config import ConfigValue +from pyomo.common.dependencies import attempt_import +from pyomo.contrib.fbbt.fbbt import compute_bounds_on_expr +from pyomo.contrib.solver.common.config import BranchAndBoundConfig +from pyomo.contrib.solver.common.solution_loader import SolutionLoader +from pyomo.contrib.solver.common.util import NoSolutionError +from pyomo.core.base.expression import ExpressionData, ScalarExpression +from pyomo.core.base.param import ParamData, ScalarParam +from pyomo.core.base.units_container import _PyomoUnit +from pyomo.core.base.var import VarData, ScalarVar +from pyomo.core.expr.numvalue import is_constant, NumericConstant +from pyomo.core.expr.numeric_expr import ( + NegationExpression, + PowExpression, + ProductExpression, + MonomialTermExpression, + DivisionExpression, + SumExpression, + LinearExpression, + UnaryFunctionExpression, + NPV_NegationExpression, + NPV_PowExpression, + NPV_ProductExpression, + NPV_DivisionExpression, + NPV_SumExpression, + NPV_UnaryFunctionExpression, +) +from pyomo.core.expr.relational_expr import ( + EqualityExpression, + InequalityExpression, + RangedExpression, +) +from pyomo.core.expr.visitor import StreamBasedExpressionVisitor +from pyomo.gdp.disjunct import AutoLinkedBinaryVar + +logger = logging.getLogger(__name__) + +scip, scip_available = attempt_import('pyscipopt') + + +class ScipConfig(BranchAndBoundConfig): + def __init__( + self, + description=None, + doc=None, + implicit=False, + implicit_domain=None, + visibility=0, + ): + BranchAndBoundConfig.__init__( + self, + description=description, + doc=doc, + implicit=implicit, + implicit_domain=implicit_domain, + visibility=visibility, + ) + self.warmstart_discrete_vars: bool = self.declare( + 'warmstart_discrete_vars', + ConfigValue( + default=False, + domain=bool, + description="If True, the current values of the integer variables " + "will be passed to Scip.", + ), + ) + + +def _handle_var(node, data, opt, visitor): + if node not in opt._pyomo_var_to_solver_var_map: + scip_var = opt._add_var(node) + else: + scip_var = opt._pyomo_var_to_solver_var_map[node] + return scip_var + + +def _handle_param(node, data, opt, visitor): + # for the persistent interface, we create scip variables in place + # of parameters. However, this makes things complicated for range + # constraints because scip does not allow variables in the + # lower and upper parts of range constraints + if visitor.in_range: + return node.value + if not opt.is_persistent(): + return node.value + if node.is_constant(): + return node.value + if node not in opt._pyomo_param_to_solver_param_map: + scip_param = opt._add_param(node) + else: + scip_param = opt._pyomo_param_to_solver_param_map[node] + return scip_param + + +def _handle_constant(node, data, opt, visitor): + return node.value + + +def _handle_float(node, data, opt, visitor): + return float(node) + + +def _handle_negation(node, data, opt, visitor): + return -data[0] + + +def _handle_pow(node, data, opt, visitor): + x, y = data # x ** y = exp(log(x**y)) = exp(y*log(x)) + if is_constant(node.args[1]): + return x**y + else: + xlb, xub = compute_bounds_on_expr(node.args[0]) + if xlb > 0: + return scip.exp(y * scip.log(x)) + else: + return x**y # scip will probably raise an error here + + +def _handle_product(node, data, opt, visitor): + assert len(data) == 2 + return data[0] * data[1] + + +def _handle_division(node, data, opt, visitor): + return data[0] / data[1] + + +def _handle_sum(node, data, opt, visitor): + return sum(data) + + +def _handle_exp(node, data, opt, visitor): + return scip.exp(data[0]) + + +def _handle_log(node, data, opt, visitor): + return scip.log(data[0]) + + +def _handle_log10(node, data, opt, visitor): + return scip.log(data[0]) / math.log(10) + + +def _handle_sin(node, data, opt, visitor): + return scip.sin(data[0]) + + +def _handle_cos(node, data, opt, visitor): + return scip.cos(data[0]) + + +def _handle_sqrt(node, data, opt, visitor): + return scip.sqrt(data[0]) + + +def _handle_abs(node, data, opt, visitor): + return abs(data[0]) + + +def _handle_tan(node, data, opt, visitor): + return scip.sin(data[0]) / scip.cos(data[0]) + + +def _handle_tanh(node, data, opt, visitor): + x = data[0] + _exp = scip.exp + return (_exp(x) - _exp(-x)) / (_exp(x) + _exp(-x)) + + +_unary_map = { + 'exp': _handle_exp, + 'log': _handle_log, + 'sin': _handle_sin, + 'cos': _handle_cos, + 'sqrt': _handle_sqrt, + 'abs': _handle_abs, + 'tan': _handle_tan, + 'log10': _handle_log10, + 'tanh': _handle_tanh, +} + + +def _handle_unary(node, data, opt, visitor): + if node.getname() in _unary_map: + return _unary_map[node.getname()](node, data, opt, visitor) + else: + raise NotImplementedError(f'unable to handle unary expression: {str(node)}') + + +def _handle_equality(node, data, opt, visitor): + return data[0] == data[1] + + +def _handle_ranged(node, data, opt, visitor): + # note that the lower and upper parts of the + # range constraint cannot have variables + return data[0] <= (data[1] <= data[2]) + + +def _handle_inequality(node, data, opt, visitor): + return data[0] <= data[1] + + +def _handle_named_expression(node, data, opt, visitor): + return data[0] + + +def _handle_unit(node, data, opt, visitor): + return node.value + + +_operator_map = { + NegationExpression: _handle_negation, + PowExpression: _handle_pow, + ProductExpression: _handle_product, + MonomialTermExpression: _handle_product, + DivisionExpression: _handle_division, + SumExpression: _handle_sum, + LinearExpression: _handle_sum, + UnaryFunctionExpression: _handle_unary, + NPV_NegationExpression: _handle_negation, + NPV_PowExpression: _handle_pow, + NPV_ProductExpression: _handle_product, + NPV_DivisionExpression: _handle_division, + NPV_SumExpression: _handle_sum, + NPV_UnaryFunctionExpression: _handle_unary, + EqualityExpression: _handle_equality, + RangedExpression: _handle_ranged, + InequalityExpression: _handle_inequality, + ScalarExpression: _handle_named_expression, + ExpressionData: _handle_named_expression, + VarData: _handle_var, + ScalarVar: _handle_var, + ParamData: _handle_param, + ScalarParam: _handle_param, + float: _handle_float, + int: _handle_float, + AutoLinkedBinaryVar: _handle_var, + _PyomoUnit: _handle_unit, + NumericConstant: _handle_constant, +} + + +class _PyomoToScipVisitor(StreamBasedExpressionVisitor): + def __init__(self, solver, **kwds): + super().__init__(**kwds) + self.solver = solver + self.in_range = False + + def initializeWalker(self, expr): + self.in_range = False + return True, None + + def exitNode(self, node, data): + nt = type(node) + if nt in _operator_map: + return _operator_map[nt](node, data, self.solver, self) + elif nt in native_numeric_types: + _operator_map[nt] = _handle_float + return _handle_float(node, data, self.solver, self) + else: + raise NotImplementedError(f'unrecognized expression type: {nt}') + + def enterNode(self, node): + if type(node) is RangedExpression: + self.in_range = True + return None, [] + + +class ScipSolutionLoader(SolutionLoader): + def __init__(self, solver_model, var_map, con_map, pyomo_model, opt) -> None: + super().__init__() + self._solver_model = solver_model + self._var_map = var_map + self._con_map = con_map + self._pyomo_model = pyomo_model + # make sure the scip model does not get freed until the solution loader is garbage collected + self._opt = opt + self._active_solution_id = 0 + + def _set_solution_id(self, solution_id: int) -> int: + if solution_id is None: + solution_id = 0 + previous_id = self._active_solution_id + self._active_solution_id = solution_id + return previous_id + + def get_number_of_solutions(self) -> int: + return self._solver_model.getNSols() + + def get_solution_ids(self) -> List: + return list(range(self.get_number_of_solutions())) + + def get_vars( + self, vars_to_load: Optional[Sequence[VarData]] = None + ) -> Mapping[VarData, float]: + if self.get_number_of_solutions() == 0: + raise NoSolutionError() + if vars_to_load is None: + vars_to_load = list(self._var_map.keys()) + sol = self._solver_model.getSols()[self._active_solution_id] + res = ComponentMap() + for v in vars_to_load: + sv = self._var_map[v] + res[v] = sol[sv] + return res diff --git a/pyomo/contrib/solver/solvers/scip/scip_direct.py b/pyomo/contrib/solver/solvers/scip/scip_direct.py new file mode 100644 index 00000000000..58281218b8e --- /dev/null +++ b/pyomo/contrib/solver/solvers/scip/scip_direct.py @@ -0,0 +1,404 @@ +# ____________________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2026 National Technology and Engineering Solutions of Sandia, LLC +# Under the terms of Contract DE-NA0003525 with National Technology and Engineering +# Solutions of Sandia, LLC, the U.S. Government retains certain rights in this +# software. This software is distributed under the 3-clause BSD License. +# ____________________________________________________________________________________ + +from __future__ import annotations + +import datetime +import io +import logging +import math +from typing import Tuple, List + +from pyomo.common.collections import ComponentMap +from pyomo.common.errors import InfeasibleConstraintException +from pyomo.common.timing import HierarchicalTimer +from pyomo.common.tee import capture_output, TeeStream + +from pyomo.core.base.block import BlockData +from pyomo.core.base.constraint import Constraint, ConstraintData +from pyomo.core.base.sos import SOSConstraint, SOSConstraintData +from pyomo.core.kernel.objective import minimize, maximize +from pyomo.core.staleflag import StaleFlagManager + +from pyomo.contrib.solver.common.base import SolverBase, Availability +from pyomo.contrib.solver.common.solution_loader import NoSolutionSolutionLoader +from pyomo.contrib.solver.common.results import ( + Results, + SolutionStatus, + TerminationCondition, +) +from pyomo.contrib.solver.common.util import ( + NoFeasibleSolutionError, + NoOptimalSolutionError, + get_objective, +) + +from pyomo.contrib.solver.solvers.scip.base import ( + scip, + scip_available, + ScipConfig, + _PyomoToScipVisitor, + ScipSolutionLoader, +) + +logger = logging.getLogger(__name__) + + +class ScipDirect(SolverBase): + + _available = None + _tc_map = None + _minimum_version = (5, 5, 0) # this is probably conservative + + CONFIG = ScipConfig() + + def __init__(self, **kwds): + super().__init__(**kwds) + self._solver_model = None + self._pyomo_var_to_solver_var_map = ComponentMap() + self._pyomo_con_to_solver_con_map = {} + self._pyomo_param_to_solver_param_map = ( + ComponentMap() + ) # param to scip var with equal bounds + self._pyomo_sos_to_solver_sos_map = {} + self._expr_visitor = _PyomoToScipVisitor(self) + self._objective = None # pyomo objective + self._obj_var = ( + None # a scip variable because the objective cannot be nonlinear + ) + self._obj_con = None # a scip constraint (obj_var >= obj_expr) + + def _clear(self): + self._solver_model = None + self._pyomo_var_to_solver_var_map = ComponentMap() + self._pyomo_con_to_solver_con_map = {} + self._pyomo_param_to_solver_param_map = ComponentMap() + self._pyomo_sos_to_solver_sos_map = {} + self._objective = None + self._obj_var = None + self._obj_con = None + + def available(self) -> Availability: + if self._available is not None: + return self._available + + if not scip_available: + ScipDirect._available = Availability.NotFound + elif self.version() < self._minimum_version: + ScipDirect._available = Availability.BadVersion + else: + ScipDirect._available = Availability.FullLicense + + return self._available + + def version(self) -> Tuple: + return tuple(int(i) for i in scip.__version__.split('.')) + + def solve(self, model: BlockData, **kwds) -> Results: + start_timestamp = datetime.datetime.now(datetime.timezone.utc) + config = self.config(value=kwds, preserve_implicit=True) + + StaleFlagManager.mark_all_as_stale() + + if config.timer is None: + config.timer = HierarchicalTimer() + timer = config.timer + + ostreams = [io.StringIO()] + config.tee + + scip_model, solution_loader, has_obj = self._create_solver_model(model, config) + + scip_model.hideOutput(quiet=False) + if config.threads is not None: + scip_model.setParam('lp/threads', config.threads) + if config.time_limit is not None: + scip_model.setParam('limits/time', config.time_limit) + if config.rel_gap is not None: + scip_model.setParam('limits/gap', config.rel_gap) + if config.abs_gap is not None: + scip_model.setParam('limits/absgap', config.abs_gap) + + if config.warmstart_discrete_vars: + self._mipstart() + + for key, option in config.solver_options.items(): + scip_model.setParam(key, option) + + timer.start('optimize') + with capture_output(TeeStream(*ostreams), capture_fd=True): + scip_model.optimize() + timer.stop('optimize') + + results = self._populate_results(scip_model, solution_loader, has_obj, config) + + results.solver_log = ostreams[0].getvalue() + end_timestamp = datetime.datetime.now(datetime.timezone.utc) + results.timing_info.start_timestamp = start_timestamp + results.timing_info.wall_time = ( + end_timestamp - start_timestamp + ).total_seconds() + results.timing_info.timer = timer + return results + + def _get_tc_map(self): + if ScipDirect._tc_map is None: + tc = TerminationCondition + ScipDirect._tc_map = { + "unknown": tc.unknown, + "userinterrupt": tc.interrupted, + "nodelimit": tc.iterationLimit, + "totalnodelimit": tc.iterationLimit, + "stallnodelimit": tc.iterationLimit, + "timelimit": tc.maxTimeLimit, + "memlimit": tc.unknown, + "gaplimit": tc.convergenceCriteriaSatisfied, # TODO: check this + "primallimit": tc.objectiveLimit, + "duallimit": tc.objectiveLimit, + "sollimit": tc.unknown, + "bestsollimit": tc.unknown, + "restartlimit": tc.unknown, + "optimal": tc.convergenceCriteriaSatisfied, + "infeasible": tc.provenInfeasible, + "unbounded": tc.unbounded, + "inforunbd": tc.infeasibleOrUnbounded, + "terminate": tc.unknown, + } + return ScipDirect._tc_map + + def _scip_lb_ub_from_var(self, var): + if var.is_fixed(): + val = var.value + return val, val + + lb, ub = var.bounds + + if lb is None: + lb = -self._solver_model.infinity() + if ub is None: + ub = self._solver_model.infinity() + + return lb, ub + + def _add_var(self, var): + vtype = self._scip_vtype_from_var(var) + lb, ub = self._scip_lb_ub_from_var(var) + + scip_var = self._solver_model.addVar(lb=lb, ub=ub, vtype=vtype) + + self._pyomo_var_to_solver_var_map[var] = scip_var + return scip_var + + def _add_param(self, p): + vtype = "C" + lb = ub = p.value + scip_var = self._solver_model.addVar(lb=lb, ub=ub, vtype=vtype) + self._pyomo_param_to_solver_param_map[p] = scip_var + return scip_var + + def _add_constraints(self, cons: List[ConstraintData]): + for con in cons: + self._add_constraint(con) + + def _add_sos_constraints(self, cons: List[SOSConstraintData]): + for con in cons: + self._add_sos_constraint(con) + + def _create_solver_model(self, model, config): + timer = config.timer + timer.start('create scip model') + self._clear() + self._solver_model = scip.Model() + timer.start('collect constraints') + cons = list( + model.component_data_objects(Constraint, descend_into=True, active=True) + ) + timer.stop('collect constraints') + timer.start('translate constraints') + self._add_constraints(cons) + timer.stop('translate constraints') + timer.start('sos') + sos = list( + model.component_data_objects(SOSConstraint, descend_into=True, active=True) + ) + self._add_sos_constraints(sos) + timer.stop('sos') + timer.start('get objective') + obj = get_objective(model) + timer.stop('get objective') + timer.start('translate objective') + self._set_objective(obj) + timer.stop('translate objective') + has_obj = obj is not None + solution_loader = ScipSolutionLoader( + solver_model=self._solver_model, + var_map=self._pyomo_var_to_solver_var_map, + con_map=self._pyomo_con_to_solver_con_map, + pyomo_model=model, + opt=self, + ) + timer.stop('create scip model') + return self._solver_model, solution_loader, has_obj + + def _add_constraint(self, con): + scip_expr = self._expr_visitor.walk_expression(con.expr) + scip_con = self._solver_model.addCons(scip_expr) + self._pyomo_con_to_solver_con_map[con] = scip_con + + def _add_sos_constraint(self, con): + level = con.level + if level not in [1, 2]: + raise ValueError( + f"{self.name} does not support SOS level {level} constraints" + ) + + scip_vars = [] + weights = [] + + for v, w in con.get_items(): + if v not in self._pyomo_var_to_solver_var_map: + self._add_var(v) + scip_vars.append(self._pyomo_var_to_solver_var_map[v]) + weights.append(w) + + if level == 1: + scip_cons = self._solver_model.addConsSOS1(scip_vars, weights=weights) + else: + scip_cons = self._solver_model.addConsSOS2(scip_vars, weights=weights) + self._pyomo_con_to_solver_con_map[con] = scip_cons + + def _scip_vtype_from_var(self, var): + """ + This function takes a pyomo variable and returns the appropriate SCIP variable type + + Parameters + ---------- + var: pyomo.core.base.var.Var + The pyomo variable that we want to retrieve the SCIP vtype of + + Returns + ------- + vtype: str + B for Binary, I for Integer, or C for Continuous + """ + if var.is_binary(): + vtype = "B" + elif var.is_integer(): + vtype = "I" + elif var.is_continuous(): + vtype = "C" + else: + raise ValueError(f"Variable domain type is not recognized for {var.domain}") + return vtype + + def _set_objective(self, obj): + if self._obj_var is None: + self._obj_var = self._solver_model.addVar( + lb=-self._solver_model.infinity(), + ub=self._solver_model.infinity(), + vtype="C", + ) + + if self._obj_con is not None: + self._solver_model.delCons(self._obj_con) + + if obj is None: + scip_expr = 0 + sense = "minimize" + else: + scip_expr = self._expr_visitor.walk_expression(obj.expr) + if obj.sense == minimize: + sense = "minimize" + elif obj.sense == maximize: + sense = "maximize" + else: + raise ValueError(f"Objective sense is not recognized: {obj.sense}") + + if sense == "minimize": + self._obj_con = self._solver_model.addCons(self._obj_var >= scip_expr) + else: + self._obj_con = self._solver_model.addCons(self._obj_var <= scip_expr) + + self._solver_model.setObjective(self._obj_var, sense=sense) + self._objective = obj + + def _populate_results(self, scip_model, solution_loader, has_obj, config): + results = Results() + results.solution_loader = solution_loader + results.timing_info.scip_time = scip_model.getSolvingTime() + results.termination_condition = self._get_tc_map().get( + scip_model.getStatus(), TerminationCondition.unknown + ) + + if solution_loader.get_number_of_solutions() > 0: + if ( + results.termination_condition + == TerminationCondition.convergenceCriteriaSatisfied + ): + results.solution_status = SolutionStatus.optimal + else: + results.solution_status = SolutionStatus.feasible + else: + results.solution_status = SolutionStatus.noSolution + + if ( + results.termination_condition + != TerminationCondition.convergenceCriteriaSatisfied + and config.raise_exception_on_nonoptimal_result + ): + raise NoOptimalSolutionError() + + if has_obj: + try: + if ( + scip_model.getNSols() > 0 + and scip_model.getObjVal() < scip_model.infinity() + ): + results.incumbent_objective = scip_model.getObjVal() + else: + results.incumbent_objective = None + except: + results.incumbent_objective = None + try: + results.objective_bound = scip_model.getDualbound() + if results.objective_bound <= -scip_model.infinity(): + results.objective_bound = -math.inf + if results.objective_bound >= scip_model.infinity(): + results.objective_bound = math.inf + except: + if self._objective.sense == minimize: + results.objective_bound = -math.inf + else: + results.objective_bound = math.inf + else: + results.incumbent_objective = None + results.objective_bound = None + + config.timer.start('load solution') + if config.load_solutions: + if solution_loader.get_number_of_solutions() > 0: + solution_loader.load_solution() + else: + raise NoFeasibleSolutionError() + config.timer.stop('load solution') + + results.extra_info['NNodes'] = scip_model.getNNodes() + results.solver_config = config + results.solver_name = self.name + results.solver_version = self.version() + + return results + + def _mipstart(self): + # TODO: it is also possible to specify continuous variables, but + # I think we should have a different option for that + sol = self._solver_model.createPartialSol() + for pyomo_var, scip_var in self._pyomo_var_to_solver_var_map.items(): + if pyomo_var.is_integer(): + sol[scip_var] = pyomo_var.value + self._solver_model.addSol(sol) diff --git a/pyomo/contrib/solver/solvers/scip/scip_persistent.py b/pyomo/contrib/solver/solvers/scip/scip_persistent.py new file mode 100644 index 00000000000..663d26f1675 --- /dev/null +++ b/pyomo/contrib/solver/solvers/scip/scip_persistent.py @@ -0,0 +1,418 @@ +# ____________________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2026 National Technology and Engineering Solutions of Sandia, LLC +# Under the terms of Contract DE-NA0003525 with National Technology and Engineering +# Solutions of Sandia, LLC, the U.S. Government retains certain rights in this +# software. This software is distributed under the 3-clause BSD License. +# ____________________________________________________________________________________ + +from __future__ import annotations + +from typing import List, Mapping, Optional + +from pyomo.common.collections import ComponentMap +from pyomo.common.timing import HierarchicalTimer + +from pyomo.core.base.constraint import ConstraintData +from pyomo.core.base.objective import ObjectiveData +from pyomo.core.base.param import ParamData +from pyomo.core.base.sos import SOSConstraintData +from pyomo.core.base.var import VarData + +from pyomo.contrib.observer.model_observer import ( + Observer, + ModelChangeDetector, + AutoUpdateConfig, + Reason, +) +from pyomo.contrib.solver.common.base import PersistentSolverBase +from pyomo.contrib.solver.common.results import Results + +import pyomo.contrib.solver.solvers.scip.base as scip_base +import pyomo.contrib.solver.solvers.scip.scip_direct as scip_direct + + +class ScipPersistentSolutionLoader(scip_base.ScipSolutionLoader): + def __init__(self, solver_model, var_map, con_map, pyomo_model, opt) -> None: + super().__init__(solver_model, var_map, con_map, pyomo_model, opt) + self._valid = True + + def invalidate(self): + self._valid = False + + def _assert_solution_still_valid(self): + if not self._valid: + raise RuntimeError('The results in the solver are no longer valid.') + + def load_vars(self, vars_to_load: List[VarData] | None = None) -> None: + self._assert_solution_still_valid() + return super().load_vars(vars_to_load) + + def get_vars( + self, vars_to_load: List[VarData] | None = None + ) -> Mapping[VarData, float]: + self._assert_solution_still_valid() + return super().get_vars(vars_to_load) + + def get_number_of_solutions(self) -> int: + self._assert_solution_still_valid() + return super().get_number_of_solutions() + + def get_solution_ids(self) -> List: + self._assert_solution_still_valid() + return super().get_solution_ids() + + def load_import_suffixes(self): + self._assert_solution_still_valid() + super().load_import_suffixes() + + def _set_solution_id(self, solution_id: int) -> int: + self._assert_solution_still_valid() + return super()._set_solution_id(solution_id) + + +class ScipPersistentConfig(scip_base.ScipConfig): + def __init__( + self, + description=None, + doc=None, + implicit=False, + implicit_domain=None, + visibility=0, + ): + scip_base.ScipConfig.__init__( + self, + description=description, + doc=doc, + implicit=implicit, + implicit_domain=implicit_domain, + visibility=visibility, + ) + self.auto_updates: bool = self.declare('auto_updates', AutoUpdateConfig()) + + +class ScipPersistent(scip_direct.ScipDirect, PersistentSolverBase, Observer): + _minimum_version = (5, 5, 0) # this is probably conservative + CONFIG = ScipPersistentConfig() + + def __init__(self, **kwds): + super().__init__(**kwds) + self._pyomo_model = None + self._change_detector = None + self._last_results_object: Optional[Results] = None + self._needs_reopt = False + self._range_constraints = set() + + def _clear(self): + super()._clear() + self._pyomo_model = None + self._change_detector = None + self._needs_reopt = False + self._range_constraints = set() + + def _check_reopt(self): + if self._needs_reopt: + # self._solver_model.freeReoptSolve() # when is it safe to use this one??? + self._solver_model.freeTransform() + self._needs_reopt = False + + def _create_solver_model(self, pyomo_model, config): + if pyomo_model is self._pyomo_model: + self.update(**config) + else: + self.set_instance(pyomo_model, **config) + + solution_loader = ScipPersistentSolutionLoader( + solver_model=self._solver_model, + var_map=self._pyomo_var_to_solver_var_map, + con_map=self._pyomo_con_to_solver_con_map, + pyomo_model=pyomo_model, + opt=self, + ) + + has_obj = self._objective is not None + return self._solver_model, solution_loader, has_obj + + def solve(self, model, **kwds) -> Results: + res = super().solve(model, **kwds) + self._last_results_object = res + self._needs_reopt = True + return res + + def update(self, **kwds): + config = self.config(value=kwds, preserve_implicit=True) + if config.timer is None: + timer = HierarchicalTimer() + else: + timer = config.timer + if self._pyomo_model is None: + raise RuntimeError('must call set_instance or solve before update') + timer.start('update') + self._change_detector.update(timer=timer, **config.auto_updates) + timer.stop('update') + + def set_instance(self, pyomo_model, **kwds): + config = self.config(value=kwds, preserve_implicit=True) + if config.timer is None: + timer = HierarchicalTimer() + else: + timer = config.timer + self._clear() + self._pyomo_model = pyomo_model + self._solver_model = scip_base.scip.Model() + timer.start('set_instance') + self._change_detector = ModelChangeDetector( + model=self._pyomo_model, observers=[self], **config.auto_updates + ) + timer.stop('set_instance') + + def _invalidate_last_results(self): + if self._last_results_object is not None: + self._last_results_object.solution_loader.invalidate() + + def _update_variables(self, variables: Mapping[VarData, Reason]): + new_vars = [] + old_vars = [] + mod_vars = [] + for v, reason in variables.items(): + if reason & Reason.added: + new_vars.append(v) + elif reason & Reason.removed: + old_vars.append(v) + else: + mod_vars.append(v) + + if new_vars: + self._add_variables(new_vars) + if old_vars: + self._remove_variables(old_vars) + if mod_vars: + self._update_vars_for_real(mod_vars) + + def _update_parameters(self, params: Mapping[ParamData, Reason]): + new_params = [] + old_params = [] + mod_params = [] + for p, reason in params.items(): + if reason & Reason.added: + new_params.append(p) + elif reason & Reason.removed: + old_params.append(p) + else: + mod_params.append(p) + + if new_params: + self._add_parameters(new_params) + if old_params: + self._remove_parameters(old_params) + if mod_params: + self._update_params_for_real(mod_params) + + def _update_constraints(self, cons: Mapping[ConstraintData, Reason]): + new_cons = [] + old_cons = [] + for c, reason in cons.items(): + if reason & Reason.added: + new_cons.append(c) + elif reason & Reason.removed: + old_cons.append(c) + elif reason & Reason.expr: + old_cons.append(c) + new_cons.append(c) + + if old_cons: + self._remove_constraints(old_cons) + if new_cons: + self._add_constraints(new_cons) + + def _update_sos_constraints(self, cons: Mapping[SOSConstraintData, Reason]): + new_cons = [] + old_cons = [] + for c, reason in cons.items(): + if reason & Reason.added: + new_cons.append(c) + elif reason & Reason.removed: + old_cons.append(c) + elif reason & Reason.sos_items: + old_cons.append(c) + new_cons.append(c) + + if old_cons: + self._remove_sos_constraints(old_cons) + if new_cons: + self._add_sos_constraints(new_cons) + + def _update_objectives(self, objs: Mapping[ObjectiveData, Reason]): + new_objs = [] + old_objs = [] + for obj, reason in objs.items(): + if reason & Reason.added: + new_objs.append(obj) + elif reason & Reason.removed: + old_objs.append(obj) + elif reason & (Reason.expr | Reason.sense): + old_objs.append(obj) + new_objs.append(obj) + + if old_objs: + self._remove_objectives(old_objs) + if new_objs: + self._add_objectives(new_objs) + + def _add_variables(self, variables: List[VarData]): + self._check_reopt() + self._invalidate_last_results() + for v in variables: + self._add_var(v) + + def _add_parameters(self, params: List[ParamData]): + self._check_reopt() + self._invalidate_last_results() + for p in params: + self._add_param(p) + + def _add_constraints(self, cons: List[ConstraintData]): + self._check_reopt() + self._invalidate_last_results() + for con in cons: + if type(con.expr) is scip_base.RangedExpression: + self._range_constraints.add(con) + super()._add_constraints(cons) + + def _add_sos_constraints(self, cons: List[SOSConstraintData]): + self._check_reopt() + self._invalidate_last_results() + return super()._add_sos_constraints(cons) + + def _add_objectives(self, objs: List[ObjectiveData]): + self._check_reopt() + if len(objs) > 1: + raise NotImplementedError( + 'the persistent interface to scip currently ' + f'only supports single-objective problems; got {len(objs)}: ' + f'{[str(i) for i in objs]}' + ) + + if len(objs) == 0: + return + + obj = objs[0] + + if self._objective is not None: + raise NotImplementedError( + 'the persistent interface to scip currently ' + 'only supports single-objective problems; tried to add ' + f'an objective ({str(obj)}), but there is already an ' + f'active objective ({str(self._objective)})' + ) + + self._invalidate_last_results() + self._set_objective(obj) + + def _remove_objectives(self, objs: List[ObjectiveData]): + self._check_reopt() + for obj in objs: + if obj is not self._objective: + raise RuntimeError( + 'tried to remove an objective that has not been added: ' + f'{str(obj)}' + ) + else: + self._invalidate_last_results() + self._set_objective(None) + + def _remove_constraints(self, cons: List[ConstraintData]): + self._check_reopt() + self._invalidate_last_results() + for con in cons: + scip_con = self._pyomo_con_to_solver_con_map.pop(con) + self._solver_model.delCons(scip_con) + self._range_constraints.discard(con) + + def _remove_sos_constraints(self, cons: List[SOSConstraintData]): + self._check_reopt() + self._invalidate_last_results() + for con in cons: + scip_con = self._pyomo_con_to_solver_con_map.pop(con) + self._solver_model.delCons(scip_con) + + def _remove_variables(self, variables: List[VarData]): + self._check_reopt() + self._invalidate_last_results() + for v in variables: + scip_var = self._pyomo_var_to_solver_var_map.pop(v) + self._solver_model.delVar(scip_var) + + def _remove_parameters(self, params: List[ParamData]): + self._check_reopt() + self._invalidate_last_results() + for p in params: + scip_var = self._pyomo_param_to_solver_param_map.pop(p) + self._solver_model.delVar(scip_var) + + def _update_vars_for_real(self, variables: List[VarData]): + self._check_reopt() + self._invalidate_last_results() + for v in variables: + scip_var = self._pyomo_var_to_solver_var_map[v] + vtype = self._scip_vtype_from_var(v) + lb, ub = self._scip_lb_ub_from_var(v) + self._solver_model.chgVarLb(scip_var, lb) + self._solver_model.chgVarUb(scip_var, ub) + self._solver_model.chgVarType(scip_var, vtype) + + def _update_params_for_real(self, params: List[ParamData]): + self._check_reopt() + self._invalidate_last_results() + for p in params: + scip_var = self._pyomo_param_to_solver_param_map[p] + lb = ub = p.value + self._solver_model.chgVarLb(scip_var, lb) + self._solver_model.chgVarUb(scip_var, ub) + impacted_vars = self._change_detector.get_variables_impacted_by_param(p) + if impacted_vars: + impacted_vars_mapping = ComponentMap( + (v, Reason.bounds) for v in impacted_vars + ) + self._update_variables(impacted_vars_mapping) + impacted_cons = self._change_detector.get_constraints_impacted_by_param(p) + for con in impacted_cons: + if con in self._range_constraints: + self._remove_constraints([con]) + self._add_constraints([con]) + + def add_constraints(self, cons): + if self._change_detector is None: + raise RuntimeError('call set_instance first') + self._change_detector.add_constraints(cons) + + def add_sos_constraints(self, cons): + if self._change_detector is None: + raise RuntimeError('call set_instance first') + self._change_detector.add_sos_constraints(cons) + + def set_objective(self, obj: ObjectiveData): + if self._change_detector is None: + raise RuntimeError('call set_instance first') + self._change_detector.add_objectives([obj]) + + def remove_constraints(self, cons): + if self._change_detector is None: + raise RuntimeError('call set_instance first') + self._change_detector.remove_constraints(cons) + + def remove_sos_constraints(self, cons): + if self._change_detector is None: + raise RuntimeError('call set_instance first') + self._change_detector.remove_sos_constraints(cons) + + def update_variables(self, variables): + if self._change_detector is None: + raise RuntimeError('call set_instance first') + self._change_detector.update_variables(variables) + + def update_parameters(self, params): + if self._change_detector is None: + raise RuntimeError('call set_instance first') + self._change_detector.update_parameters(params) diff --git a/pyomo/contrib/solver/tests/solvers/test_scip_direct.py b/pyomo/contrib/solver/tests/solvers/test_scip_direct.py new file mode 100644 index 00000000000..221039db764 --- /dev/null +++ b/pyomo/contrib/solver/tests/solvers/test_scip_direct.py @@ -0,0 +1,353 @@ +# ____________________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2026 National Technology and Engineering Solutions of Sandia, LLC +# Under the terms of Contract DE-NA0003525 with National Technology and Engineering +# Solutions of Sandia, LLC, the U.S. Government retains certain rights in this +# software. This software is distributed under the 3-clause BSD License. +# ____________________________________________________________________________________ + +import datetime + +import pyomo.common.unittest as unittest +import pyomo.environ as pyo + +from pyomo.common.config import ConfigDict +from pyomo.common.timing import HierarchicalTimer +from pyomo.contrib.solver.common.base import Availability +from pyomo.contrib.solver.common.results import SolutionStatus, TerminationCondition +from pyomo.contrib.solver.common.util import ( + NoFeasibleSolutionError, + NoOptimalSolutionError, + NoSolutionError, +) +from pyomo.contrib.solver.solvers.scip.base import ScipConfig +from pyomo.contrib.solver.solvers.scip.scip_direct import ScipDirect +from pyomo.contrib.solver.tests.solvers.test_gurobi_persistent import ( + create_pmedian_model, +) + +scip_available = ScipDirect().available() + + +@unittest.pytest.mark.solver("scip_direct") +class TestScipDirectConfig(unittest.TestCase): + def test_default_instantiation(self): + config = ScipConfig() + self.assertIsNone(config._description) + self.assertEqual(config._visibility, 0) + self.assertFalse(config.tee) + self.assertTrue(config.load_solutions) + self.assertTrue(config.raise_exception_on_nonoptimal_result) + self.assertFalse(config.symbolic_solver_labels) + self.assertIsNone(config.timer) + self.assertIsNone(config.threads) + self.assertIsNone(config.time_limit) + self.assertIsNone(config.rel_gap) + self.assertIsNone(config.abs_gap) + self.assertFalse(config.warmstart_discrete_vars) + self.assertIsInstance(config.solver_options, ConfigDict) + + def test_custom_instantiation(self): + config = ScipConfig(description="A description") + config.tee = True + config.warmstart_discrete_vars = True + self.assertTrue(config.tee) + self.assertEqual(config._description, "A description") + self.assertTrue(config.warmstart_discrete_vars) + + +@unittest.pytest.mark.solver("scip_direct") +class TestScipDirectInterface(unittest.TestCase): + def test_class_member_list(self): + opt = ScipDirect() + expected_list = [ + 'CONFIG', + 'available', + 'config', + 'api_version', + 'is_persistent', + 'name', + 'solve', + 'version', + ] + method_list = [method for method in dir(opt) if not method.startswith('_')] + self.assertEqual(sorted(expected_list), sorted(method_list)) + + def test_default_instantiation(self): + opt = ScipDirect() + self.assertFalse(opt.is_persistent()) + self.assertEqual(opt.name, 'scip_direct') + self.assertEqual(opt.CONFIG, opt.config) + self.assertIn( + opt.available(), + {Availability.NotFound, Availability.BadVersion, Availability.FullLicense}, + ) + + def test_context_manager(self): + with ScipDirect() as opt: + self.assertFalse(opt.is_persistent()) + self.assertEqual(opt.name, 'scip_direct') + self.assertEqual(opt.CONFIG, opt.config) + + def test_version(self): + opt = ScipDirect() + if opt.available() == Availability.FullLicense: + ver = opt.version() + self.assertIsInstance(ver, tuple) + self.assertGreaterEqual(len(ver), 3) + self.assertTrue(all(isinstance(_, int) for _ in ver)) + + def test_get_tc_map(self): + opt = ScipDirect() + tc_map = opt._get_tc_map() + self.assertEqual( + tc_map["optimal"], TerminationCondition.convergenceCriteriaSatisfied + ) + self.assertEqual(tc_map["timelimit"], TerminationCondition.maxTimeLimit) + self.assertEqual(tc_map["infeasible"], TerminationCondition.provenInfeasible) + self.assertEqual(tc_map["unbounded"], TerminationCondition.unbounded) + self.assertEqual( + tc_map["inforunbd"], TerminationCondition.infeasibleOrUnbounded + ) + + def test_scip_vtype_from_var(self): + m = pyo.ConcreteModel() + m.b = pyo.Var(within=pyo.Binary) + m.i = pyo.Var(within=pyo.Integers) + m.c = pyo.Var(within=pyo.Reals) + + opt = ScipDirect() + self.assertEqual(opt._scip_vtype_from_var(m.b), "B") + self.assertEqual(opt._scip_vtype_from_var(m.i), "I") + self.assertEqual(opt._scip_vtype_from_var(m.c), "C") + + +@unittest.skipIf(not scip_available, "SCIP is not available") +@unittest.pytest.mark.solver("scip_direct") +class TestScipDirect(unittest.TestCase): + def create_lp_model(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, None), initialize=0) + m.y = pyo.Var(bounds=(0, None), initialize=0) + m.obj = pyo.Objective(expr=m.x + 2 * m.y) + m.c = pyo.Constraint(expr=m.x + m.y >= 1) + return m + + def create_feasible_model_no_objective(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, None), initialize=0) + m.c = pyo.Constraint(expr=m.x >= 1) + return m + + def create_infeasible_model(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, 1), initialize=0) + m.obj = pyo.Objective(expr=m.x) + m.c = pyo.Constraint(expr=m.x >= 2) + return m + + def create_sos_model(self, level): + m = pyo.ConcreteModel() + m.I = pyo.RangeSet(3) + m.x = pyo.Var(m.I, bounds=(0, 1)) + m.obj = pyo.Objective(expr=sum(i * m.x[i] for i in m.I)) + m.c = pyo.Constraint(expr=sum(m.x[i] for i in m.I) == 1) + m.sos = pyo.SOSConstraint(var=m.x, sos=level) + return m + + def test_solve(self): + m = self.create_lp_model() + opt = ScipDirect() + res = opt.solve(m) + + self.assertEqual( + res.termination_condition, TerminationCondition.convergenceCriteriaSatisfied + ) + self.assertEqual(res.solution_status, SolutionStatus.optimal) + self.assertAlmostEqual(m.x.value, 1) + self.assertAlmostEqual(m.y.value, 0) + self.assertAlmostEqual(res.incumbent_objective, 1) + self.assertIsNotNone(res.objective_bound) + self.assertEqual(res.solver_name, 'scip_direct') + self.assertIsInstance(res.solver_version, tuple) + self.assertIsNotNone(res.solver_log) + self.assertIsNotNone(res.timing_info.scip_time) + self.assertEqual(res.timing_info.start_timestamp.tzinfo, datetime.timezone.utc) + self.assertGreaterEqual(res.timing_info.wall_time, 0) + self.assertIn('NNodes', res.extra_info) + + def test_solve_load_solutions_false(self): + m = self.create_lp_model() + opt = ScipDirect() + res = opt.solve(m, load_solutions=False) + + self.assertEqual( + res.termination_condition, TerminationCondition.convergenceCriteriaSatisfied + ) + self.assertEqual(res.solution_status, SolutionStatus.optimal) + + self.assertAlmostEqual(m.x.value, 0) + self.assertAlmostEqual(m.y.value, 0) + + self.assertEqual(res.solution_loader.get_number_of_solutions(), 1) + self.assertEqual(res.solution_loader.get_solution_ids(), [0]) + + vals = res.solution_loader.get_vars() + self.assertAlmostEqual(vals[m.x], 1) + self.assertAlmostEqual(vals[m.y], 0) + + res.solution_loader.load_vars() + self.assertAlmostEqual(m.x.value, 1) + self.assertAlmostEqual(m.y.value, 0) + + def test_no_objective(self): + m = self.create_feasible_model_no_objective() + opt = ScipDirect() + res = opt.solve(m) + + self.assertEqual( + res.termination_condition, TerminationCondition.convergenceCriteriaSatisfied + ) + self.assertEqual(res.solution_status, SolutionStatus.optimal) + self.assertIsNone(res.incumbent_objective) + self.assertIsNone(res.objective_bound) + + def test_infeasible_no_exception(self): + m = self.create_infeasible_model() + opt = ScipDirect() + res = opt.solve( + m, load_solutions=False, raise_exception_on_nonoptimal_result=False + ) + + self.assertEqual( + res.termination_condition, TerminationCondition.provenInfeasible + ) + self.assertEqual(res.solution_status, SolutionStatus.noSolution) + self.assertIsNone(res.incumbent_objective) + self.assertEqual(res.solution_loader.get_number_of_solutions(), 0) + with self.assertRaises(NoSolutionError): + res.solution_loader.get_vars() + + def test_infeasible_raises_no_optimal_solution_error(self): + m = self.create_infeasible_model() + opt = ScipDirect() + with self.assertRaises(NoOptimalSolutionError): + opt.solve(m, load_solutions=False) + + def test_infeasible_raises_no_feasible_solution_error(self): + m = self.create_infeasible_model() + opt = ScipDirect() + with self.assertRaises(NoFeasibleSolutionError): + opt.solve( + m, load_solutions=True, raise_exception_on_nonoptimal_result=False + ) + + def test_timer(self): + m = self.create_lp_model() + timer = HierarchicalTimer() + opt = ScipDirect() + res = opt.solve(m, timer=timer) + self.assertIs(res.timing_info.timer, timer) + + def test_fixed_var(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(-10, 10), initialize=0) + m.y = pyo.Var(bounds=(-10, 10), initialize=0) + m.x.fix(2) + m.obj = pyo.Objective(expr=m.y) + m.c = pyo.Constraint(expr=m.y >= m.x + 1) + + opt = ScipDirect() + res = opt.solve(m) + + self.assertEqual(res.solution_status, SolutionStatus.optimal) + self.assertAlmostEqual(m.x.value, 2) + self.assertAlmostEqual(m.y.value, 3) + + def test_sos1(self): + m = self.create_sos_model(level=1) + opt = ScipDirect() + res = opt.solve(m) + + self.assertEqual( + res.termination_condition, TerminationCondition.convergenceCriteriaSatisfied + ) + self.assertEqual(res.solution_status, SolutionStatus.optimal) + + def test_sos2(self): + m = self.create_sos_model(level=2) + opt = ScipDirect() + res = opt.solve(m) + + self.assertEqual( + res.termination_condition, TerminationCondition.convergenceCriteriaSatisfied + ) + self.assertEqual(res.solution_status, SolutionStatus.optimal) + + def test_bad_sos_level(self): + m = self.create_sos_model(level=1) + m.del_component(m.sos) + m.sos = pyo.SOSConstraint(var=m.x, sos=3) + + opt = ScipDirect() + with self.assertRaisesRegex(ValueError, "does not support SOS level 3"): + opt.solve(m) + + def test_warmstart_discrete_vars(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(within=pyo.Binary, initialize=1) + m.y = pyo.Var(within=pyo.Binary, initialize=0) + m.obj = pyo.Objective(expr=m.x + 2 * m.y, sense=pyo.maximize) + m.c = pyo.Constraint(expr=m.x + m.y <= 1) + + opt = ScipDirect() + res = opt.solve(m, warmstart_discrete_vars=True) + + self.assertEqual( + res.termination_condition, TerminationCondition.convergenceCriteriaSatisfied + ) + self.assertEqual(res.solution_status, SolutionStatus.optimal) + self.assertAlmostEqual(m.x.value, 0) + self.assertAlmostEqual(m.y.value, 1) + + def test_multiple_solutions(self): + m = create_pmedian_model() + + # The solutions found by scip may change from version to version. + # Let's warmstart scip with a suboptimal solution to ensure we + # have at least 2 solutions. + + init_sol = {1, 2, 3} + for k, y in m.y.items(): + if k in init_sol: + y.value = 1 + else: + y.value = 0 + + opt = ScipDirect() + opt.config.warmstart_discrete_vars = True + opt.config.solver_options['limits/maxsol'] = 100000 + opt.config.solver_options['heuristics/completesol/maxunknownrate'] = 1.0 + res = opt.solve(m, load_solutions=True) + num_solutions = res.solution_loader.get_number_of_solutions() + self.assertGreaterEqual(num_solutions, 2) + + # the best solution + self.assertAlmostEqual(pyo.value(m.obj.expr), 6.431184939357673) + sol = {3, 6, 9} + for k, v in m.y.items(): + if k in sol: + self.assertAlmostEqual(v.value, 1) + else: + self.assertAlmostEqual(v.value, 0) + + # the worst solution that we used to warmstart + res.solution_loader.solution(num_solutions - 1).load_vars() + self.assertAlmostEqual(pyo.value(m.obj.expr), 7.607295680844689) + sol = {1, 2, 3} + for k, v in m.y.items(): + if k in sol: + self.assertAlmostEqual(v.value, 1) + else: + self.assertAlmostEqual(v.value, 0) diff --git a/pyomo/contrib/solver/tests/solvers/test_scip_persistent.py b/pyomo/contrib/solver/tests/solvers/test_scip_persistent.py new file mode 100644 index 00000000000..aa17d5c3495 --- /dev/null +++ b/pyomo/contrib/solver/tests/solvers/test_scip_persistent.py @@ -0,0 +1,291 @@ +# ____________________________________________________________________________________ +# +# Pyomo: Python Optimization Modeling Objects +# Copyright (c) 2008-2026 National Technology and Engineering Solutions of Sandia, LLC +# Under the terms of Contract DE-NA0003525 with National Technology and Engineering +# Solutions of Sandia, LLC, the U.S. Government retains certain rights in this +# software. This software is distributed under the 3-clause BSD License. +# ____________________________________________________________________________________ + +import pyomo.common.unittest as unittest +import pyomo.environ as pyo + +from pyomo.common.timing import HierarchicalTimer +from pyomo.contrib.solver.common.base import Availability +from pyomo.contrib.solver.common.results import SolutionStatus, TerminationCondition +from pyomo.contrib.solver.common.util import ( + NoFeasibleSolutionError, + NoOptimalSolutionError, +) +from pyomo.contrib.solver.solvers.scip.scip_persistent import ( + ScipPersistent, + ScipPersistentConfig, +) + +scip_available = ScipPersistent().available() + + +@unittest.pytest.mark.solver("scip_persistent") +class TestScipPersistentConfig(unittest.TestCase): + def test_default_instantiation(self): + config = ScipPersistentConfig() + self.assertIsNone(config._description) + self.assertEqual(config._visibility, 0) + self.assertFalse(config.tee) + self.assertTrue(config.load_solutions) + self.assertTrue(config.raise_exception_on_nonoptimal_result) + self.assertFalse(config.symbolic_solver_labels) + self.assertIsNone(config.timer) + self.assertIsNone(config.threads) + self.assertIsNone(config.time_limit) + self.assertIsNone(config.rel_gap) + self.assertIsNone(config.abs_gap) + self.assertFalse(config.warmstart_discrete_vars) + self.assertTrue(hasattr(config, 'auto_updates')) + + def test_custom_instantiation(self): + config = ScipPersistentConfig(description="A description") + config.tee = True + config.warmstart_discrete_vars = True + self.assertTrue(config.tee) + self.assertEqual(config._description, "A description") + self.assertTrue(config.warmstart_discrete_vars) + + +@unittest.pytest.mark.solver("scip_persistent") +class TestScipPersistentInterface(unittest.TestCase): + def test_default_instantiation(self): + opt = ScipPersistent() + self.assertTrue(opt.is_persistent()) + self.assertEqual(opt.name, 'scip_persistent') + self.assertEqual(opt.CONFIG, opt.config) + self.assertIn( + opt.available(), + {Availability.NotFound, Availability.BadVersion, Availability.FullLicense}, + ) + + def test_context_manager(self): + with ScipPersistent() as opt: + self.assertTrue(opt.is_persistent()) + self.assertEqual(opt.name, 'scip_persistent') + self.assertEqual(opt.CONFIG, opt.config) + + def test_update_before_set_instance_raises(self): + opt = ScipPersistent() + with self.assertRaisesRegex( + RuntimeError, 'must call set_instance or solve before update' + ): + opt.update() + + def test_add_constraints_before_set_instance_raises(self): + opt = ScipPersistent() + m = pyo.ConcreteModel() + m.x = pyo.Var() + m.c = pyo.Constraint(expr=m.x >= 1) + with self.assertRaisesRegex(RuntimeError, 'call set_instance first'): + opt.add_constraints([m.c]) + + +@unittest.skipIf(not scip_available, "SCIP is not available") +@unittest.pytest.mark.solver("scip_persistent") +class TestScipPersistent(unittest.TestCase): + def create_lp_model(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, None), initialize=0) + m.y = pyo.Var(bounds=(0, None), initialize=0) + m.obj = pyo.Objective(expr=m.x + 2 * m.y) + m.c = pyo.Constraint(expr=m.x + m.y >= 1) + return m + + def create_range_model(self): + m = pyo.ConcreteModel() + m.x = pyo.Var() + m.xl = pyo.Param(initialize=-1, mutable=True) + m.xu = pyo.Param(initialize=1, mutable=True) + m.c = pyo.Constraint(expr=pyo.inequality(m.xl, m.x, m.xu)) + m.obj = pyo.Objective(expr=m.x) + return m + + def test_set_instance_and_solve(self): + m = self.create_lp_model() + opt = ScipPersistent() + opt.set_instance(m) + res = opt.solve(m) + + self.assertEqual( + res.termination_condition, TerminationCondition.convergenceCriteriaSatisfied + ) + self.assertEqual(res.solution_status, SolutionStatus.optimal) + self.assertAlmostEqual(m.x.value, 1) + self.assertAlmostEqual(m.y.value, 0) + + def test_solve_twice_same_instance(self): + m = self.create_lp_model() + opt = ScipPersistent() + + res = opt.solve(m) + self.assertEqual(res.solution_status, SolutionStatus.optimal) + self.assertAlmostEqual(m.x.value, 1) + self.assertAlmostEqual(m.y.value, 0) + + m.c.set_value(m.x + m.y >= 2) + res = opt.solve(m) + self.assertEqual(res.solution_status, SolutionStatus.optimal) + self.assertAlmostEqual(m.x.value, 2) + self.assertAlmostEqual(m.y.value, 0) + + def test_load_solutions_false(self): + m = self.create_lp_model() + opt = ScipPersistent() + res = opt.solve(m, load_solutions=False) + + self.assertEqual(res.solution_status, SolutionStatus.optimal) + self.assertAlmostEqual(m.x.value, 0) + self.assertAlmostEqual(m.y.value, 0) + + vals = res.solution_loader.get_vars() + self.assertAlmostEqual(vals[m.x], 1) + self.assertAlmostEqual(vals[m.y], 0) + + res.solution_loader.load_vars() + self.assertAlmostEqual(m.x.value, 1) + self.assertAlmostEqual(m.y.value, 0) + + def test_solution_loader_invalidated_after_update(self): + m = self.create_lp_model() + opt = ScipPersistent() + res = opt.solve(m, load_solutions=False) + + m.c.set_value(m.x + m.y >= 2) + opt.update() + + with self.assertRaisesRegex( + RuntimeError, 'The results in the solver are no longer valid.' + ): + res.solution_loader.get_vars() + + def test_range_constraint_mutable_params(self): + m = self.create_range_model() + opt = ScipPersistent() + opt.set_instance(m) + + res = opt.solve(m) + self.assertEqual(res.solution_status, SolutionStatus.optimal) + self.assertAlmostEqual(m.x.value, -1) + + m.xl.value = -3 + res = opt.solve(m) + self.assertAlmostEqual(m.x.value, -3) + + del m.obj + m.obj = pyo.Objective(expr=m.x, sense=pyo.maximize) + res = opt.solve(m) + self.assertAlmostEqual(m.x.value, 1) + + m.xu.value = 3 + res = opt.solve(m) + self.assertAlmostEqual(m.x.value, 3) + + def test_add_remove_constraints(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(-10, 10)) + m.y = pyo.Var() + m.obj = pyo.Objective(expr=m.y) + m.c1 = pyo.Constraint(expr=m.y >= 2 * m.x + 1) + + opt = ScipPersistent() + opt.set_instance(m) + + res = opt.solve(m) + self.assertAlmostEqual(m.x.value, -10) + self.assertAlmostEqual(m.y.value, -19) + + m.c2 = pyo.Constraint(expr=m.y >= -m.x + 1) + opt.add_constraints([m.c2]) + + res = opt.solve(m) + self.assertAlmostEqual(m.x.value, 0) + self.assertAlmostEqual(m.y.value, 1) + + opt.remove_constraints([m.c2]) + m.del_component(m.c2) + + res = opt.solve(m) + self.assertAlmostEqual(m.x.value, -10) + self.assertAlmostEqual(m.y.value, -19) + + def test_update_variables_manual(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(-1, 1)) + m.obj = pyo.Objective(expr=m.x) + + opt = ScipPersistent() + opt.config.auto_updates.update_vars = False + opt.set_instance(m) + + res = opt.solve(m) + self.assertAlmostEqual(m.x.value, -1) + + m.x.setlb(-3) + opt.update_variables([m.x]) + + res = opt.solve(m) + self.assertAlmostEqual(m.x.value, -3) + + def test_update_parameters_manual(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(-10, 10)) + m.p = pyo.Param(initialize=1, mutable=True) + m.obj = pyo.Objective(expr=m.x) + m.c = pyo.Constraint(expr=m.x >= m.p) + + opt = ScipPersistent() + opt.config.auto_updates.update_parameters = False + opt.set_instance(m) + + res = opt.solve(m) + self.assertAlmostEqual(m.x.value, 1) + + m.p.value = 3 + opt.update_parameters([m.p]) + + res = opt.solve(m) + self.assertAlmostEqual(m.x.value, 3) + + def test_timer(self): + m = self.create_lp_model() + timer = HierarchicalTimer() + opt = ScipPersistent() + res = opt.solve(m, timer=timer) + self.assertIs(res.timing_info.timer, timer) + + def test_infeasible_no_exception(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, 1)) + m.obj = pyo.Objective(expr=m.x) + m.c = pyo.Constraint(expr=m.x >= 2) + + opt = ScipPersistent() + res = opt.solve( + m, load_solutions=False, raise_exception_on_nonoptimal_result=False + ) + + self.assertEqual( + res.termination_condition, TerminationCondition.provenInfeasible + ) + self.assertEqual(res.solution_status, SolutionStatus.noSolution) + + def test_infeasible_raises(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, 1)) + m.obj = pyo.Objective(expr=m.x) + m.c = pyo.Constraint(expr=m.x >= 2) + + opt = ScipPersistent() + with self.assertRaises(NoOptimalSolutionError): + opt.solve(m, load_solutions=False) + + with self.assertRaises(NoFeasibleSolutionError): + opt.solve( + m, load_solutions=True, raise_exception_on_nonoptimal_result=False + ) diff --git a/pyomo/contrib/solver/tests/solvers/test_solvers.py b/pyomo/contrib/solver/tests/solvers/test_solvers.py index 686c45adb74..723aa39433c 100644 --- a/pyomo/contrib/solver/tests/solvers/test_solvers.py +++ b/pyomo/contrib/solver/tests/solvers/test_solvers.py @@ -25,6 +25,8 @@ SolutionStatus, TerminationCondition, ) +from pyomo.contrib.solver.solvers.scip.scip_direct import ScipDirect +from pyomo.contrib.solver.solvers.scip.scip_persistent import ScipPersistent from pyomo.contrib.solver.common.util import ( NoDualsError, NoReducedCostsError, @@ -74,6 +76,8 @@ def param_as_standalone_func(cls, p, func, name): ('gurobi_direct_minlp', GurobiDirectMINLP), ('ipopt', Ipopt), ('highs', Highs), + ('scip_direct', ScipDirect), + ('scip_persistent', ScipPersistent), ('gams', GAMS), ('knitro_direct', KnitroDirectSolver), ] @@ -82,27 +86,47 @@ def param_as_standalone_func(cls, p, func, name): ('gurobi_direct', GurobiDirect), ('gurobi_direct_minlp', GurobiDirectMINLP), ('highs', Highs), + ('scip_direct', ScipDirect), + ('scip_persistent', ScipPersistent), ('knitro_direct', KnitroDirectSolver), ] nlp_solvers = [ ('gurobi_direct_minlp', GurobiDirectMINLP), ('ipopt', Ipopt), + ('scip_direct', ScipDirect), + ('scip_persistent', ScipPersistent), ('knitro_direct', KnitroDirectSolver), ] qcp_solvers = [ ('gurobi_persistent', GurobiPersistent), ('gurobi_direct_minlp', GurobiDirectMINLP), ('ipopt', Ipopt), + ('scip_direct', ScipDirect), + ('scip_persistent', ScipPersistent), ('knitro_direct', KnitroDirectSolver), ] qp_solvers = qcp_solvers + [("highs", Highs)] miqcqp_solvers = [ ('gurobi_direct_minlp', GurobiDirectMINLP), ('gurobi_persistent', GurobiPersistent), + ('scip_direct', ScipDirect), + ('scip_persistent', ScipPersistent), ('knitro_direct', KnitroDirectSolver), ] nl_solvers = [('ipopt', Ipopt)] nl_solvers_set = {i[0] for i in nl_solvers} +dual_solvers = [ + ('gurobi_persistent', GurobiPersistent), + ('gurobi_direct', GurobiDirect), + ('gurobi_direct_minlp', GurobiDirectMINLP), + ('ipopt', Ipopt), + ('highs', Highs), +] +sos_solvers = [ + ('gurobi_persistent', GurobiPersistent), + ('scip_direct', ScipDirect), + ('scip_persistent', ScipPersistent), +] def _load_tests(solver_list): @@ -129,7 +153,7 @@ def test_all_solvers_list(): class TestDualSignConvention(unittest.TestCase): - @mark_parameterized.expand(input=_load_tests(all_solvers)) + @mark_parameterized.expand(input=_load_tests(dual_solvers)) def test_equality(self, name: str, opt_class: Type[SolverBase], use_presolve: bool): opt: SolverBase = opt_class() if not opt.available(): @@ -181,7 +205,7 @@ def test_equality(self, name: str, opt_class: Type[SolverBase], use_presolve: bo self.assertAlmostEqual(duals[m.c1], 0) self.assertAlmostEqual(duals[m.c2], -1) - @mark_parameterized.expand(input=_load_tests(all_solvers)) + @mark_parameterized.expand(input=_load_tests(dual_solvers)) def test_inequality( self, name: str, opt_class: Type[SolverBase], use_presolve: bool ): @@ -243,7 +267,7 @@ def test_inequality( self.assertAlmostEqual(duals[m.c1], 0.5) self.assertAlmostEqual(duals[m.c2], 0.5) - @mark_parameterized.expand(input=_load_tests(all_solvers)) + @mark_parameterized.expand(input=_load_tests(dual_solvers)) def test_bounds(self, name: str, opt_class: Type[SolverBase], use_presolve: bool): opt: SolverBase = opt_class() if not opt.available(): @@ -298,7 +322,7 @@ def test_bounds(self, name: str, opt_class: Type[SolverBase], use_presolve: bool rc = res.solution_loader.get_reduced_costs() self.assertAlmostEqual(rc[m.x], -1) - @mark_parameterized.expand(input=_load_tests(all_solvers)) + @mark_parameterized.expand(input=_load_tests(dual_solvers)) def test_range(self, name: str, opt_class: Type[SolverBase], use_presolve: bool): opt: SolverBase = opt_class() if not opt.available(): @@ -350,7 +374,7 @@ def test_range(self, name: str, opt_class: Type[SolverBase], use_presolve: bool) self.assertAlmostEqual(duals[m.c1], -0.5) self.assertAlmostEqual(duals[m.c2], -0.5) - @mark_parameterized.expand(input=_load_tests(all_solvers)) + @mark_parameterized.expand(input=_load_tests(dual_solvers)) def test_equality_max( self, name: str, opt_class: Type[SolverBase], use_presolve: bool ): @@ -404,7 +428,7 @@ def test_equality_max( self.assertAlmostEqual(duals[m.c1], 0) self.assertAlmostEqual(duals[m.c2], 1) - @mark_parameterized.expand(input=_load_tests(all_solvers)) + @mark_parameterized.expand(input=_load_tests(dual_solvers)) def test_inequality_max( self, name: str, opt_class: Type[SolverBase], use_presolve: bool ): @@ -466,7 +490,7 @@ def test_inequality_max( self.assertAlmostEqual(duals[m.c1], -0.5) self.assertAlmostEqual(duals[m.c2], -0.5) - @mark_parameterized.expand(input=_load_tests(all_solvers)) + @mark_parameterized.expand(input=_load_tests(dual_solvers)) def test_bounds_max( self, name: str, opt_class: Type[SolverBase], use_presolve: bool ): @@ -523,7 +547,7 @@ def test_bounds_max( rc = res.solution_loader.get_reduced_costs() self.assertAlmostEqual(rc[m.x], 1) - @mark_parameterized.expand(input=_load_tests(all_solvers)) + @mark_parameterized.expand(input=_load_tests(dual_solvers)) def test_range_max( self, name: str, opt_class: Type[SolverBase], use_presolve: bool ): @@ -745,16 +769,18 @@ def test_range_constraint( res = opt.solve(m) self.assertEqual(res.solution_status, SolutionStatus.optimal) self.assertAlmostEqual(m.x.value, -1) - duals = res.solution_loader.get_duals() - self.assertAlmostEqual(duals[m.c], 1) + if (name, opt_class) in dual_solvers: + duals = res.solution_loader.get_duals() + self.assertAlmostEqual(duals[m.c], 1) m.obj.sense = pyo.maximize res = opt.solve(m) self.assertEqual(res.solution_status, SolutionStatus.optimal) self.assertAlmostEqual(m.x.value, 1) - duals = res.solution_loader.get_duals() - self.assertAlmostEqual(duals[m.c], 1) + if (name, opt_class) in dual_solvers: + duals = res.solution_loader.get_duals() + self.assertAlmostEqual(duals[m.c], 1) - @mark_parameterized.expand(input=_load_tests(all_solvers)) + @mark_parameterized.expand(input=_load_tests(dual_solvers)) def test_reduced_costs( self, name: str, opt_class: Type[SolverBase], use_presolve: bool ): @@ -783,7 +809,7 @@ def test_reduced_costs( self.assertAlmostEqual(rc[m.x], -3) self.assertAlmostEqual(rc[m.y], -4) - @mark_parameterized.expand(input=_load_tests(all_solvers)) + @mark_parameterized.expand(input=_load_tests(dual_solvers)) def test_reduced_costs2( self, name: str, opt_class: Type[SolverBase], use_presolve: bool ): @@ -849,9 +875,10 @@ def test_param_changes( else: bound = res.objective_bound self.assertTrue(bound <= m.y.value) - duals = res.solution_loader.get_duals() - self.assertAlmostEqual(duals[m.c1], (1 + a1 / (a2 - a1))) - self.assertAlmostEqual(duals[m.c2], a1 / (a2 - a1)) + if (name, opt_class) in dual_solvers: + duals = res.solution_loader.get_duals() + self.assertAlmostEqual(duals[m.c1], (1 + a1 / (a2 - a1))) + self.assertAlmostEqual(duals[m.c2], a1 / (a2 - a1)) @mark_parameterized.expand(input=_load_tests(all_solvers)) def test_immutable_param( @@ -896,9 +923,10 @@ def test_immutable_param( else: bound = res.objective_bound self.assertTrue(bound <= m.y.value) - duals = res.solution_loader.get_duals() - self.assertAlmostEqual(duals[m.c1], (1 + a1 / (a2 - a1))) - self.assertAlmostEqual(duals[m.c2], a1 / (a2 - a1)) + if (name, opt_class) in dual_solvers: + duals = res.solution_loader.get_duals() + self.assertAlmostEqual(duals[m.c1], (1 + a1 / (a2 - a1))) + self.assertAlmostEqual(duals[m.c2], a1 / (a2 - a1)) @mark_parameterized.expand(input=_load_tests(all_solvers)) def test_equality(self, name: str, opt_class: Type[SolverBase], use_presolve: bool): @@ -912,6 +940,8 @@ def test_equality(self, name: str, opt_class: Type[SolverBase], use_presolve: bo check_duals = False else: opt.config.writer_config.linear_presolve = False + if (name, opt_class) not in dual_solvers: + check_duals = False m = pyo.ConcreteModel() m.x = pyo.Var() m.y = pyo.Var() @@ -1003,6 +1033,8 @@ def test_no_objective( opt.config.writer_config.linear_presolve = True else: opt.config.writer_config.linear_presolve = False + if (name, opt_class) not in dual_solvers: + check_duals = False m = pyo.ConcreteModel() m.x = pyo.Var() m.y = pyo.Var() @@ -1064,9 +1096,10 @@ def test_add_remove_cons( else: bound = res.objective_bound self.assertTrue(bound <= m.y.value) - duals = res.solution_loader.get_duals() - self.assertAlmostEqual(duals[m.c1], -(1 + a1 / (a2 - a1))) - self.assertAlmostEqual(duals[m.c2], a1 / (a2 - a1)) + if (name, opt_class) in dual_solvers: + duals = res.solution_loader.get_duals() + self.assertAlmostEqual(duals[m.c1], -(1 + a1 / (a2 - a1))) + self.assertAlmostEqual(duals[m.c2], a1 / (a2 - a1)) m.c3 = pyo.Constraint(expr=m.y >= a3 * m.x + b3) res = opt.solve(m) @@ -1075,10 +1108,11 @@ def test_add_remove_cons( self.assertAlmostEqual(m.y.value, a1 * (b3 - b1) / (a1 - a3) + b1) self.assertAlmostEqual(res.incumbent_objective, m.y.value) self.assertTrue(res.objective_bound is None or res.objective_bound <= m.y.value) - duals = res.solution_loader.get_duals() - self.assertAlmostEqual(duals[m.c1], -(1 + a1 / (a3 - a1))) - self.assertAlmostEqual(duals[m.c2], 0) - self.assertAlmostEqual(duals[m.c3], a1 / (a3 - a1)) + if (name, opt_class) in dual_solvers: + duals = res.solution_loader.get_duals() + self.assertAlmostEqual(duals[m.c1], -(1 + a1 / (a3 - a1))) + self.assertAlmostEqual(duals[m.c2], 0) + self.assertAlmostEqual(duals[m.c3], a1 / (a3 - a1)) del m.c3 res = opt.solve(m) @@ -1087,9 +1121,10 @@ def test_add_remove_cons( self.assertAlmostEqual(m.y.value, a1 * (b2 - b1) / (a1 - a2) + b1) self.assertAlmostEqual(res.incumbent_objective, m.y.value) self.assertTrue(res.objective_bound is None or res.objective_bound <= m.y.value) - duals = res.solution_loader.get_duals() - self.assertAlmostEqual(duals[m.c1], -(1 + a1 / (a2 - a1))) - self.assertAlmostEqual(duals[m.c2], a1 / (a2 - a1)) + if (name, opt_class) in dual_solvers: + duals = res.solution_loader.get_duals() + self.assertAlmostEqual(duals[m.c1], -(1 + a1 / (a2 - a1))) + self.assertAlmostEqual(duals[m.c2], a1 / (a2 - a1)) @mark_parameterized.expand(input=_load_tests(all_solvers)) def test_results_infeasible( @@ -1138,14 +1173,16 @@ def test_results_infeasible( NoSolutionError, '.*does not currently have a valid solution.*' ): res.solution_loader.load_vars() - with self.assertRaisesRegex( - NoDualsError, '.*does not currently have valid duals.*' - ): - res.solution_loader.get_duals() - with self.assertRaisesRegex( - NoReducedCostsError, '.*does not currently have valid reduced costs.*' - ): - res.solution_loader.get_reduced_costs() + if (name, opt_class) in dual_solvers: + with self.assertRaisesRegex( + NoDualsError, '.*does not currently have valid duals.*' + ): + res.solution_loader.get_duals() + with self.assertRaisesRegex( + NoReducedCostsError, + '.*does not currently have valid reduced costs.*', + ): + res.solution_loader.get_reduced_costs() @mark_parameterized.expand(input=_load_tests(all_solvers)) def test_trivial_constraints( @@ -1206,7 +1243,7 @@ def test_trivial_constraints( self.assertIn(res.termination_condition, acceptable_termination_conditions) self.assertIsNone(res.incumbent_objective) - @mark_parameterized.expand(input=_load_tests(all_solvers)) + @mark_parameterized.expand(input=_load_tests(dual_solvers)) def test_duals(self, name: str, opt_class: Type[SolverBase], use_presolve: bool): opt: SolverBase = opt_class() if not opt.available(): @@ -1255,13 +1292,13 @@ def test_mutable_quadratic_coefficient( m.c = pyo.Constraint(expr=m.y >= (m.a * m.x + m.b) ** 2) res = opt.solve(m) - self.assertAlmostEqual(m.x.value, 0.41024548525899274, 4) - self.assertAlmostEqual(m.y.value, 0.34781038127030117, 4) + self.assertAlmostEqual(m.x.value, 0.41024548525899274, 3) + self.assertAlmostEqual(m.y.value, 0.34781038127030117, 3) m.a.value = 2 m.b.value = -0.5 res = opt.solve(m) - self.assertAlmostEqual(m.x.value, 0.10256137418973625, 4) - self.assertAlmostEqual(m.y.value, 0.0869525991355825, 4) + self.assertAlmostEqual(m.x.value, 0.10256137418973625, 3) + self.assertAlmostEqual(m.y.value, 0.0869525991355825, 3) @mark_parameterized.expand(input=_load_tests(qcp_solvers)) def test_mutable_quadratic_objective_qcp( @@ -1286,14 +1323,14 @@ def test_mutable_quadratic_objective_qcp( m.ccon = pyo.Constraint(expr=m.y >= (m.a * m.x + m.b) ** 2) res = opt.solve(m) - self.assertAlmostEqual(m.x.value, 0.2719178742733325, 4) - self.assertAlmostEqual(m.y.value, 0.5301035741688002, 4) + self.assertAlmostEqual(m.x.value, 0.2719178742733325, 3) + self.assertAlmostEqual(m.y.value, 0.5301035741688002, 3) m.c.value = 3.5 m.d.value = -1 res = opt.solve(m) - self.assertAlmostEqual(m.x.value, 0.6962249634573562, 4) - self.assertAlmostEqual(m.y.value, 0.09227926676152151, 4) + self.assertAlmostEqual(m.x.value, 0.6962249634573562, 3) + self.assertAlmostEqual(m.y.value, 0.09227926676152151, 3) @mark_parameterized.expand(input=_load_tests(qp_solvers)) def test_mutable_quadratic_objective_qp( @@ -1587,9 +1624,10 @@ def test_mutable_param_with_range( res.objective_bound is None or res.objective_bound <= m.y.value + 1e-12 ) - duals = res.solution_loader.get_duals() - self.assertAlmostEqual(duals[m.con1], (1 + a1 / (a2 - a1)), 6) - self.assertAlmostEqual(duals[m.con2], -a1 / (a2 - a1), 6) + if (name, opt_class) in dual_solvers: + duals = res.solution_loader.get_duals() + self.assertAlmostEqual(duals[m.con1], (1 + a1 / (a2 - a1)), 6) + self.assertAlmostEqual(duals[m.con2], -a1 / (a2 - a1), 6) else: self.assertAlmostEqual(m.x.value, (c2 - c1) / (a1 - a2), 6) self.assertAlmostEqual(m.y.value, a1 * (c2 - c1) / (a1 - a2) + c1, 6) @@ -1598,9 +1636,10 @@ def test_mutable_param_with_range( res.objective_bound is None or res.objective_bound >= m.y.value - 1e-12 ) - duals = res.solution_loader.get_duals() - self.assertAlmostEqual(duals[m.con1], (1 + a1 / (a2 - a1)), 6) - self.assertAlmostEqual(duals[m.con2], -a1 / (a2 - a1), 6) + if (name, opt_class) in dual_solvers: + duals = res.solution_loader.get_duals() + self.assertAlmostEqual(duals[m.con1], (1 + a1 / (a2 - a1)), 6) + self.assertAlmostEqual(duals[m.con2], -a1 / (a2 - a1), 6) @mark_parameterized.expand(input=_load_tests(all_solvers)) def test_add_and_remove_vars( @@ -1688,8 +1727,8 @@ def test_log(self, name: str, opt_class: Type[SolverBase], use_presolve: bool): m.obj = pyo.Objective(expr=m.x**2 + m.y**2) m.c1 = pyo.Constraint(expr=m.y <= pyo.log(m.x)) res = opt.solve(m) - self.assertAlmostEqual(m.x.value, 0.6529186341994245) - self.assertAlmostEqual(m.y.value, -0.42630274815985264) + self.assertAlmostEqual(m.x.value, 0.6529186341994245, 3) + self.assertAlmostEqual(m.y.value, -0.42630274815985264, 3) @mark_parameterized.expand(input=_load_tests(all_solvers)) def test_with_numpy( @@ -1799,24 +1838,25 @@ def test_solution_loader( self.assertNotIn(m.x, primals) self.assertIn(m.y, primals) self.assertAlmostEqual(primals[m.y], 1) - reduced_costs = res.solution_loader.get_reduced_costs() - self.assertIn(m.x, reduced_costs) - self.assertIn(m.y, reduced_costs) - self.assertAlmostEqual(reduced_costs[m.x], 1) - self.assertAlmostEqual(reduced_costs[m.y], 0) - reduced_costs = res.solution_loader.get_reduced_costs([m.y]) - self.assertNotIn(m.x, reduced_costs) - self.assertIn(m.y, reduced_costs) - self.assertAlmostEqual(reduced_costs[m.y], 0) - duals = res.solution_loader.get_duals() - self.assertIn(m.c1, duals) - self.assertIn(m.c2, duals) - self.assertAlmostEqual(duals[m.c1], 1) - self.assertAlmostEqual(duals[m.c2], 0) - duals = res.solution_loader.get_duals([m.c1]) - self.assertNotIn(m.c2, duals) - self.assertIn(m.c1, duals) - self.assertAlmostEqual(duals[m.c1], 1) + if (name, opt_class) in dual_solvers: + reduced_costs = res.solution_loader.get_reduced_costs() + self.assertIn(m.x, reduced_costs) + self.assertIn(m.y, reduced_costs) + self.assertAlmostEqual(reduced_costs[m.x], 1) + self.assertAlmostEqual(reduced_costs[m.y], 0) + reduced_costs = res.solution_loader.get_reduced_costs([m.y]) + self.assertNotIn(m.x, reduced_costs) + self.assertIn(m.y, reduced_costs) + self.assertAlmostEqual(reduced_costs[m.y], 0) + duals = res.solution_loader.get_duals() + self.assertIn(m.c1, duals) + self.assertIn(m.c2, duals) + self.assertAlmostEqual(duals[m.c1], 1) + self.assertAlmostEqual(duals[m.c2], 0) + duals = res.solution_loader.get_duals([m.c1]) + self.assertNotIn(m.c2, duals) + self.assertIn(m.c1, duals) + self.assertAlmostEqual(duals[m.c1], 1) @mark_parameterized.expand(input=_load_tests(all_solvers)) def test_time_limit( @@ -2297,6 +2337,8 @@ def test_scaling(self, name: str, opt_class: Type[SolverBase], use_presolve: boo opt.config.writer_config.linear_presolve = True else: opt.config.writer_config.linear_presolve = False + if (name, opt_class) not in dual_solvers: + check_duals = False m = pyo.ConcreteModel() m.x = pyo.Var() @@ -2377,6 +2419,43 @@ def test_external_function( res = opt.solve(model) self.assertAlmostEqual(pyo.value(model.o), 0.885603194411, 7) + @mark_parameterized.expand(input=_load_tests(sos_solvers)) + def test_sos(self, name: str, opt_class: Type[SolverBase], use_presolve: bool): + opt: SolverBase = opt_class() + if not opt.available(): + raise unittest.SkipTest(f'Solver {opt.name} not available.') + + m = pyo.ConcreteModel() + m.a = pyo.Set(initialize=[0, 1, 2, 3]) + m.x = pyo.Var(m.a, within=pyo.Binary) + m.obj = pyo.Objective( + expr=sum((i + 1) * m.x[i] for i in range(4)), sense=pyo.maximize + ) + m.c = pyo.SOSConstraint(var=m.x, sos=1) + + res = opt.solve(m) + self.assertEqual(res.solution_status, SolutionStatus.optimal) + self.assertAlmostEqual(res.incumbent_objective, 4) + for i in range(3): + self.assertAlmostEqual(m.x[i].value, 0) + self.assertAlmostEqual(m.x[3].value, 1) + + del m.c + res = opt.solve(m) + self.assertEqual(res.solution_status, SolutionStatus.optimal) + self.assertAlmostEqual(res.incumbent_objective, 10) + for i in range(4): + self.assertAlmostEqual(m.x[i].value, 1) + + m.c = pyo.SOSConstraint(var=m.x, sos=2) + res = opt.solve(m) + self.assertEqual(res.solution_status, SolutionStatus.optimal) + self.assertAlmostEqual(res.incumbent_objective, 7) + for i in range(2): + self.assertAlmostEqual(m.x[i].value, 0) + self.assertAlmostEqual(m.x[2].value, 1) + self.assertAlmostEqual(m.x[3].value, 1) + class TestLegacySolverInterface(unittest.TestCase): @mark_parameterized.expand(input=all_solvers) @@ -2394,7 +2473,8 @@ def test_param_updates(self, name: str, opt_class: Type[SolverBase]): m.obj = pyo.Objective(expr=m.y) m.c1 = pyo.Constraint(expr=(0, m.y - m.a1 * m.x - m.b1, None)) m.c2 = pyo.Constraint(expr=(None, -m.y + m.a2 * m.x + m.b2, 0)) - m.dual = pyo.Suffix(direction=pyo.Suffix.IMPORT) + if (name, opt_class) in dual_solvers: + m.dual = pyo.Suffix(direction=pyo.Suffix.IMPORT) params_to_test = [(1, -1, 2, 1), (1, -2, 2, 1), (1, -1, 3, 1)] for a1, a2, b1, b2 in params_to_test: @@ -2406,8 +2486,9 @@ def test_param_updates(self, name: str, opt_class: Type[SolverBase]): pyo.assert_optimal_termination(res) self.assertAlmostEqual(m.x.value, (b2 - b1) / (a1 - a2)) self.assertAlmostEqual(m.y.value, a1 * (b2 - b1) / (a1 - a2) + b1) - self.assertAlmostEqual(m.dual[m.c1], (1 + a1 / (a2 - a1))) - self.assertAlmostEqual(m.dual[m.c2], a1 / (a2 - a1)) + if (name, opt_class) in dual_solvers: + self.assertAlmostEqual(m.dual[m.c1], (1 + a1 / (a2 - a1))) + self.assertAlmostEqual(m.dual[m.c2], a1 / (a2 - a1)) @mark_parameterized.expand(input=all_solvers) def test_load_solutions(self, name: str, opt_class: Type[SolverBase]): @@ -2418,11 +2499,14 @@ def test_load_solutions(self, name: str, opt_class: Type[SolverBase]): m.x = pyo.Var() m.obj = pyo.Objective(expr=m.x) m.c = pyo.Constraint(expr=(-1, m.x, 1)) - m.dual = pyo.Suffix(direction=pyo.Suffix.IMPORT) + if (name, opt_class) in dual_solvers: + m.dual = pyo.Suffix(direction=pyo.Suffix.IMPORT) res = opt.solve(m, load_solutions=False) pyo.assert_optimal_termination(res) self.assertIsNone(m.x.value) - self.assertNotIn(m.c, m.dual) + if (name, opt_class) in dual_solvers: + self.assertNotIn(m.c, m.dual) m.solutions.load_from(res) self.assertAlmostEqual(m.x.value, -1) - self.assertAlmostEqual(m.dual[m.c], 1) + if (name, opt_class) in dual_solvers: + self.assertAlmostEqual(m.dual[m.c], 1)