diff --git a/doc/OnlineDocs/reference/topical/solvers/index.rst b/doc/OnlineDocs/reference/topical/solvers/index.rst index 400032df076..092841c6e52 100644 --- a/doc/OnlineDocs/reference/topical/solvers/index.rst +++ b/doc/OnlineDocs/reference/topical/solvers/index.rst @@ -8,4 +8,5 @@ Solver Interfaces cplex_persistent.rst gurobi_direct.rst gurobi_persistent.rst + xpress.rst xpress_persistent.rst diff --git a/doc/OnlineDocs/reference/topical/solvers/xpress.rst b/doc/OnlineDocs/reference/topical/solvers/xpress.rst new file mode 100644 index 00000000000..01b3a224d20 --- /dev/null +++ b/doc/OnlineDocs/reference/topical/solvers/xpress.rst @@ -0,0 +1,200 @@ +Xpress (New Interface) +====================== + +.. currentmodule:: pyomo.contrib.solver.solvers.xpress + +Pyomo provides two solver interfaces to the FICO Xpress solver: +:class:`XpressDirect` for one-shot solves, and :class:`XpressPersistent` +for workflows that solve a model repeatedly with small modifications +between solves. + +Both interfaces support the complete range of problem classes that Xpress +handles: LP, MIP, QP, MIQP, NLP, MINLP, second-order cone programs, +and SOS Type 1 and 2 constraints. + +Expression Walker +----------------- + +:class:`XpressDirect` uses a custom expression walker that translates the +full Pyomo expression tree (linear, quadratic, or nonlinear) directly +into an equivalent Xpress expression object, avoiding further intermediate +Python transformations, and handing off to the Xpress C library as directly +as possible. Quadratic terms arising from Cartesian-product expansions are +expanded on the C side. The result is a lean, single-path translation +with no additional overhead for more complex expression types. + +:class:`XpressPersistent` takes a slightly different approach. +Pyomo's ``generate_standard_repn`` runs first: it decomposes each +constraint into its linear and quadratic parts and, crucially, provides +symbolic (non-evaluated) coefficients that are used to register the +mutable-parameter update helpers driving the targeted ``chgMCoef`` / +``chgRHS`` / ``chgQRowCoeff`` calls between solves. If a nonlinear +subexpression remains after that decomposition, the same walker handles +it, producing an Xpress nonlinear expression. + +XpressDirect +------------ + +:class:`XpressDirect` builds a fresh Xpress problem from the Pyomo model +on every call to :meth:`~XpressDirect.solve`. Use it for one-shot solves +or exploratory modeling. + +.. code-block:: python + + from pyomo.contrib.solver.solvers.xpress import XpressDirect + import pyomo.environ as pyo + + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, 10)) + m.c = pyo.Constraint(expr=m.x >= 3) + m.obj = pyo.Objective(expr=m.x) + + res = XpressDirect().solve(m) + +XpressPersistent +---------------- + +:class:`XpressPersistent` keeps the Xpress problem in memory between +solves and uses Pyomo's model-change notification framework to apply +only the minimal set of solver API calls required to reflect each change. + +Mutable :class:`~pyomo.environ.Param` components are tracked +automatically. Updating a parameter value before the next +:meth:`~XpressPersistent.solve` call triggers targeted coefficient or +bound updates (``chgMCoef``, ``chgRHS``, ``chgQRowCoeff``) rather than a +full model rebuild. + +.. code-block:: python + + from pyomo.contrib.solver.solvers.xpress import XpressPersistent + import pyomo.environ as pyo + + m = pyo.ConcreteModel() + m.cost = pyo.Param(mutable=True, initialize=2.0) + m.x = pyo.Var(bounds=(0, 10)) + m.c = pyo.Constraint(expr=m.x >= 3) + m.obj = pyo.Objective(expr=m.cost * m.x) + + opt = XpressPersistent() + opt.solve(m) # full build + m.cost.set_value(5.0) + opt.solve(m) # incremental: only the objective coefficient is updated + +Incremental operations +^^^^^^^^^^^^^^^^^^^^^^ + +Between solves the persistent interface supports: + +- **LP/QP coefficient and bound updates** without row removal, using + the Xpress ``chgMCoef`` / ``chgRHS`` / ``chgQRowCoeff`` API. +- **NLP constraint updates** via row removal and re-insertion (required + when the nonlinear structure changes). +- **Variable fixing and unfixing** through bound updates only. No + constraints are removed or rebuilt as a result of fixing; Xpress + folds fixed variables natively during the solve. Fixing all integer + variables in a MINLP therefore reduces to a sequence of bound calls, + after which Xpress can treat the problem as continuous without any + structural modification to the Pyomo model. +- **Structural modifications**: add and remove constraints, variables, + SOS constraints, and sub-blocks. + +Configuration +------------- + +Common configuration options (time limits, thread count, MIP gaps, +symbolic labels, raw solver options, etc.) are documented on +:class:`~pyomo.contrib.solver.common.config.BranchAndBoundConfig`, which +both interfaces accept as keyword arguments to :meth:`~XpressDirect.solve`. + +Xpress-specific options: + +.. list-table:: + :header-rows: 1 + :widths: 25 75 + + * - Option + - Description + * - ``warmstart`` + - Pass variable values as a MIP start hint (default ``True``). + * - ``pool_solutions`` + - Collect multiple MIP solutions during branch-and-bound. + ``N > 0``: keep a rolling window of the last ``N`` solutions + found. + +Any Xpress control name accepted by ``prob.controls.`` can be passed: + +.. code-block:: python + + res = opt.solve(m, + solver_options={ + 'outputlog': 0, # suppress solver output + 'maxnode': 500, # B&B node limit + 'feastol': 1e-8, # primal feasibility tolerance + } + ) + +Results +------- + +Every :meth:`~XpressDirect.solve` call returns a +:class:`~pyomo.contrib.solver.common.results.Results` object: + +.. code-block:: python + + res = opt.solve(m) + print(res.termination_condition) # e.g. convergenceCriteriaSatisfied + print(res.solution_status) # e.g. optimal + print(res.incumbent_objective) # objective value at the best solution + +See :class:`~pyomo.contrib.solver.common.results.Results` for the full set +of attributes. For NLP problems solved via Xpress SLP, +:attr:`~pyomo.contrib.solver.common.results.Results.solution_status` will +be ``feasible`` rather than ``optimal``, reflecting the local convergence +nature of the algorithm. + +Solution Pool +------------- + +:class:`XpressDirect` and :class:`XpressPersistent` can collect multiple +feasible MIP solutions found during branch-and-bound via the +``pool_solutions`` configuration option. Setting ``pool_solutions=N`` +(N > 0) keeps a rolling window of the last ``N`` solutions found: once +the window is full, the oldest entry is evicted on each new solution, +so the pool always contains the N most recently discovered feasible +solutions. + +.. code-block:: python + + res = opt.solve(m, pool_solutions=5) + loader = res.solution_loader + print(loader.get_number_of_solutions()) # up to 6 (incumbent + pool) + + # Load the incumbent (solution 0) into the model + loader.solution(0).load_vars() + + # Inspect pool entry 1 without modifying the model permanently + with loader.solution(1): + loader.load_vars() + print(m.x.value) + # After the with-block the active solution reverts to the incumbent + +NLP and Nonlinear Expressions +------------------------------ + +All standard Pyomo nonlinear operators, trigonometric and hyperbolic +functions, and user-defined Python callback functions +(``pyo.ExternalFunction``) are supported. + +``pyo.floor`` and ``pyo.ceil`` are not currently supported and raise +:class:`~pyomo.contrib.solver.common.util.IncompatibleModelError`. +These operations must be reformulated by introducing an auxiliary integer +variable together with two linear inequality constraints that encode the +floor or ceil relationship. Adding an integer variable to a continuous +NLP produces a MINLP. + +Testing +------- + +The interface ships with a test suite covering LP, MIP, QP, QCP, NLP, +MINLP, SOS, mutable parameter tracking, incremental structural updates, +and the solution pool. diff --git a/pyomo/contrib/solver/plugins.py b/pyomo/contrib/solver/plugins.py index 4fa18dc9694..af68dd9722f 100644 --- a/pyomo/contrib/solver/plugins.py +++ b/pyomo/contrib/solver/plugins.py @@ -10,6 +10,8 @@ from .common.factory import SolverFactory from .solvers.ipopt import Ipopt, LegacyIpoptSolver +from .solvers.xpress.xpress_direct import XpressDirect +from .solvers.xpress.xpress_persistent import XpressPersistent from .solvers.gurobi.gurobi_direct import GurobiDirect from .solvers.gurobi.gurobi_persistent import GurobiPersistent from .solvers.gurobi.gurobi_direct_minlp import GurobiDirectMINLP @@ -58,3 +60,13 @@ def load(): legacy_name='scip_persistent', doc='Persistent interface pyscipopt', )(ScipPersistent) + SolverFactory.register( + name="xpress_direct", + legacy_name="xpress_direct_v2", + doc="Direct interface to Xpress", + )(XpressDirect) + SolverFactory.register( + name="xpress_persistent", + legacy_name="xpress_persistent_v2", + doc="Persistent interface to Xpress", + )(XpressPersistent) diff --git a/pyomo/contrib/solver/solvers/xpress/__init__.py b/pyomo/contrib/solver/solvers/xpress/__init__.py new file mode 100644 index 00000000000..b60aae37ddc --- /dev/null +++ b/pyomo/contrib/solver/solvers/xpress/__init__.py @@ -0,0 +1,11 @@ +# ____________________________________________________________________________________ +# +# 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 pyomo.contrib.solver.solvers.xpress.xpress_direct import XpressDirect +from pyomo.contrib.solver.solvers.xpress.xpress_persistent import XpressPersistent diff --git a/pyomo/contrib/solver/solvers/xpress/xpress_base.py b/pyomo/contrib/solver/solvers/xpress/xpress_base.py new file mode 100644 index 00000000000..ae66f414017 --- /dev/null +++ b/pyomo/contrib/solver/solvers/xpress/xpress_base.py @@ -0,0 +1,786 @@ +# ____________________________________________________________________________________ +# +# 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. +# ____________________________________________________________________________________ + +"""Xpress interface: expression walker, solution loaders, and build helpers.""" + +import datetime +import io +import math +import os +import time +from typing import Any, Iterator, Mapping, Optional, Sequence, cast + +from pyomo.common.collections import ComponentMap +from pyomo.common.dependencies import attempt_import +from pyomo.common.errors import InfeasibleConstraintException +from pyomo.common.tee import capture_output, TeeStream +from pyomo.common.timing import HierarchicalTimer +from pyomo.common.enums import ObjectiveSense +from pyomo.core.staleflag import StaleFlagManager +from pyomo.core.base.block import BlockData +from pyomo.core.base.var import VarData +from pyomo.core.base.constraint import ConstraintData +from pyomo.core.base.external import PythonCallbackFunction +from pyomo.core.base.sos import SOSConstraintData +from pyomo.core.expr.numeric_expr import ( + AbsExpression, + DivisionExpression, + ExpressionBase, + Expr_ifExpression, + ExternalFunctionExpression, + LinearExpression, + MaxExpression, + MinExpression, + MonomialTermExpression, + NegationExpression, + ProductExpression, + PowExpression, + SumExpression, + UnaryFunctionExpression, +) +from pyomo.core.expr.visitor import StreamBasedExpressionVisitor +from pyomo.core.base import Expression as NamedExpressionComponent +from pyomo.core.expr.numvalue import value +from pyomo.common.gc_manager import PauseGC +from pyomo.repn.util import BeforeChildDispatcher + +from pyomo.common.config import ConfigValue, NonNegativeInt +from pyomo.contrib.solver.common.base import SolverBase, Availability +from pyomo.contrib.solver.common.config import BranchAndBoundConfig +from pyomo.contrib.solver.common.results import ( + Results, + SolutionStatus, + TerminationCondition, + get_infeasible_results, +) +from pyomo.contrib.solver.common.solution_loader import SolutionLoader +from pyomo.contrib.solver.common.util import ( + IncompatibleModelError, + NoDualsError, + NoReducedCostsError, + NoSolutionError, + NoFeasibleSolutionError, + NoOptimalSolutionError, +) + +# L: lower bound; B: upper bound. +_BOUND_TYPE_CODES = ['L', 'U'] + +# Element: (is_binary, is_integer) -> var type char as integer +_VAR_TYPE_CODES: dict[tuple[bool, bool], int] = { + (True, True): 66, # 'B': binary + (False, True): 73, # 'I': integer + (False, False): 67, # 'C': continuous +} + +# Xpress dependent maps, filled by _init_xp_maps +_CON_TYPE_MAP: dict = {} # (is_range, is_equality, has_ub) -> xpress constr type +_OBJ_SENSE_MAP: dict = {} # ObjectiveSense -> xpress.ObjSense +_STOP_TYPE_MAP: dict = {} # xpress.StopType -> TerminationCondition +_SOL_STATUS_MAP: dict = {} # xp.SolStatus -> (TerminationCondition, SolutionStatus) +_VAR_XP_TYPE_MAP: dict = {} # (is_binary, is_integer) -> xpress var type +_XP_FUNCTION_MAP: dict = {} # pyomo fn name -> xpress fn object + + +class _ExitHandlerMap(dict): + """Dict with MRO fallback for unregistered expression subclasses.""" + + def __missing__(self, key): + for cls in key.__mro__: + if cls in self: + self[key] = self[cls] # cache for subsequent lookups + return self[cls] + raise IncompatibleModelError( + f"Expression type '{type(key).__name__}' is not supported by Xpress." + ) + + +_EXIT_HANDLERS: _ExitHandlerMap = _ExitHandlerMap() + + +def _exit_unary(visitor: 'XpressExpressionWalker', node, arg) -> Any: + fn = _XP_FUNCTION_MAP.get(node.getname()) + if fn is None: + raise IncompatibleModelError( + f"Unsupported function '{node.getname()}' in expression. " + "Xpress does not support this function natively." + ) + return fn(arg) + + +def _exit_named_expression(visitor: 'XpressExpressionWalker', node, arg) -> Any: + visitor.subexpression_cache[id(node)] = arg + return arg + + +def _exit_external_function(visitor: 'XpressExpressionWalker', node, *data) -> Any: + """Handle ExternalFunctionExpression: data is (xp_arg1, ..., xp_argN, fcn_id_int).""" + pyo_fcn = node._fcn + if not isinstance(pyo_fcn, PythonCallbackFunction): + raise IncompatibleModelError( + f"ExternalFunction of type '{type(pyo_fcn).__name__}' is not supported; " + "only PythonCallbackFunction is supported." + ) + xp_args = data[:-1] # strip fcn_id (last element is always a plain int) + fcn_id = data[-1] + # No public accessor for gradient capability -> we test private members + has_grad = pyo_fcn._grad is not None or pyo_fcn._fgh is not None + fgh = 1 if has_grad else 0 + + def xp_cb(*vals): + # evaluate_fgh() needs fcn_id and its derivative appended. + f, g, _ = pyo_fcn.evaluate_fgh([*vals, fcn_id], fixed=None, fgh=fgh) + return (f, *g[:-1]) if g is not None else f + + # Each ExternalFunction must have a unique __name__ in Xpress. + xp_cb.__name__ = f'xp_user_{id(pyo_fcn)}' + return xp.user(xp_cb, *xp_args, derivatives="always" if has_grad else "never") + + +def _init_xpress(xp, xpress_available): + """Populate all Pyomo-Xpress entity maps.""" + if not xpress_available: + return + + _VAR_XP_TYPE_MAP.update( + { + (True, True): xp.binary, + (False, True): xp.integer, + (False, False): xp.continuous, + } + ) + _OBJ_SENSE_MAP.update( + { + ObjectiveSense.minimize: xp.ObjSense.MINIMIZE, + ObjectiveSense.maximize: xp.ObjSense.MAXIMIZE, + } + ) + TC = TerminationCondition + SS = SolutionStatus + _SOL_STATUS_MAP.update( + { + xp.SolStatus.OPTIMAL: (TC.convergenceCriteriaSatisfied, SS.optimal), + xp.SolStatus.FEASIBLE: (TC.convergenceCriteriaSatisfied, SS.feasible), + xp.SolStatus.INFEASIBLE: (TC.provenInfeasible, SS.infeasible), + xp.SolStatus.UNBOUNDED: (TC.unbounded, SS.unknown), + } + ) + _STOP_TYPE_MAP.update( + { + xp.StopType.TIMELIMIT: TC.maxTimeLimit, + xp.StopType.NODELIMIT: TC.iterationLimit, + xp.StopType.ITERLIMIT: TC.iterationLimit, + xp.StopType.WORKLIMIT: TC.iterationLimit, + xp.StopType.MIPGAP: TC.convergenceCriteriaSatisfied, + xp.StopType.CTRLC: TC.interrupted, + xp.StopType.USER: TC.interrupted, + xp.StopType.SOLLIMIT: TC.objectiveLimit, + xp.StopType.GENERICERROR: TC.error, + xp.StopType.MEMORYERROR: TC.error, + xp.StopType.NUMERICALERROR: TC.error, + } + ) + + # (is_range, is_equality, has_ub) -> xp constraint type + _CON_TYPE_MAP.update( + { + (True, False, True): xp.rng, + (False, True, True): xp.eq, + (False, True, False): xp.eq, + (False, False, True): xp.leq, + (False, False, False): xp.geq, + } + ) + + # Pyomo nodes -> Xpress operators map + _EXIT_HANDLERS.update( + { + NegationExpression: lambda v, n, a: -a, + SumExpression: lambda v, n, *a: xp.Sum(list(a)), + ProductExpression: lambda v, n, a, b: a * b, + MonomialTermExpression: lambda v, n, a, b: a * b, + DivisionExpression: lambda v, n, a, b: a / b, + PowExpression: lambda v, n, a, b: a**b, + MaxExpression: lambda v, n, *a: xp.max(list(a)), + MinExpression: lambda v, n, *a: xp.min(list(a)), + AbsExpression: lambda v, n, a: xp.abs(a), + UnaryFunctionExpression: _exit_unary, + NamedExpressionComponent: _exit_named_expression, + ExternalFunctionExpression: _exit_external_function, + } + ) + + # Second level dispatch for unary operators (class + name attribute). + _XP_FUNCTION_MAP.update( + { + 'sin': xp.sin, + 'cos': xp.cos, + 'tan': xp.tan, + 'asin': xp.asin, + 'acos': xp.acos, + 'atan': xp.atan, + 'exp': xp.exp, + 'log': xp.log, + 'log10': xp.log10, + 'sqrt': xp.sqrt, + 'abs': xp.abs, + 'sinh': lambda e: 0.5 * (xp.exp(e) - xp.exp(-e)), + 'cosh': lambda e: 0.5 * (xp.exp(e) + xp.exp(-e)), + 'tanh': lambda e: (xp.exp(e) - xp.exp(-e)) / (xp.exp(e) + xp.exp(-e)), + 'asinh': lambda e: xp.log(e + xp.sqrt(e * e + 1)), + 'acosh': lambda e: xp.log(e + xp.sqrt(e * e - 1)), + 'atanh': lambda e: 0.5 * xp.log((1 + e) / (1 - e)), + # TODO: we could support ceil/floor adding 1 auxiliary int var and constraint + } + ) + + +xp, xpress_available = attempt_import('xpress', callback=_init_xpress) + + +def _register_pool_collector(prob, pool_limit: int) -> 'list[list[float]]': + """Collect MIP solutions; pool[k] holds solution k+1 in discovery order. + + pool_limit > 0: size of the rolling window of the last N solutions found. + """ + pool: list[list[float]] = [] + + def _intsol_cb(cbprob, cbdata): + if len(pool) == pool_limit: + pool.pop(0) # Evict oldest when window is full + pool.append(cbprob.getCallbackSolution()) + + prob.controls.serializepreintsol = 1 # Deterministic solution order + prob.addIntsolCallback(_intsol_cb) + return pool + + +class XpressConfig(BranchAndBoundConfig): + """Configuration options shared by Xpress interfaces.""" + + def __init__( + self, + description=None, + doc=None, + implicit=False, + implicit_domain=None, + visibility=0, + ): + super().__init__( + description=description, + doc=doc, + implicit=implicit, + implicit_domain=implicit_domain, + visibility=visibility, + ) + self.pool_solutions: int = self.declare( + 'pool_solutions', + ConfigValue( + default=0, + domain=NonNegativeInt, + description=( + 'MIP solution pool size (0 = disabled). ' + 'N > 0: keep a rolling window of the last N solutions found.' + ), + ), + ) + self.warmstart: bool = self.declare( + 'warmstart', + ConfigValue( + default=True, + domain=bool, + description='Pass current integer variable values as a MIP warm start.', + ), + ) + + +class EntityMaps: + """Stable handle maps: vars by id(VarData); cons and sos by entity.""" + + def __init__(self, vars: dict, cons: dict, sos: dict): + self.vars = vars # id(VarData) -> xp.var + self.cons = cons # ConstraintData -> xp.constraint + self.sos = sos # SOSConstraintData -> xp.sos + + +class XpressSolutionLoaderBase(SolutionLoader): + """Solution loader for Xpress solvers (direct and persistent).""" + + def __init__( + self, + xp_prob, + pyomo_model: BlockData, + variables: list[VarData], + maps: EntityMaps, + pool_solutions: Optional[list] = None, + ) -> None: + super().__init__() + self._xp_prob = xp_prob + self._pyomo_model = pyomo_model + self._vars = variables + self._maps = maps + self._pool: list = pool_solutions if pool_solutions is not None else [] + self._active_id: int = 0 # 0=incumbent (default), k>0=pool[k-1] + + def _query_vars( + self, variables: Optional[Sequence[VarData]], fn, exc_type: type[Exception] + ) -> Iterator[tuple[VarData, float]]: + """Query variable-valued attribute; yields (VarData, value) pairs.""" + try: + if variables is None: + return zip(self._vars, fn()) + xp_vars = [self._maps.vars[id(var)] for var in variables] + return zip(variables, fn(xp_vars)) + except xp.ModelError as e: + raise exc_type() from e + + def get_number_of_solutions(self) -> int: + return 1 + len(self._pool) + + def get_solution_ids(self) -> list: + return list(range(1 + len(self._pool))) + + def _set_solution_id(self, solution_id: Optional[int]) -> Optional[int]: + prev = self._active_id + self._active_id = solution_id + return prev + + def _get_solution_vals( + self, vars_to_load: Optional[Sequence[VarData]] + ) -> Iterator[tuple[VarData, float]]: + """Yield (VarData, float) pairs for the active solution (0=incumbent, k>0=pool).""" + sid = self._active_id + if not sid: + prob = self._xp_prob + return self._query_vars(vars_to_load, prob.getSolution, NoSolutionError) + k = sid - 1 + if k >= len(self._pool): + raise NoSolutionError( + f'Solution {sid} not available: pool contains {len(self._pool)} solutions.' + ) + sol = self._pool[k] + return self._query_vars( + vars_to_load, + lambda vs=None: sol if vs is None else [sol[v.index] for v in vs], + NoSolutionError, + ) + + def load_vars(self, vars_to_load: Sequence[VarData] | None = None) -> None: + for var, val in self._get_solution_vals(vars_to_load): + var.set_value(val, skip_validation=True) + StaleFlagManager.mark_all_as_stale(delayed=True) + + def get_vars( + self, vars_to_load: Sequence[VarData] | None = None + ) -> Mapping[VarData, float]: + return ComponentMap(self._get_solution_vals(vars_to_load)) + + def get_reduced_costs( + self, vars_to_load: Sequence[VarData] | None = None + ) -> Mapping[VarData, float]: + if self._active_id != 0: + raise NoReducedCostsError( + 'Reduced costs available only for incumbent (solution_id=0).' + ) + prob = self._xp_prob + pairs = self._query_vars(vars_to_load, prob.getRedCosts, NoReducedCostsError) + return ComponentMap(pairs) + + def get_duals( + self, cons_to_load: Sequence[ConstraintData] | None = None + ) -> dict[ConstraintData, float]: + if self._active_id != 0: + raise NoDualsError('Duals available only for incumbent (solution_id=0).') + if cons_to_load is None: + # Explicit keys needed: Xpress cons order != maps.cons order + cons_to_load = list(self._maps.cons.keys()) + xp_cons = list(self._maps.cons.values()) + else: + xp_cons = [self._maps.cons[c] for c in cons_to_load] + try: + vals = self._xp_prob.getDuals(xp_cons) + except xp.ModelError as e: + raise NoDualsError() from e + return {con: float(vals[i]) for i, con in enumerate(cons_to_load)} + + +class XpressSolverMixin(SolverBase): + """Shared solver logic (direct and persistent).""" + + _available = None + _version = None + _xpress_available = xpress_available + + def available(self) -> Availability: + if self._available is None: + if not self._xpress_available: + type(self)._available = Availability.NotFound + else: + try: + xp.problem() + type(self)._available = Availability.FullLicense + except Exception: + type(self)._available = Availability.BadLicense + assert self._available is not None + return self._available + + def version(self) -> tuple: + if not xpress_available: + return tuple() + if XpressSolverMixin._version is None: + XpressSolverMixin._version = tuple( + getattr(xp, 'getVersionNumbers', xp.getversionnumbers)() + ) + return XpressSolverMixin._version + + @staticmethod + def _var_bounds(var: VarData) -> tuple: + """Return variable bounds respecting fixed status.""" + if var.fixed: + if var.value is None: + raise ValueError(f"Variable '{var.name}' is fixed but has no value.") + val = value(var.value) + return val, val + vlb, vub = var.bounds + inf = xp.infinity + return (-inf if vlb is None else value(vlb), inf if vub is None else value(vub)) + + @staticmethod + def _set_var_types(prob, pyo_vars: list[VarData], xp_vars) -> None: + """Set column types in bulk.""" + ctypes = [_VAR_TYPE_CODES[v.is_binary(), v.is_integer()] for v in pyo_vars] + prob.chgColType(xp_vars, ctypes) + + def _set_var_bounds(self, prob, pyo_vars: list[VarData], xp_vars) -> None: + """Set variable bounds in bulk.""" + n = len(pyo_vars) + cbounds = [b for var in pyo_vars for b in self._var_bounds(var)] + cols = [v for v in xp_vars for _ in range(2)] + prob.chgBounds(cols, _BOUND_TYPE_CODES * n, cbounds) + + def _add_vars_impl(self, prob, pyo_vars: list[VarData], symbolic_labels: bool): + """Add columns and set types/bounds. Return xp.var array.""" + n = len(pyo_vars) + if n == 0: + return [] + ncol = prob.attributes.cols + xp_vars = prob.addVariables(len(pyo_vars), name='') + self._set_var_types(prob, pyo_vars, xp_vars) + self._set_var_bounds(prob, pyo_vars, xp_vars) + if symbolic_labels: + names = [v.name for v in pyo_vars] + prob.addNames(xp.Namespaces.COLUMN, names, ncol, ncol + n - 1) + return xp_vars + + @staticmethod + def _add_cons_impl( + prob, + pyo_cons: list[ConstraintData], + walker: 'XpressExpressionWalker', + symbolic_labels: bool, + ) -> list: + """Walk pyomo constraints and build xp.constraint objects.""" + if len(pyo_cons) == 0: + return [] + _walk_expr = walker.walk_expression + xp_cons = [] + with PauseGC(): + for c in pyo_cons: + lb, body, ub = c.to_bounded_expression() + if type(body) is LinearExpression: + result = _before_linear(walker, body)[1] + else: + result = _walk_expr(body) + name = c.name if symbolic_labels else None + vlb, vub = value(lb), value(ub) + xp_cons.append(xp.constraint(body=result, lb=vlb, ub=vub, name=name)) + prob.addConstraint(xp_cons) + return xp_cons + + @staticmethod + def _add_sos_impl( + prob, pyo_sos: list, var_map: dict[int, Any], symbolic_labels: bool = False + ): + """Add SOS sets. Return xp.sos handle array.""" + n = len(pyo_sos) + if n == 0: + return [] + settype: list[int] = [] + setstart: list[int] = [0] + setind: list[int] = [] + refval: list[float] = [] + for con in pyo_sos: + setind.extend(var_map[id(var)].index for var in con.variables) + refval.extend(float(w) for _, w in con.get_items()) + settype.append(ord('1' if con.level == 1 else '2')) + setstart.append(len(setind)) + nsos = prob.attributes.sets + prob.addSets(settype, setstart, setind, refval) + if symbolic_labels: + names = [con.name for con in pyo_sos] + prob.addNames(xp.Namespaces.SET, names, nsos, nsos + n - 1) + return prob.getSOS(first=nsos, last=nsos + n - 1) + + def _warmstart(self, prob, vars: list[VarData], entind: list[int]) -> None: + ws_vals: list[float] = [] + ws_cols: list[int] = [] + for j in entind: + var = vars[j] + if var.value is not None: + ws_vals.append(var.value) + ws_cols.append(j) + if ws_vals: + prob.addMipSol(ws_vals, ws_cols) + + def _apply_solver_controls(self, prob, config: BranchAndBoundConfig) -> None: + if config.time_limit is not None: + prob.controls.timelimit = float(config.time_limit) + if config.threads is not None: + prob.controls.threads = config.threads + if config.rel_gap is not None: + prob.controls.miprelstop = config.rel_gap + if config.abs_gap is not None: + prob.controls.mipabsstop = config.abs_gap + for key, val in config.solver_options.items(): + setattr(prob.controls, key, val) + + def _create_xpress_model( + self, _m: BlockData, _c: BranchAndBoundConfig, _t: HierarchicalTimer + ) -> tuple: + raise NotImplementedError + + def solve(self, model: BlockData, **kwds) -> Results: + start_timestamp = datetime.datetime.now(datetime.timezone.utc) + tick = time.perf_counter() + + config = cast( + BranchAndBoundConfig, self.config(value=kwds, preserve_implicit=True) + ) + if config.timer is None: + config.timer = HierarchicalTimer() + timer = config.timer + + StaleFlagManager.mark_all_as_stale() + log_stream = io.StringIO() + ostreams = [log_stream] + config.tee + + orig_cwd = None + if config.working_dir is not None: + orig_cwd = os.getcwd() + os.chdir(str(config.working_dir)) + + try: + with capture_output(TeeStream(*ostreams), capture_fd=False): + prob, solution_loader, has_obj = self._create_xpress_model( + model, config, timer + ) + self._apply_solver_controls(prob, config) + timer.start('optimize') + prob.optimize() + timer.stop('optimize') + res = self._populate_results(prob, solution_loader, has_obj, config) + + except InfeasibleConstraintException as err: + res = get_infeasible_results( + model=model, + solver=self, + config=config, + err_msg=( + 'The problem was proven infeasible during compilation:\n' f'\t{err}' + ), + ) + finally: + if orig_cwd is not None: + os.chdir(orig_cwd) + + res.solver_log = log_stream.getvalue() + tock = time.perf_counter() + res.timing_info.start_timestamp = start_timestamp + res.timing_info.wall_time = tock - tick + res.timing_info.timer = timer + return res + + def _populate_results( + self, + prob, + solution_loader: XpressSolutionLoaderBase, + has_obj: bool, + config: BranchAndBoundConfig, + ) -> Results: + sv = prob.attributes.solvestatus + ss = prob.attributes.solstatus + st = prob.attributes.stopstatus + + TC = TerminationCondition + SS = SolutionStatus + + if sv == xp.SolveStatus.COMPLETED: + tc, sol_status = _SOL_STATUS_MAP.get(ss, (TC.unknown, SS.noSolution)) + elif sv == xp.SolveStatus.STOPPED: + sol_status = SS.feasible if ss == xp.SolStatus.FEASIBLE else SS.noSolution + tc = _STOP_TYPE_MAP.get(st, TC.unknown) + elif sv == xp.SolveStatus.FAILED: + tc = TC.error + sol_status = SS.noSolution + else: # UNSTARTED + tc = TC.unknown + sol_status = SS.noSolution + + results = Results() + results.termination_condition = tc + results.solution_status = sol_status + results.solution_loader = solution_loader + results.solver_name = self.name + results.solver_version = self.version() + results.solver_config = config + + has_solution = sol_status in (SS.optimal, SS.feasible) + + if has_obj and has_solution: + try: + obj_val = float(prob.attributes.objval) + results.incumbent_objective = ( + None if not math.isfinite(obj_val) else obj_val + ) + except (xp.ModelError, AttributeError): + results.incumbent_objective = None + try: + results.objective_bound = float(prob.attributes.bestbound) + except (xp.ModelError, AttributeError): + results.objective_bound = ( + -math.inf if sol_status == SS.optimal else math.inf + ) + else: + results.incumbent_objective = None + results.objective_bound = None + + results.timing_info.xpress_time = prob.attributes.time + results.extra_info.simplex_iterations = prob.attributes.simplexiter + results.extra_info.barrier_iterations = prob.attributes.bariter + results.extra_info.node_count = prob.attributes.nodes + results.extra_info.mip_solutions_found = prob.attributes.mipsols + + if ( + tc != TC.convergenceCriteriaSatisfied + and config.raise_exception_on_nonoptimal_result + ): + raise NoOptimalSolutionError() + + if config.load_solutions: + if has_solution: + solution_loader.load_solution() + else: + raise NoFeasibleSolutionError() + + return results + + +# ---- Expression Walker ------------------------------------------------------ + + +def _register_variable(visitor: 'XpressExpressionWalker', pyo_var: VarData): + """Register (or get) pyo_var in prob. Return its xp.var.""" + xp_var = visitor.var_map.get(id(pyo_var)) + if xp_var is not None: + return xp_var + lb, ub = XpressSolverMixin._var_bounds(pyo_var) + vtype = _VAR_XP_TYPE_MAP[pyo_var.is_binary(), pyo_var.is_integer()] + vname = pyo_var.name if visitor.use_names else None + xp_var = visitor.prob.addVariable(lb=lb, ub=ub, vartype=vtype, name=vname) + visitor.var_map[id(pyo_var)] = xp_var + if visitor.registered_vars is not None: + visitor.registered_vars.append(pyo_var) + return xp_var + + +def _before_monomial(visitor: 'XpressExpressionWalker', child: MonomialTermExpression): + coef, var = child.args + coef_val = value(coef) + xp_var = _register_variable(visitor, var) + return False, coef_val * xp_var + + +def _before_linear(visitor: 'XpressExpressionWalker', child: LinearExpression): + terms = [] + for arg in child.args: + if isinstance(arg, VarData): + terms.append(_register_variable(visitor, arg)) + elif type(arg) is MonomialTermExpression: + coef, var = arg.args + terms.append(value(coef) * _register_variable(visitor, var)) + else: + terms.append(value(arg)) + return False, xp.Sum(terms) + + +def _before_incompatible(_v, child: ExpressionBase): + raise IncompatibleModelError( + f"Expression '{child}' of type '{type(child).__name__}' " + "is not supported by the Xpress solver." + ) + + +class XpressBeforeChildDispatcher(BeforeChildDispatcher): + __slots__ = () + + def __init__(self) -> None: + super().__init__() + self[MonomialTermExpression] = _before_monomial + self[LinearExpression] = _before_linear + self[Expr_ifExpression] = _before_incompatible + + @staticmethod + def _before_var(visitor: 'XpressExpressionWalker', child): + return False, _register_variable(visitor, child) + + @staticmethod + def _before_named_expression(visitor: 'XpressExpressionWalker', child) -> tuple: + cached = visitor.subexpression_cache.get(id(child), None) + return cached is None, cached + + @staticmethod + def _before_leaf(visitor, child): + return False, value(child) + + _before_npv = _before_leaf + _before_param = _before_leaf + _before_native_numeric = _before_leaf + _before_native_logical = _before_leaf + _before_complex = _before_incompatible + _before_string = _before_incompatible + _before_invalid = _before_incompatible + + +class XpressExpressionWalker(StreamBasedExpressionVisitor): + """Pyomo expression tree -> Xpress expression objects.""" + + before_child_dispatcher = XpressBeforeChildDispatcher() + + def __init__( + self, + var_map: dict[int, Any], + prob, + use_names: bool = False, + registered_vars: 'Optional[list[VarData]]' = None, + ) -> None: + super().__init__() + self.var_map = var_map + self.prob = prob + self.use_names = use_names + self.registered_vars = registered_vars + self.subexpression_cache: dict = {} + + def initializeWalker(self, expr) -> tuple: + return self.beforeChild(None, expr, 0) + + def beforeChild(self, _n, child, _c: int) -> tuple: + return self.before_child_dispatcher[type(child)](self, child) + + def exitNode(self, node: ExpressionBase, data: list) -> Any: + return _EXIT_HANDLERS[type(node)](self, node, *data) diff --git a/pyomo/contrib/solver/solvers/xpress/xpress_direct.py b/pyomo/contrib/solver/solvers/xpress/xpress_direct.py new file mode 100644 index 00000000000..5d31e2422b7 --- /dev/null +++ b/pyomo/contrib/solver/solvers/xpress/xpress_direct.py @@ -0,0 +1,94 @@ +# ____________________________________________________________________________________ +# +# 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 typing import Any + +from pyomo.common.timing import HierarchicalTimer +from pyomo.core.base.block import BlockData +from pyomo.core.base.constraint import Constraint +from pyomo.core.base.objective import Objective +from pyomo.core.base.sos import SOSConstraint +from pyomo.core.base.var import VarData + +from pyomo.contrib.solver.common.config import BranchAndBoundConfig +from pyomo.contrib.solver.common.util import IncompatibleModelError +from .xpress_base import ( + XpressSolverMixin, + XpressSolutionLoaderBase, + EntityMaps, + XpressExpressionWalker, + XpressConfig, + _OBJ_SENSE_MAP, + _register_variable, + _register_pool_collector, + xp, +) + + +class XpressDirectSolutionLoader(XpressSolutionLoaderBase): + """Solution loader for the non-persistent XpressDirect solver.""" + + +class XpressDirect(XpressSolverMixin): + CONFIG = XpressConfig() + + def _create_xpress_model( + self, model: BlockData, config: BranchAndBoundConfig, timer: HierarchicalTimer + ) -> tuple[Any, XpressDirectSolutionLoader, bool]: + timer.start('compile_model') + pyo_objs = list(model.component_data_objects(Objective, active=True)) + pyo_cons = list(model.component_data_objects(Constraint, active=True)) + pyo_sos = list(model.component_data_objects(SOSConstraint, active=True)) + if len(pyo_objs) > 1: + raise IncompatibleModelError( + f'Xpress supports at most one objective (received {len(pyo_objs)}).' + ) + timer.stop('compile_model') + + timer.start('load_model') + xp_prob = xp.problem() + use_names = config.symbolic_solver_labels + var_map: dict[int, Any] = {} + pyo_vars: list[VarData] = [] + walker = XpressExpressionWalker( + var_map, prob=xp_prob, use_names=use_names, registered_vars=pyo_vars + ) + xp_cons = self._add_cons_impl(xp_prob, pyo_cons, walker, use_names) + cons_map = dict(zip(pyo_cons, xp_cons)) + if len(pyo_objs) > 0: + obj_result = walker.walk_expression(pyo_objs[0].expr) + sense = _OBJ_SENSE_MAP[pyo_objs[0].sense] + xp_prob.setObjective(obj_result, sense=sense) + # Register SOS variables if not already in constraints/objective. + for sos_con in pyo_sos: + for var in sos_con.variables: + _register_variable(walker, var) + xp_sos = self._add_sos_impl(xp_prob, pyo_sos, var_map, use_names) + sos_map = dict(zip(pyo_sos, xp_sos)) + timer.stop('load_model') + + maps = EntityMaps(vars=var_map, cons=cons_map, sos=sos_map) + # Warm start only if problem has MIP entities or SOS sets. + if config.warmstart and ( + xp_prob.attributes.mipents > 0 or xp_prob.attributes.sets > 0 + ): + entind = [i for i, var in enumerate(pyo_vars) if not var.is_continuous()] + self._warmstart(xp_prob, pyo_vars, entind) + + pool: list = [] + if config.pool_solutions > 0: + pool = _register_pool_collector(xp_prob, config.pool_solutions) + + return ( + xp_prob, + XpressDirectSolutionLoader( + xp_prob, model, pyo_vars, maps, pool_solutions=pool + ), + len(pyo_objs) > 0, + ) diff --git a/pyomo/contrib/solver/solvers/xpress/xpress_persistent.py b/pyomo/contrib/solver/solvers/xpress/xpress_persistent.py new file mode 100644 index 00000000000..cd59480b3a0 --- /dev/null +++ b/pyomo/contrib/solver/solvers/xpress/xpress_persistent.py @@ -0,0 +1,922 @@ +# ____________________________________________________________________________________ +# +# 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 typing import Any, Collection, Mapping, Optional, Sequence + +from pyomo.common.timing import HierarchicalTimer +from pyomo.core.base.block import BlockData +from pyomo.core.base.constraint import Constraint, ConstraintData +from pyomo.core.base.objective import ObjectiveData +from pyomo.core.base.param import ParamData +from pyomo.core.base.sos import SOSConstraint, SOSConstraintData +from pyomo.core.base.var import VarData +from pyomo.core.expr.numvalue import value +from pyomo.common.collections import ComponentMap +from pyomo.repn import generate_standard_repn + +from pyomo.contrib.solver.common.base import PersistentSolverBase +from pyomo.contrib.solver.common.config import BranchAndBoundConfig +from pyomo.contrib.solver.common.util import ( + IncompatibleModelError, + NoSolutionError, + get_objective, +) +from pyomo.contrib.observer.model_observer import ( + AutoUpdateConfig, + ModelChangeDetector, + Observer, + Reason, +) + +from .xpress_base import ( + XpressSolverMixin, + XpressSolutionLoaderBase, + EntityMaps, + XpressExpressionWalker, + XpressConfig, + _OBJ_SENSE_MAP, + _CON_TYPE_MAP, + _register_pool_collector, + xp, +) + + +def _is_constant(expr) -> bool: + """True if expr has no mutable params (safe to treat as a fixed coefficient).""" + return getattr(expr, 'is_constant', lambda: True)() + + +def _collect_fixed(vars: list) -> ComponentMap: + """Collect and unfix all fixed variables; return ComponentMap for re-fixing.""" + fixed = ComponentMap() + for var in vars: + if var.is_fixed(): + fixed[var] = var.value + var.unfix() + return fixed + + +def _refix(fixed: ComponentMap) -> None: + """Re-fix variables that were temporarily unfixed before repn computation.""" + for var, val in fixed.items(): + var.fix(val) + + +# --------------------------------------------------------------------------- +# Mutable parameter helpers +# --------------------------------------------------------------------------- + + +class _UpdateBatch: + """Accumulates constraint updates. + + Applies LP/QP chgMCoef/chgRHS first, then NL delConstraint+addConstraint. + """ + + def __init__(self): + # chgMCoef(rows, cols, vals) + self.coef_rows = [] + self.coef_cols = [] + self.coef_vals = [] + # chgRHS(rows, vals) + self.rhs_rows = [] + self.rhs_vals = [] + # chgRHSRange(rows, vals) + self.rng_rows = [] + self.rng_vals = [] + # chgQRowCoeff(row, col1, col2, val) + self.quad_updates = [] + # delConstraint(old_cons), addConstraint(new_cons) + self.nl_old_cons = [] + self.nl_new_cons = [] + + def flush(self, prob) -> None: + if self.coef_rows: + prob.chgMCoef(self.coef_rows, self.coef_cols, self.coef_vals) + if self.rhs_rows: + prob.chgRHS(self.rhs_rows, self.rhs_vals) + if self.rng_rows: + prob.chgRHSRange(self.rng_rows, self.rng_vals) + for row, c1, c2, val in self.quad_updates: + prob.chgQRowCoeff(row, c1, c2, val) + if self.nl_old_cons: + prob.delConstraint(self.nl_old_cons) + prob.addConstraint(self.nl_new_cons) + + +class _MutableConstraint: + """Updates mutable params in constraint rows. + + LP/QP: targeted chgMCoef/chgRHS/chgQRowCoeff. + NL: partial rebuild via re-walk of nl_expr. + """ + + __slots__ = ( + '_con', + '_xp_con', + '_lin_vars', + '_lin_coefs', + '_quad_v1s', + '_quad_v2s', + '_quad_coefs', + '_rhs_expr', + '_rng_expr', + '_type', + '_stable_xp', + '_nl_expr', + ) + + def __init__( + self, + con, + xp_con, + lin_vars: list, + lin_coefs: list, + quad_v1s: list, + quad_v2s: list, + quad_coefs: list, + rhs_expr: Any, + rng_expr: Any, + con_type, + stable_xp: Any, + nl_expr: Any, + ) -> None: + self._con = con + self._xp_con = xp_con + self._lin_vars = lin_vars + self._lin_coefs = lin_coefs + self._quad_v1s = quad_v1s + self._quad_v2s = quad_v2s + self._quad_coefs = quad_coefs + self._rhs_expr = rhs_expr + self._rng_expr = rng_expr + self._type = con_type + self._stable_xp = stable_xp + self._nl_expr = nl_expr + + def collect(self, batch: _UpdateBatch, walker, use_names: bool, maps) -> None: + if self._nl_expr is not None: + batch.nl_old_cons.append(self._xp_con) + self._xp_con = self._make_xp_con(walker, use_names) + maps.cons[self._con] = self._xp_con + batch.nl_new_cons.append(self._xp_con) + return + row = self._xp_con + if self._lin_vars: + batch.coef_rows.extend([row] * len(self._lin_vars)) + batch.coef_cols.extend(self._lin_vars) + batch.coef_vals.extend(value(c) for c in self._lin_coefs) + if self._rhs_expr is not None: + batch.rhs_rows.append(row) + batch.rhs_vals.append(value(self._rhs_expr)) + if self._rng_expr is not None: + batch.rng_rows.append(row) + batch.rng_vals.append(value(self._rng_expr)) + if self._quad_v1s: + batch.quad_updates.extend( + (row, v1, v2, value(c)) + for v1, v2, c in zip(self._quad_v1s, self._quad_v2s, self._quad_coefs) + ) + + def _make_xp_con(self, walker, use_names: bool) -> Any: + nl_expr = ( + walker.walk_expression(self._nl_expr) if self._nl_expr is not None else None + ) + return xp.constraint( + body=_assemble_xp_expr( + stable_xp=self._stable_xp, + lin_vars=self._lin_vars, + lin_coefs=self._lin_coefs, + quad_v1s=self._quad_v1s, + quad_v2s=self._quad_v2s, + quad_coefs=self._quad_coefs, + constant=None, + nl_expr=nl_expr, + ), + rhs=value(self._rhs_expr), + rhsrange=value(self._rng_expr), + type=self._type, + name=self._con.name if use_names else None, + ) + + +def _assemble_xp_expr( + stable_xp, + lin_vars: list, + lin_coefs: list, + quad_v1s: list, + quad_v2s: list, + quad_coefs: list, + constant, + nl_expr, +) -> Any: + """Assemble Xpress expression from precomputed parts.""" + parts = [] + if stable_xp is not None: + parts.append(stable_xp) + parts.extend(value(c) * v for c, v in zip(lin_coefs, lin_vars)) + parts.extend( + value(c) * v1 * v2 for c, v1, v2 in zip(quad_coefs, quad_v1s, quad_v2s) + ) + if constant is not None: + parts.append(value(constant)) + if nl_expr is not None: + parts.append(nl_expr) + return xp.Sum(parts) if len(parts) > 0 else 0.0 + + +class _MutableObjective: + """Updates mutable params in the objective. + + LP/QP: chgObj + chgMQObj. + NL: full setObjective with re-walked nl_expr. + """ + + __slots__ = ( + '_lin_vars', + '_lin_coefs', + '_quad_v1s', + '_quad_v2s', + '_quad_coefs', + '_constant', + '_stable_xp', + '_nl_expr', + ) + + def __init__( + self, + lin_vars: list, + lin_coefs: list, + quad_v1s: list, + quad_v2s: list, + quad_coefs: list, + constant, + stable_xp: Any, + nl_expr: Any, + ) -> None: + self._lin_vars = lin_vars + self._lin_coefs = lin_coefs + self._quad_v1s = quad_v1s + self._quad_v2s = quad_v2s + self._quad_coefs = quad_coefs + self._constant = constant + self._stable_xp = stable_xp + self._nl_expr = nl_expr + + def update(self, prob, walker: 'XpressExpressionWalker') -> None: + if self._nl_expr is None: + if self._constant is not None: + # -1 is the objective index. chgObj negates objective value internally. + prob.chgObj([-1], [-value(self._constant)]) + if len(self._lin_vars) > 0: + vals = [value(c) for c in self._lin_coefs] + prob.chgObj(self._lin_vars, vals) + if len(self._quad_v1s) > 0: + # Xpress QP objective: (1/2)*x'Qx; chgMQObj expects Hessian-scaled values. + # Diagonal entries must be doubled. Off-diagonal appear twice in symmetric Hessian. + prob.chgMQObj( + self._quad_v1s, + self._quad_v2s, + [ + value(c) * (2 if v1 is v2 else 1) + for v1, v2, c in zip( + self._quad_v1s, self._quad_v2s, self._quad_coefs + ) + ], + ) + else: + prob.setObjective( + _assemble_xp_expr( + stable_xp=self._stable_xp, + lin_vars=self._lin_vars, + lin_coefs=self._lin_coefs, + quad_v1s=self._quad_v1s, + quad_v2s=self._quad_v2s, + quad_coefs=self._quad_coefs, + constant=self._constant, + nl_expr=walker.walk_expression(self._nl_expr), + ) + ) + + +class XpressPersistentSolutionLoader(XpressSolutionLoaderBase): + """Solution loader. Invalidated before re-solve.""" + + def __init__( + self, + prob, + pyomo_model: BlockData, + variables: list[VarData], + maps: EntityMaps, + pool_solutions: Optional[list] = None, + ) -> None: + super().__init__(prob, pyomo_model, variables, maps, pool_solutions) + self._valid = True + + def invalidate(self) -> None: + self._valid = False + + def _assert_valid(self) -> None: + if not self._valid: + raise NoSolutionError( + 'The results from the previous solve are no longer valid because ' + 'the model has been modified.' + ) + + def load_vars(self, vars_to_load: Sequence[VarData] | None = None) -> None: + self._assert_valid() + return super().load_vars(vars_to_load) + + def get_vars( + self, vars_to_load: Sequence[VarData] | None = None + ) -> Mapping[VarData, float]: + self._assert_valid() + return super().get_vars(vars_to_load) + + def get_duals( + self, cons_to_load: Sequence[ConstraintData] | None = None + ) -> dict[ConstraintData, float]: + self._assert_valid() + return super().get_duals(cons_to_load) + + def get_reduced_costs( + self, vars_to_load: Sequence[VarData] | None = None + ) -> Mapping[VarData, float]: + self._assert_valid() + return super().get_reduced_costs(vars_to_load) + + +class XpressPersistentConfig(XpressConfig): + """Config for XpressPersistent.""" + + def __init__( + self, + description=None, + doc=None, + implicit=False, + implicit_domain=None, + visibility=0, + ): + XpressConfig.__init__( + self, + description=description, + doc=doc, + implicit=implicit, + implicit_domain=implicit_domain, + visibility=visibility, + ) + self.auto_updates = self.declare('auto_updates', AutoUpdateConfig()) + + +class XpressPersistent(XpressSolverMixin, PersistentSolverBase, Observer): + """Persistent Xpress solver Interface.""" + + CONFIG = XpressPersistentConfig() + + def __init__(self, **kwds): + super().__init__(**kwds) + self._xp_prob = None # xpress.problem() object + self._pyomo_model: Optional[BlockData] = None + self._vars: Optional[list[VarData]] = None + self._maps: Optional[EntityMaps] = None # Pyomo -> Xpress entities maps + self._objective: Optional[ObjectiveData] = None # Active Pyomo objective + self._last_solution_loader: Optional[XpressPersistentSolutionLoader] = None + self._change_detector: Optional[ModelChangeDetector] = None + self._mutable_helpers: dict[ConstraintData, _MutableConstraint] = {} + self._mutable_objective: Optional[_MutableObjective] = None + self._use_names: bool = False # symbolic_solver_labels config + self._walker: XpressExpressionWalker | None = None + + def _clear(self): + self._xp_prob = None + self._pyomo_model = None + self._vars = None + self._maps = None + self._objective = None + self._last_solution_loader = None + self._change_detector = None + self._mutable_helpers = {} + self._mutable_objective = None + self._use_names = False + self._walker = None + + def _invalidate_last_results(self): + if self._last_solution_loader is not None: + self._last_solution_loader.invalidate() + self._last_solution_loader = None + + def _create_xpress_model( + self, model: BlockData, config: BranchAndBoundConfig, timer: HierarchicalTimer + ) -> tuple[Any, XpressPersistentSolutionLoader, bool]: + self._invalidate_last_results() + + if model is self._pyomo_model: + timer.start('update') + self.update(timer=timer, auto_updates=config.auto_updates) + timer.stop('update') + else: + timer.start('set_instance') + self.set_instance( + model, + _t=timer, + use_names=config.symbolic_solver_labels, + auto_updates=config.auto_updates, + ) + timer.stop('set_instance') + + assert self._vars is not None + assert self._maps is not None + has_obj = self._objective is not None + xp_prob = self._xp_prob + vars = self._vars + + pool: list = [] + if config.pool_solutions > 0: + pool = _register_pool_collector(xp_prob, config.pool_solutions) + + self._last_solution_loader = XpressPersistentSolutionLoader( + xp_prob, model, vars, self._maps, pool + ) + + if config.warmstart and ( + xp_prob.attributes.mipents > 0 or xp_prob.attributes.sets > 0 + ): + entind = [i for i, var in enumerate(vars) if not var.is_continuous()] + self._warmstart(self._xp_prob, vars, entind) + + return xp_prob, self._last_solution_loader, has_obj + + def set_instance(self, model, _t=None, use_names=False, auto_updates=None): + self._clear() + self._pyomo_model = model + self._use_names = use_names + self._xp_prob = xp.problem() + self._maps = EntityMaps(vars={}, cons={}, sos={}) + self._vars = [] + + detector_kwds = {} if auto_updates is None else auto_updates + self._change_detector = ModelChangeDetector( + model=model, observers=[self], **detector_kwds + ) + + def update(self, timer=None, auto_updates=None): + if self._pyomo_model is None: + raise RuntimeError('set_instance must be called before update') + assert self._change_detector is not None + detector_kwds = {} if auto_updates is None else auto_updates + self._change_detector.update(timer=timer, **detector_kwds) + + def _add_variables(self, pyo_vars: list[VarData]) -> None: + assert self._maps is not None + assert self._vars is not None + xp_vars = self._add_vars_impl(self._xp_prob, pyo_vars, self._use_names) + self._maps.vars.update((id(pv), xv) for pv, xv in zip(pyo_vars, xp_vars)) + self._vars.extend(pyo_vars) + + def _remove_variables(self, pyo_vars: list[VarData]) -> None: + assert self._maps is not None + assert self._vars is not None + self._xp_prob.delVariable([self._maps.vars.pop(id(var)) for var in pyo_vars]) + removed = {id(var) for var in pyo_vars} + self._vars = [var for var in self._vars if id(var) not in removed] + + def _update_variables(self, variables: Mapping[VarData, Reason]) -> None: + assert self._maps is not None + self._invalidate_last_results() + new_vars = [] + del_vars = [] + mod_vars = [] + for var, reason in variables.items(): + if reason & Reason.added: + new_vars.append(var) + elif reason & Reason.removed: + del_vars.append(var) + else: + mod_vars.append(var) + if len(new_vars) > 0: + self._add_variables(new_vars) + if len(del_vars) > 0: + self._remove_variables(del_vars) + if len(mod_vars) > 0: + xp_vars = [self._maps.vars[id(var)] for var in mod_vars] + self._set_var_bounds(self._xp_prob, mod_vars, xp_vars) + self._set_var_types(self._xp_prob, mod_vars, xp_vars) + + # ----------------------------------------------------------------------- + # repn-based xp expression builder + # ----------------------------------------------------------------------- + + def _build_xp_from_repn(self, repn) -> Any: + """Build Xpress expression from StandardRepn.""" + var_map = self._maps.vars + nl_expr = None + if repn.nonlinear_expr is not None: + if self._walker is None: + self._walker = XpressExpressionWalker( + var_map, self._xp_prob, use_names=self._use_names + ) + nl_expr = self._walker.walk_expression(repn.nonlinear_expr) + + return _assemble_xp_expr( + stable_xp=None, + lin_vars=[var_map[id(v)] for v in repn.linear_vars], + lin_coefs=repn.linear_coefs, + quad_v1s=[var_map[id(v1)] for v1, _ in repn.quadratic_vars], + quad_v2s=[var_map[id(v2)] for _, v2 in repn.quadratic_vars], + quad_coefs=repn.quadratic_coefs, + constant=repn.constant, + nl_expr=nl_expr, + ) + + # ----------------------------------------------------------------------- + # Mutable helper registration + # ----------------------------------------------------------------------- + + def _register_mutable_constraint( + self, con: ConstraintData, xp_con, repn, lb, ub + ) -> None: + has_ub = ub is not None + is_range = has_ub and lb is not None and not con.equality + is_nl = repn.nonlinear_expr is not None + + rhs_ref = ub if has_ub else lb + const_mutable = not _is_constant(repn.constant) + mut_rhs = not _is_constant(rhs_ref) + + # NL always needs rhs/rng for rebuild; LP/QP only if mutable. + if rhs_ref is None: + rhs_expr = None + elif is_nl or const_mutable or mut_rhs: + rhs_expr = rhs_ref - repn.constant + if not is_nl and _is_constant(rhs_expr): + rhs_expr = None + else: + rhs_expr = None + rng_expr = ( + (ub - lb) + if is_range and (is_nl or not _is_constant(ub) or not _is_constant(lb)) + else None + ) + + if is_nl: + stable_xp = 0.0 + con_type = _CON_TYPE_MAP[is_range, con.equality, has_ub] + else: + stable_xp = con_type = None + + lin_vars, lin_coefs = [], [] + var_map = self._maps.vars + for coef, var in zip(repn.linear_coefs, repn.linear_vars): + if not _is_constant(coef): + lin_vars.append(var_map[id(var)]) + lin_coefs.append(coef) + elif is_nl: + stable_xp += value(coef) * var_map[id(var)] + + quad_v1s, quad_v2s, quad_coefs = [], [], [] + for coef, (v1, v2) in zip(repn.quadratic_coefs, repn.quadratic_vars): + if not _is_constant(coef): + quad_v1s.append(var_map[id(v1)]) + quad_v2s.append(var_map[id(v2)]) + quad_coefs.append(coef) + elif is_nl: + stable_xp += value(coef) * var_map[id(v1)] * var_map[id(v2)] + + if ( + is_nl + or len(lin_vars) > 0 + or len(quad_v1s) > 0 + or rhs_expr is not None + or rng_expr is not None + ): + self._mutable_helpers[con] = _MutableConstraint( + con=con, + xp_con=xp_con, + lin_vars=lin_vars, + lin_coefs=lin_coefs, + quad_v1s=quad_v1s, + quad_v2s=quad_v2s, + quad_coefs=quad_coefs, + rhs_expr=rhs_expr, + rng_expr=rng_expr, + con_type=con_type, + stable_xp=stable_xp, + nl_expr=repn.nonlinear_expr, + ) + + # ----------------------------------------------------------------------- + # Constraint management + # ----------------------------------------------------------------------- + + def _add_constraints(self, pyo_cons: list[ConstraintData]) -> None: + assert self._maps is not None + if len(pyo_cons) == 0: + return + xp_cons = [] + fixed_vars = _collect_fixed(self._vars) + try: + for con in pyo_cons: + lb, body, ub = con.to_bounded_expression() + repn = generate_standard_repn(body, compute_values=False) + name = con.name if self._use_names else None + xp_expr = self._build_xp_from_repn(repn) + + vlb, vub = value(lb), value(ub) + xp_con = xp.constraint(body=xp_expr, lb=vlb, ub=vub, name=name) + xp_cons.append(xp_con) + self._maps.cons[con] = xp_con + self._register_mutable_constraint(con, xp_con, repn, lb, ub) + finally: + _refix(fixed_vars) + self._xp_prob.addConstraint(xp_cons) + + def _remove_constraints(self, cons: list[ConstraintData]) -> None: + assert self._maps is not None + self._xp_prob.delConstraint([self._maps.cons.pop(con) for con in cons]) + for con in cons: + self._mutable_helpers.pop(con, None) + + def _add_sos_constraints(self, cons: list[SOSConstraintData]) -> None: + assert self._maps is not None + xp_sos = self._add_sos_impl( + self._xp_prob, cons, self._maps.vars, self._use_names + ) + self._maps.sos.update(zip(cons, xp_sos)) + + def _remove_sos_constraints(self, cons: list[SOSConstraintData]) -> None: + assert self._maps is not None + xp_sos = [self._maps.sos.pop(con) for con in cons] + self._xp_prob.delSOS(xp_sos) + + # ----------------------------------------------------------------------- + # Objective management + # ----------------------------------------------------------------------- + + def _clear_objective(self): + self._xp_prob.delObj(0) + + def _set_objective(self, obj: ObjectiveData | None) -> None: + assert self._maps is not None + self._objective = obj + self._mutable_objective = None + self._clear_objective() + if obj is None: + return + + fixed_vars = _collect_fixed(self._vars) + repn = generate_standard_repn(obj.expr, compute_values=False) + _refix(fixed_vars) + + sense = _OBJ_SENSE_MAP[obj.sense] + self._xp_prob.setObjective(self._build_xp_from_repn(repn), sense=sense) + + is_nl = repn.nonlinear_expr is not None + lin_vars, lin_coefs, stable_lin = [], [], [] + for coef, var in zip(repn.linear_coefs, repn.linear_vars): + xp_var = self._maps.vars[id(var)] + if not _is_constant(coef): + lin_vars.append(xp_var) + lin_coefs.append(coef) + elif is_nl: + stable_lin.append(value(coef) * xp_var) + + quad_v1s, quad_v2s, quad_coefs, stable_quad = [], [], [], [] + for coef, (v1, v2) in zip(repn.quadratic_coefs, repn.quadratic_vars): + xp_v1 = self._maps.vars[id(v1)] + xp_v2 = self._maps.vars[id(v2)] + if not _is_constant(coef): + quad_v1s.append(xp_v1) + quad_v2s.append(xp_v2) + quad_coefs.append(coef) + elif is_nl: + stable_quad.append(value(coef) * xp_v1 * xp_v2) + + stable_xp = ( + xp.Sum(stable_lin + stable_quad) + if is_nl and (stable_lin or stable_quad) + else 0.0 + ) + constant = repn.constant if is_nl or not _is_constant(repn.constant) else None + + if lin_vars or quad_v1s or constant is not None or is_nl: + self._mutable_objective = _MutableObjective( + lin_vars=lin_vars, + lin_coefs=lin_coefs, + quad_v1s=quad_v1s, + quad_v2s=quad_v2s, + quad_coefs=quad_coefs, + constant=constant, + stable_xp=stable_xp if is_nl else None, + nl_expr=repn.nonlinear_expr, + ) + + def _update_constraints(self, cons: Mapping[ConstraintData, Reason]) -> None: + self._invalidate_last_results() + old_cons = [c for c, r in cons.items() if r & (Reason.removed | Reason.expr)] + new_cons = [c for c, r in cons.items() if r & (Reason.added | Reason.expr)] + if len(old_cons) > 0: + self._remove_constraints(old_cons) + if len(new_cons) > 0: + self._add_constraints(new_cons) + + def _update_sos_constraints(self, cons: Mapping[SOSConstraintData, Reason]) -> None: + self._invalidate_last_results() + old_sos = [ + s for s, r in cons.items() if r & (Reason.removed | Reason.sos_items) + ] + new_sos = [s for s, r in cons.items() if r & (Reason.added | Reason.sos_items)] + if len(old_sos) > 0: + self._remove_sos_constraints(old_sos) + if len(new_sos) > 0: + self._add_sos_constraints(new_sos) + + def _update_objectives(self, objs: Mapping[ObjectiveData, Reason]) -> None: + assert self._pyomo_model is not None + self._invalidate_last_results() + any_added = False + for obj, reason in objs.items(): + if reason & (Reason.added | Reason.expr): + self._set_objective(obj) + any_added = True + elif reason & Reason.removed: + if obj is self._objective: + self._set_objective(None) + elif reason & Reason.sense: + self._xp_prob.chgObjSense(_OBJ_SENSE_MAP[obj.sense]) + if any_added: + try: + get_objective(self._pyomo_model) + except ValueError as e: + raise IncompatibleModelError( + 'Xpress supports at most one active objective. ' + 'Deactivate extras with obj.deactivate().' + ) from e + + def _update_parameters(self, params: Mapping[ParamData, Reason]) -> None: + if self._change_detector is None: + return + assert self._maps is not None + self._invalidate_last_results() + cd = self._change_detector + affected_cons = set() + affected_vars = {} + affected_obj = False + for p, reason in params.items(): + if not (reason & Reason.value): # type: ignore[operator] + continue + affected_cons.update(cd.get_constraints_impacted_by_param(p)) + if not affected_obj and len(cd.get_objectives_impacted_by_param(p)) > 0: + affected_obj = True + for var in cd.get_variables_impacted_by_param(p): + affected_vars[id(var)] = var + if len(affected_vars) > 0: + av = list(affected_vars.values()) + self._set_var_bounds( + self._xp_prob, av, [self._maps.vars[id(var)] for var in av] + ) + prob = self._xp_prob + if len(affected_cons) > 0: + batch = _UpdateBatch() + walker, use_names, maps = self._walker, self._use_names, self._maps + for con in affected_cons: + if con in self._mutable_helpers: + self._mutable_helpers[con].collect(batch, walker, use_names, maps) + batch.flush(prob) + if affected_obj and self._objective is not None: + if self._mutable_objective is not None: + self._mutable_objective.update(prob, self._walker) + else: + self._set_objective(self._objective) + + def add_variables(self, variables: list[VarData]) -> None: + assert self._maps is not None + self._add_variables(variables) + + def remove_variables(self, variables: list[VarData]) -> None: + assert self._maps is not None + self._remove_variables(variables) + + def update_variables(self, variables: list[VarData]) -> None: + assert self._change_detector is not None + self._change_detector.update_variables(variables) + + def add_constraints(self, cons: list[ConstraintData]) -> None: + assert self._change_detector is not None + self._change_detector.add_constraints(cons) + + def remove_constraints(self, cons: list[ConstraintData]) -> None: + assert self._change_detector is not None + self._change_detector.remove_constraints(cons) + + def add_sos_constraints(self, cons: list[SOSConstraintData]) -> None: + assert self._change_detector is not None + self._change_detector.add_sos_constraints(cons) + + def remove_sos_constraints(self, cons: list[SOSConstraintData]) -> None: + assert self._change_detector is not None + self._change_detector.remove_sos_constraints(cons) + + def set_objective(self, obj: ObjectiveData) -> None: + assert self._change_detector is not None + self._change_detector.add_objectives([obj]) + + def update_parameters(self, params: Optional[Collection[ParamData]] = None) -> None: + assert self._change_detector is not None + self._change_detector.update_parameters(params) + + def add_block(self, block: BlockData) -> None: + assert self._change_detector is not None + new_cons = list( + block.component_data_objects(Constraint, descend_into=True, active=True) + ) + new_sos = list( + block.component_data_objects(SOSConstraint, descend_into=True, active=True) + ) + if len(new_cons) > 0: + self._change_detector.add_constraints(new_cons) + if len(new_sos) > 0: + self._change_detector.add_sos_constraints(new_sos) + + def remove_block(self, block: BlockData) -> None: + assert self._change_detector is not None + old_cons = list(block.component_data_objects(Constraint, active=True)) + old_sos = list(block.component_data_objects(SOSConstraint, active=True)) + if len(old_cons) > 0: + self._change_detector.remove_constraints(old_cons) + if len(old_sos) > 0: + self._change_detector.remove_sos_constraints(old_sos) + + def has_instance(self) -> bool: + """True if set_instance has been called.""" + return self._pyomo_model is not None + + def get_xpress_problem(self) -> Any: + """Return underlying xp.problem for direct Xpress API access.""" + assert self._xp_prob is not None + return self._xp_prob + + def get_xpress_control(self, *args) -> Any: + assert self._xp_prob is not None + return self._xp_prob.getControl(*args) + + def set_xpress_control(self, *args) -> None: + assert self._xp_prob is not None + self._xp_prob.setControl(*args) + + def get_xpress_attribute(self, *args) -> Any: + assert self._xp_prob is not None + return self._xp_prob.getAttrib(*args) + + def get_xpress_var(self, var: VarData) -> Any: + """Return xp.var handle for a Pyomo variable.""" + assert self._maps is not None + return self._maps.vars[id(var)] + + def get_xpress_constraint(self, con: ConstraintData) -> Any: + """Return xp.constraint handle for a Pyomo constraint.""" + assert self._maps is not None + return self._maps.cons[con] + + def get_xpress_sos(self, con: SOSConstraintData) -> Any: + """Return xp.sos handle for a Pyomo SOS constraint.""" + assert self._maps is not None + return self._maps.sos[con] + + def release(self) -> None: + """Drop xp.problem and release solver resources.""" + self._clear() + + def reset(self) -> None: + """Clear model data from xp.problem, then drop it.""" + if self._xp_prob is not None: + self._xp_prob.reset() + self._clear() + + def write(self, filename: str, flags: str = '') -> None: + """Write loaded Xpress problem to file.""" + self._xp_prob.writeProb(filename, flags) + + def write_iis(self, filename: str) -> str: + """Compute IIS and write to filename in LP format. Must follow infeasible solve.""" + assert self._xp_prob is not None + self._xp_prob.firstIIS(1) + self._xp_prob.writeIIS(1, 0, filename, 'l') + return filename + + def get_iis(self) -> dict: + """Compute IIS and return conflicting Pyomo objects.""" + assert self._xp_prob is not None and self._maps is not None + self._xp_prob.firstIIS(1) + rowind, colind, *_ = self._xp_prob.getIISData(1) + col_to_var = {i: pv for i, pv in enumerate(self._vars)} + row_to_con = {xc.index: pc for pc, xc in self._maps.cons.items()} + return { + 'constraints': [row_to_con[r] for r in rowind if r in row_to_con], + 'variables': [col_to_var[c] for c in colind if c in col_to_var], + } diff --git a/pyomo/contrib/solver/tests/solvers/_xpress_test_utils.py b/pyomo/contrib/solver/tests/solvers/_xpress_test_utils.py new file mode 100644 index 00000000000..196579c3704 --- /dev/null +++ b/pyomo/contrib/solver/tests/solvers/_xpress_test_utils.py @@ -0,0 +1,115 @@ +# ____________________________________________________________________________________ +# +# 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. +# ____________________________________________________________________________________ + +"""Shared test utilities for the Xpress interface test suite.""" + +import pyomo.environ as pyo + +from typing import TypedDict + +from pyomo.contrib.solver.common.results import SolutionStatus, TerminationCondition + + +def _simple_lp(): + m = pyo.ConcreteModel() + m.x = pyo.Var(domain=pyo.NonNegativeReals) + m.y = pyo.Var(domain=pyo.NonNegativeReals) + m.c1 = pyo.Constraint(expr=m.x + m.y <= 4) + m.c2 = pyo.Constraint(expr=2 * m.x + m.y <= 6) + m.obj = pyo.Objective(expr=-m.x - 2 * m.y) + return m + + +def _simple_mip(): + m = pyo.ConcreteModel() + m.x = pyo.Var(domain=pyo.NonNegativeIntegers) + m.y = pyo.Var(domain=pyo.NonNegativeIntegers) + m.c1 = pyo.Constraint(expr=m.x + m.y <= 4) + m.obj = pyo.Objective(expr=-m.x - 2 * m.y) + return m + + +class _SolveExpected(TypedDict, total=False): + termination: TerminationCondition # default: convergenceCriteriaSatisfied + status: SolutionStatus # default: optimal + objective: float # required when status is optimal + vars: list # required when status is optimal; [(pyo_var, float), ...] + obj_places: int # default: 6 + var_places: int # default: 6 + + +def _solve_and_check(test_case, opt, model, expected: _SolveExpected, **solve_kwargs): + """Solve model and verify termination, status, objective, and variables.""" + tc_default = TerminationCondition.convergenceCriteriaSatisfied + st_default = SolutionStatus.optimal + + termination = expected.get('termination', tc_default) + status = expected.get('status', st_default) + obj_places = expected.get('obj_places', 6) + var_places = expected.get('var_places', 6) + + if status == st_default: + assert ( + 'objective' in expected + ), "_solve_and_check: 'objective' is required when status is optimal" + assert ( + 'vars' in expected + ), "_solve_and_check: 'vars' is required when status is optimal" + num_model_vars = model.nvariables() + assert len(expected['vars']) == num_model_vars, ( + f"_solve_and_check: 'vars' must cover all {num_model_vars} active variables " + f"in the model, got {len(expected['vars'])}" + ) + + res = opt.solve(model, **solve_kwargs) + tc = test_case + tc.assertEqual(res.termination_condition, termination) + tc.assertEqual(res.solution_status, status) + if 'objective' in expected: + tc.assertAlmostEqual( + res.incumbent_objective, expected['objective'], places=obj_places + ) + if 'vars' in expected: + for var, expected_val in expected['vars']: + tc.assertAlmostEqual(pyo.value(var), expected_val, places=var_places) + return res + + +def _trivial_model(): + """Minimal model: single bounded variable, no constraints.""" + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, 1)) + m.obj = pyo.Objective(expr=m.x) + return m + + +def _solve_lp_no_load(opt): + """Solve _simple_lp() with load_solutions=False. Return (model, result).""" + m = _simple_lp() + res = opt.solve(m, load_solutions=False) + return m, res + + +def _solve_check_mutate_check( + test_case, + opt, + model, + expected_before: _SolveExpected, + param, + new_value, + expected_after: _SolveExpected, + **solve_kwargs, +): + """Solve-check, mutate param, solve-check. Return solve result before and after.""" + res_before = _solve_and_check( + test_case, opt, model, expected_before, **solve_kwargs + ) + param.set_value(new_value) + res_after = _solve_and_check(test_case, opt, model, expected_after, **solve_kwargs) + return res_before, res_after diff --git a/pyomo/contrib/solver/tests/solvers/test_solvers.py b/pyomo/contrib/solver/tests/solvers/test_solvers.py index 0b30a6eb923..705b04e668d 100644 --- a/pyomo/contrib/solver/tests/solvers/test_solvers.py +++ b/pyomo/contrib/solver/tests/solvers/test_solvers.py @@ -44,6 +44,7 @@ from pyomo.contrib.solver.solvers.ipopt import Ipopt from pyomo.contrib.solver.solvers.knitro.direct import KnitroDirectSolver +from pyomo.contrib.solver.solvers.xpress import XpressDirect, XpressPersistent from pyomo.contrib.solver.tests.solvers import instances from pyomo.core.expr.compare import assertExpressionsEqual from pyomo.core.expr.numeric_expr import LinearExpression @@ -80,6 +81,8 @@ def param_as_standalone_func(cls, p, func, name): ('scip_persistent', ScipPersistent), ('gams', GAMS), ('knitro_direct', KnitroDirectSolver), + ('xpress_direct', XpressDirect), + ('xpress_persistent', XpressPersistent), ] mip_solvers = [ ('gurobi_persistent', GurobiPersistent), @@ -89,6 +92,8 @@ def param_as_standalone_func(cls, p, func, name): ('scip_direct', ScipDirect), ('scip_persistent', ScipPersistent), ('knitro_direct', KnitroDirectSolver), + ('xpress_direct', XpressDirect), + ('xpress_persistent', XpressPersistent), ] nlp_solvers = [ ('gurobi_direct_minlp', GurobiDirectMINLP), @@ -96,6 +101,8 @@ def param_as_standalone_func(cls, p, func, name): ('scip_direct', ScipDirect), ('scip_persistent', ScipPersistent), ('knitro_direct', KnitroDirectSolver), + ('xpress_direct', XpressDirect), + ('xpress_persistent', XpressPersistent), ] qcp_solvers = [ ('gurobi_persistent', GurobiPersistent), @@ -104,6 +111,8 @@ def param_as_standalone_func(cls, p, func, name): ('scip_direct', ScipDirect), ('scip_persistent', ScipPersistent), ('knitro_direct', KnitroDirectSolver), + ('xpress_direct', XpressDirect), + ('xpress_persistent', XpressPersistent), ] qp_solvers = qcp_solvers + [("highs", Highs)] miqcqp_solvers = [ diff --git a/pyomo/contrib/solver/tests/solvers/test_xpress_direct.py b/pyomo/contrib/solver/tests/solvers/test_xpress_direct.py new file mode 100644 index 00000000000..5eb6ec509ee --- /dev/null +++ b/pyomo/contrib/solver/tests/solvers/test_xpress_direct.py @@ -0,0 +1,674 @@ +# ____________________________________________________________________________________ +# +# 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 math +import os +import tempfile +import types + +import pyomo.environ as pyo +import pyomo.common.unittest as unittest +from pyomo.common.timing import HierarchicalTimer + +from pyomo.contrib.solver.common.results import TerminationCondition +from pyomo.contrib.solver.common.util import ( + IncompatibleModelError, + NoDualsError, + NoFeasibleSolutionError, + NoReducedCostsError, + NoSolutionError, +) +from pyomo.contrib.solver.common.results import SolutionStatus +from pyomo.contrib.solver.solvers.xpress import XpressDirect +from pyomo.contrib.solver.solvers.xpress.xpress_base import _exit_external_function, xp +from pyomo.contrib.solver.tests.solvers._xpress_test_utils import ( + _simple_lp, + _simple_mip, + _solve_and_check, + _solve_lp_no_load, +) + +if not XpressDirect().available(): + raise unittest.SkipTest('Xpress not available') + + +def _infeasible(): + m = pyo.ConcreteModel() + m.x = pyo.Var() + m.c1 = pyo.Constraint(expr=m.x >= 10) + m.c2 = pyo.Constraint(expr=m.x <= 1) + m.obj = pyo.Objective(expr=m.x) + return m + + +@unittest.pytest.mark.solver("xpress_direct") +class TestXpressDirect(unittest.TestCase): + def setUp(self): + self.opt = XpressDirect() + + def test_symbolic_solver_labels_lp(self): + m = pyo.ConcreteModel() + m.distinctive_x = pyo.Var(domain=pyo.NonNegativeReals) + m.distinctive_y = pyo.Var(domain=pyo.NonNegativeReals) + m.distinctive_c1 = pyo.Constraint(expr=m.distinctive_x + m.distinctive_y <= 4) + m.distinctive_c2 = pyo.Constraint( + expr=2 * m.distinctive_x + m.distinctive_y <= 6 + ) + m.obj = pyo.Objective(expr=-m.distinctive_x - 2 * m.distinctive_y) + res = _solve_and_check( + self, + self.opt, + m, + { + 'objective': -8.0, + 'vars': [(m.distinctive_x, 0.0), (m.distinctive_y, 4.0)], + }, + symbolic_solver_labels=True, + ) + with tempfile.TemporaryDirectory() as tmp: + base = os.path.join(tmp, 'm') + res.solution_loader._xp_prob.writeProb(base + '.lp', flags='l') + with open(base + '.lp', 'r') as f: + content = f.read() + self.assertIn('distinctive_x', content) + self.assertIn('distinctive_c1', content) + + def test_symbolic_solver_labels_mip(self): + m = pyo.ConcreteModel() + m.distinctive_x = pyo.Var(domain=pyo.NonNegativeIntegers) + m.distinctive_y = pyo.Var(domain=pyo.NonNegativeIntegers) + m.distinctive_c1 = pyo.Constraint(expr=m.distinctive_x + m.distinctive_y <= 4) + m.obj = pyo.Objective(expr=-m.distinctive_x - 2 * m.distinctive_y) + res = _solve_and_check( + self, + self.opt, + m, + { + 'objective': -8.0, + 'vars': [(m.distinctive_x, 0.0), (m.distinctive_y, 4.0)], + }, + symbolic_solver_labels=True, + ) + with tempfile.TemporaryDirectory() as tmp: + base = os.path.join(tmp, 'm') + res.solution_loader._xp_prob.writeProb(base + '.lp', flags='l') + with open(base + '.lp', 'r') as f: + content = f.read() + self.assertIn('distinctive_x', content) + self.assertIn('distinctive_c1', content) + + def test_positive_time_limit(self): + m = _simple_lp() + _solve_and_check( + self, + self.opt, + m, + {'objective': -8.0, 'vars': [(m.x, 0.0), (m.y, 4.0)]}, + time_limit=60, + ) + + def test_solver_options_passthrough(self): + m = _simple_lp() + _solve_and_check( + self, + self.opt, + m, + {'objective': -8.0, 'vars': [(m.x, 0.0), (m.y, 4.0)]}, + solver_options={'outputlog': 0}, + ) + with self.assertRaises(Exception): + self.opt.solve(m, solver_options={'_invalid_control_xyz': 1}) + + def test_rel_gap(self): + m = _simple_mip() + _solve_and_check( + self, + self.opt, + m, + {'objective': -8.0, 'vars': [(m.x, 0.0), (m.y, 4.0)]}, + rel_gap=0.01, + ) + + def test_threads(self): + m = _simple_lp() + _solve_and_check( + self, + self.opt, + m, + {'objective': -8.0, 'vars': [(m.x, 0.0), (m.y, 4.0)]}, + threads=1, + ) + + def test_mip_no_duals_no_reduced_costs(self): + m = _simple_mip() + res = self.opt.solve(m, load_solutions=False) + self.assertEqual( + res.termination_condition, TerminationCondition.convergenceCriteriaSatisfied + ) + with self.assertRaises(NoDualsError): + res.solution_loader.get_duals() + with self.assertRaises(NoReducedCostsError): + res.solution_loader.get_reduced_costs() + + def test_get_vars_no_solution(self): + m = _infeasible() + res = _solve_and_check( + self, + self.opt, + m, + { + 'termination': TerminationCondition.provenInfeasible, + 'status': SolutionStatus.infeasible, + }, + raise_exception_on_nonoptimal_result=False, + load_solutions=False, + ) + with self.assertRaises(NoSolutionError): + res.solution_loader.get_vars() + + def test_warmstart(self): + m = _simple_mip() + m.x.set_value(0) + m.y.set_value(4) + res = _solve_and_check( + self, self.opt, m, {'objective': -8.0, 'vars': [(m.x, 0.0), (m.y, 4.0)]} + ) + self.assertGreaterEqual(res.extra_info.mip_solutions_found, 1) + + def test_extra_info_and_timing(self): + m = _simple_lp() + res = _solve_and_check( + self, self.opt, m, {'objective': -8.0, 'vars': [(m.x, 0.0), (m.y, 4.0)]} + ) + self.assertGreaterEqual(res.timing_info.xpress_time, 0) + self.assertGreaterEqual(res.extra_info.simplex_iterations, 1) + self.assertGreaterEqual(res.extra_info.barrier_iterations, 0) + self.assertEqual(res.extra_info.node_count, 0) + self.assertEqual(res.extra_info.mip_solutions_found, 0) + + def test_load_solutions_infeasible(self): + m = _infeasible() + with self.assertRaises(NoFeasibleSolutionError): + self.opt.solve( + m, raise_exception_on_nonoptimal_result=False, load_solutions=True + ) + + def test_load_vars_subset(self): + m, res = _solve_lp_no_load(self.opt) + m.x.set_value(99.0) + m.y.set_value(99.0) + res.solution_loader.load_vars([m.y]) + self.assertAlmostEqual(m.y.value, 4.0) + self.assertAlmostEqual(m.x.value, 99.0) + res.solution_loader.load_vars([m.x, m.y]) + self.assertAlmostEqual(m.x.value, 0.0) + self.assertAlmostEqual(m.y.value, 4.0) + + def test_get_vars_subset(self): + m, res = _solve_lp_no_load(self.opt) + result = res.solution_loader.get_vars([m.x]) + self.assertIn(m.x, result) + self.assertNotIn(m.y, result) + self.assertAlmostEqual(result[m.x], 0.0) + result = res.solution_loader.get_vars([m.x, m.y]) + self.assertAlmostEqual(result[m.x], 0.0) + self.assertAlmostEqual(result[m.y], 4.0) + result = res.solution_loader.get_vars([m.y, m.x]) + self.assertAlmostEqual(result[m.x], 0.0) + self.assertAlmostEqual(result[m.y], 4.0) + + def test_get_reduced_costs_subset(self): + m, res = _solve_lp_no_load(self.opt) + result = res.solution_loader.get_reduced_costs([m.y]) + self.assertIn(m.y, result) + self.assertNotIn(m.x, result) + self.assertAlmostEqual(result[m.y], 0.0) + result = res.solution_loader.get_reduced_costs([m.x, m.y]) + self.assertIn(m.x, result) + self.assertIn(m.y, result) + self.assertAlmostEqual(result[m.y], 0.0) + + def test_get_vars_all(self): + m, res = _solve_lp_no_load(self.opt) + result = res.solution_loader.get_vars() + self.assertIn(m.x, result) + self.assertIn(m.y, result) + self.assertAlmostEqual(result[m.x], 0.0) + self.assertAlmostEqual(result[m.y], 4.0) + + def test_multiple_objectives_raises(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, 1)) + m.obj1 = pyo.Objective(expr=m.x) + m.obj2 = pyo.Objective(expr=-m.x) + with self.assertRaises(IncompatibleModelError): + self.opt.solve(m) + + def test_sos1_direct(self): + m = pyo.ConcreteModel() + m.x = pyo.Var([1, 2, 3], domain=pyo.NonNegativeReals, bounds=(0, 1)) + m.sos1 = pyo.SOSConstraint(var=m.x, sos=1, weights={1: 1.0, 2: 2.0, 3: 3.0}) + m.obj = pyo.Objective(expr=m.x[1] + 2 * m.x[2] + 3 * m.x[3], sense=pyo.maximize) + _solve_and_check( + self, + self.opt, + m, + {'objective': 3.0, 'vars': [(m.x[1], 0.0), (m.x[2], 0.0), (m.x[3], 1.0)]}, + ) + + def test_sos1_vars_not_in_objective(self): + m = pyo.ConcreteModel() + m.x = pyo.Var([1, 2, 3], domain=pyo.NonNegativeReals, bounds=(0, 1)) + m.y = pyo.Var(bounds=(0, 10)) + m.sos1 = pyo.SOSConstraint(var=m.x, sos=1, weights={1: 1.0, 2: 2.0, 3: 3.0}) + m.obj = pyo.Objective(expr=m.y) + _solve_and_check( + self, + self.opt, + m, + { + 'objective': 0.0, + 'vars': [(m.x[1], 0.0), (m.x[2], 0.0), (m.x[3], 0.0), (m.y, 0.0)], + }, + ) + + def test_sos1_no_duplicate_columns(self): + m = pyo.ConcreteModel() + m.x = pyo.Var([1, 2, 3], domain=pyo.NonNegativeReals, bounds=(0, 1)) + m.sos1 = pyo.SOSConstraint(var=m.x, sos=1, weights={1: 1.0, 2: 2.0, 3: 3.0}) + m.obj = pyo.Objective(expr=m.x[1] + 2 * m.x[2] + 3 * m.x[3], sense=pyo.maximize) + xp_prob = self.opt._create_xpress_model( + m, + self.opt.config, + __import__( + 'pyomo.common.timing', fromlist=['HierarchicalTimer'] + ).HierarchicalTimer(), + )[0] + self.assertEqual(xp_prob.attributes.cols, 3) + + def test_get_duals_single_constraint(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, 5)) + m.c = pyo.Constraint(expr=m.x >= 2) + m.obj = pyo.Objective(expr=m.x) + res = _solve_and_check( + self, self.opt, m, {'objective': 2.0, 'vars': [(m.x, 2.0)]} + ) + duals = res.solution_loader.get_duals([m.c]) + self.assertIn(m.c, duals) + self.assertIsInstance(duals[m.c], float) + + def test_reduced_costs_value_correctness(self): + m, res = _solve_lp_no_load(self.opt) + rcs = res.solution_loader.get_reduced_costs() + self.assertAlmostEqual(rcs[m.x], 1.0, places=6) + self.assertAlmostEqual(rcs[m.y], 0.0, places=6) + + def test_duals_value_correctness(self): + m, res = _solve_lp_no_load(self.opt) + duals = res.solution_loader.get_duals() + self.assertAlmostEqual(duals[m.c1], -2.0, places=6) + self.assertAlmostEqual(duals[m.c2], 0.0, places=6) + + +@unittest.pytest.mark.solver('xpress_persistent') +class TestXpressDirectQuadratic(unittest.TestCase): + def setUp(self): + self.opt = XpressDirect() + + def test_qp_objective_direct(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(domain=pyo.NonNegativeReals) + m.y = pyo.Var(domain=pyo.NonNegativeReals) + m.c = pyo.Constraint(expr=m.x + m.y >= 1) + m.obj = pyo.Objective(expr=m.x**2 + m.y**2) + _solve_and_check( + self, self.opt, m, {'objective': 0.5, 'vars': [(m.x, 0.5), (m.y, 0.5)]} + ) + + def test_qcp_constraint_direct(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, None)) + m.y = pyo.Var(bounds=(0, None)) + m.qc = pyo.Constraint(expr=m.x**2 + m.y**2 <= 1) + m.obj = pyo.Objective(expr=-(m.x + m.y)) + _solve_and_check( + self, + self.opt, + m, + { + 'objective': -math.sqrt(2), + 'vars': [(m.x, math.sqrt(2) / 2), (m.y, math.sqrt(2) / 2)], + 'obj_places': 5, + 'var_places': 5, + }, + ) + + def test_nl_cubic_constraint_direct(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, 10)) + m.c = pyo.Constraint(expr=m.x**3 >= 1) + m.obj = pyo.Objective(expr=m.x) + _solve_and_check(self, self.opt, m, {'objective': 1.0, 'vars': [(m.x, 1.0)]}) + + +@unittest.pytest.mark.solver('xpress_persistent') +class TestXpressDirectMisc(unittest.TestCase): + def setUp(self): + self.opt = XpressDirect() + + def test_nl_cubic_objective_direct(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, 1)) + m.obj = pyo.Objective(expr=m.x**3) + _solve_and_check(self, self.opt, m, {'objective': 0.0, 'vars': [(m.x, 0.0)]}) + + def test_abs_gap_passthrough(self): + m = _simple_mip() + _solve_and_check( + self, + self.opt, + m, + {'objective': -8.0, 'vars': [(m.x, 0.0), (m.y, 4.0)]}, + abs_gap=0.5, + ) + + def test_working_dir_chdir_and_restore(self): + m = _simple_lp() + original_cwd = os.getcwd() + with tempfile.TemporaryDirectory() as tmp: + self.opt.solve(m, working_dir=tmp) + self.assertEqual(os.getcwd(), original_cwd) + + def test_empty_constraint_model(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(2, 10)) + m.obj = pyo.Objective(expr=m.x) + _solve_and_check(self, self.opt, m, {'objective': 2.0, 'vars': [(m.x, 2.0)]}) + + def test_no_objective_feasibility(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, 10)) + m.c = pyo.Constraint(expr=m.x >= 3) + res = self.opt.solve(m) + self.assertEqual( + res.termination_condition, TerminationCondition.convergenceCriteriaSatisfied + ) + self.assertIsNone(res.incumbent_objective) + + def test_constant_objective(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, 10)) + m.c = pyo.Constraint(expr=m.x >= 1) + m.obj = pyo.Objective(expr=5.0) + _solve_and_check(self, self.opt, m, {'objective': 5.0, 'vars': [(m.x, 1.0)]}) + + def test_range_constraint_lp(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(domain=pyo.NonNegativeReals) + m.y = pyo.Var(domain=pyo.NonNegativeReals) + m.c = pyo.Constraint(expr=pyo.inequality(1, m.x + m.y, 3)) + m.obj = pyo.Objective(expr=-2 * m.x - m.y) + _solve_and_check( + self, self.opt, m, {'objective': -6.0, 'vars': [(m.x, 3.0), (m.y, 0.0)]} + ) + + def test_get_duals_no_args(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, 10)) + m.c1 = pyo.Constraint(expr=m.x >= 1) + m.c2 = pyo.Constraint(expr=m.x >= 2) + m.obj = pyo.Objective(expr=m.x) + res = _solve_and_check( + self, self.opt, m, {'objective': 2.0, 'vars': [(m.x, 2.0)]} + ) + duals = res.solution_loader.get_duals() + self.assertIn(m.c1, duals) + self.assertIn(m.c2, duals) + + def test_fixed_var_without_value_raises(self): + m = pyo.ConcreteModel() + m.x = pyo.Var() + m.x.fix() + m.obj = pyo.Objective(expr=m.x) + self.assertIsNone(m.x.value) + with self.assertRaises(ValueError): + self.opt.solve(m) + + def test_controls_unit(self): + m = _simple_lp() + opt = XpressDirect() + xp_prob, _, _ = opt._create_xpress_model(m, opt.config, HierarchicalTimer()) + config = opt.config({'time_limit': 42, 'threads': 2}) + opt._apply_solver_controls(xp_prob, config) + self.assertEqual(xp_prob.controls.timelimit, 42.0) + self.assertEqual(xp_prob.controls.threads, 2) + + def test_infeasible_model_returns_infeasible_result(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(5, 3)) + m.obj = pyo.Objective(expr=m.x) + _solve_and_check( + self, + self.opt, + m, + { + 'termination': TerminationCondition.provenInfeasible, + 'status': SolutionStatus.infeasible, + }, + raise_exception_on_nonoptimal_result=False, + load_solutions=False, + ) + + def test_time_limit_zero(self): + m = _simple_lp() + opt = XpressDirect() + xp_prob, _, _ = opt._create_xpress_model(m, opt.config, HierarchicalTimer()) + config = opt.config({'time_limit': 0}) + opt._apply_solver_controls(xp_prob, config) + self.assertEqual(xp_prob.controls.timelimit, 0.0) + + def test_working_dir_restored_on_exception(self): + m = _simple_lp() + original_cwd = os.getcwd() + with tempfile.TemporaryDirectory() as tmp: + try: + self.opt.solve( + m, working_dir=tmp, solver_options={'_invalid_control_xyz': 1} + ) + except Exception: + pass + self.assertEqual(os.getcwd(), original_cwd) + + +@unittest.pytest.mark.solver('xpress_persistent') +class TestXpressDirectNLP(unittest.TestCase): + def setUp(self): + self.opt = XpressDirect() + + def test_nl_exp_objective_linear_constraints(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, 3)) + m.c = pyo.Constraint(expr=m.x >= 1) + m.obj = pyo.Objective(expr=pyo.exp(m.x)) + _solve_and_check(self, self.opt, m, {'objective': math.e, 'vars': [(m.x, 1.0)]}) + + def test_nl_sin_constraints_linear_objective(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, math.pi / 2)) + m.y = pyo.Var(bounds=(0, 1)) + m.c = pyo.Constraint(expr=pyo.sin(m.x) + m.y <= 1) + m.obj = pyo.Objective(expr=m.x + m.y) + _solve_and_check( + self, self.opt, m, {'objective': 0.0, 'vars': [(m.x, 0.0), (m.y, 0.0)]} + ) + + def test_nl_objective_nl_constraint(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, 2)) + m.c = pyo.Constraint(expr=pyo.exp(m.x) <= 2) + m.obj = pyo.Objective(expr=pyo.sin(m.x)) + _solve_and_check(self, self.opt, m, {'objective': 0.0, 'vars': [(m.x, 0.0)]}) + + def test_nl_range_constraint(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, math.pi)) + m.c = pyo.Constraint(expr=pyo.inequality(0.5, pyo.sin(m.x), 1.0)) + m.obj = pyo.Objective(expr=m.x) + _solve_and_check( + self, self.opt, m, {'objective': math.pi / 6, 'vars': [(m.x, math.pi / 6)]} + ) + self.assertAlmostEqual(pyo.sin(pyo.value(m.x)), 0.5, places=6) + + def test_fixed_variable_in_nl_constraint(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, math.pi)) + m.y = pyo.Var(bounds=(0, 10)) + m.x.fix(math.pi / 2) + m.c = pyo.Constraint(expr=pyo.sin(m.x) + m.y <= 5) + m.obj = pyo.Objective(expr=m.y) + _solve_and_check( + self, + self.opt, + m, + {'objective': 0.0, 'vars': [(m.x, math.pi / 2), (m.y, 0.0)]}, + ) + + def test_nl_abs_objective(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, 5)) + m.obj = pyo.Objective(expr=abs(m.x - 2)) + _solve_and_check(self, self.opt, m, {'objective': 0.0, 'vars': [(m.x, 2.0)]}) + + +@unittest.pytest.mark.solver('xpress_persistent') +@unittest.pytest.mark.solver('xpress_direct') +class TestXpressExternalFunction(unittest.TestCase): + def setUp(self): + self.opt = XpressDirect() + + def test_external_function_no_gradient(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, 5)) + m.y = pyo.Var(bounds=(1, 5)) + m.f = pyo.ExternalFunction(function=lambda x, y: x**2 + y) + m.obj = pyo.Objective(expr=m.f(m.x, m.y)) + _solve_and_check( + self, + self.opt, + m, + { + 'status': SolutionStatus.feasible, + 'objective': 1.0, + 'vars': [(m.x, 0.0), (m.y, 1.0)], + }, + ) + + def test_external_function_with_gradient(self): + def f(x, y): + return x**2 + y + + def grad(args, fixed): + x, _ = args + return [2 * x, 1.0] + + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, 5)) + m.y = pyo.Var(bounds=(1, 5)) + m.f = pyo.ExternalFunction(function=f, gradient=grad) + m.obj = pyo.Objective(expr=m.f(m.x, m.y)) + _solve_and_check( + self, + self.opt, + m, + { + 'status': SolutionStatus.feasible, + 'objective': 1.0, + 'vars': [(m.x, 0.0), (m.y, 1.0)], + }, + ) + + def test_non_supported_external_function_raises(self): + node = types.SimpleNamespace(_fcn=object()) + with self.assertRaises(IncompatibleModelError): + _exit_external_function(None, node) + + def test_external_function_in_constraint(self): + def g(y): + return y + + def grad(args, fixed): + return [1.0] + + m = pyo.ConcreteModel() + m.y = pyo.Var(bounds=(0, 5)) + m.g = pyo.ExternalFunction(function=g, gradient=grad) + m.c = pyo.Constraint(expr=m.g(m.y) >= 1) + m.obj = pyo.Objective(expr=m.y) + _solve_and_check( + self, + self.opt, + m, + {'status': SolutionStatus.feasible, 'objective': 1.0, 'vars': [(m.y, 1.0)]}, + ) + + def test_external_function_multiple(self): + def f1(x): + return x**2 + + def f2(y): + return y + 1.0 + + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, 5)) + m.y = pyo.Var(bounds=(0, 5)) + m.f1 = pyo.ExternalFunction(function=f1) + m.f2 = pyo.ExternalFunction(function=f2) + m.obj = pyo.Objective(expr=m.f1(m.x) + m.f2(m.y)) + if xp.featurequery("Community"): # Free community license + with self.assertRaisesRegex(Exception, r"^\?1152"): + self.opt.solve(m) + else: + _solve_and_check( + self, + self.opt, + m, + { + 'status': SolutionStatus.feasible, + 'objective': 1.0, + 'vars': [(m.x, 0.0), (m.y, 0.0)], + }, + ) + + def test_external_function_fgh_callback(self): + def fgh_func(args, fgh_flag, fixed): + x, y = args + f = x**2 + y + g = [2.0 * x, 1.0] if fgh_flag else None + return f, g, None + + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, 5)) + m.y = pyo.Var(bounds=(1, 5)) + m.f = pyo.ExternalFunction(fgh=fgh_func) + m.obj = pyo.Objective(expr=m.f(m.x, m.y)) + _solve_and_check( + self, + self.opt, + m, + { + 'status': SolutionStatus.feasible, + 'objective': 1.0, + 'vars': [(m.x, 0.0), (m.y, 1.0)], + }, + ) + + +if __name__ == '__main__': + unittest.main() diff --git a/pyomo/contrib/solver/tests/solvers/test_xpress_persistent.py b/pyomo/contrib/solver/tests/solvers/test_xpress_persistent.py new file mode 100644 index 00000000000..5a87b8d7b9b --- /dev/null +++ b/pyomo/contrib/solver/tests/solvers/test_xpress_persistent.py @@ -0,0 +1,1304 @@ +# ____________________________________________________________________________________ +# +# 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 math +import os +import tempfile + +import pyomo.common.unittest as unittest +import pyomo.environ as pyo + +from pyomo.contrib.solver.common.results import TerminationCondition +from pyomo.contrib.solver.common.util import IncompatibleModelError, NoSolutionError +from pyomo.contrib.solver.common.results import SolutionStatus +from pyomo.contrib.solver.solvers.xpress import XpressPersistent +from pyomo.contrib.solver.tests.solvers._xpress_test_utils import ( + _simple_lp, + _simple_mip, + _solve_and_check, + _solve_check_mutate_check, + _trivial_model, +) + +if not XpressPersistent().available(): + raise unittest.SkipTest('Xpress not available') + + +@unittest.pytest.mark.solver('xpress_persistent') +class TestXpressPersistentObjective(unittest.TestCase): + def setUp(self): + self.opt = XpressPersistent() + + def test_remove_objective_between_solves(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(1, 5)) + m.c = pyo.Constraint(expr=m.x >= 2) + m.obj = pyo.Objective(expr=m.x) + + _solve_and_check(self, self.opt, m, {'objective': 2.0, 'vars': [(m.x, 2.0)]}) + + del m.obj + res = self.opt.solve(m) + self.assertEqual(res.solution_status, SolutionStatus.optimal) + self.assertIsNone(res.incumbent_objective) + + def test_active_objective_toggle(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, 5)) + m.obj_min = pyo.Objective(expr=m.x, sense=pyo.minimize) + m.obj_max = pyo.Objective(expr=m.x, sense=pyo.maximize) + m.obj_max.deactivate() + + _solve_and_check(self, self.opt, m, {'objective': 0.0, 'vars': [(m.x, 0.0)]}) + + m.obj_min.deactivate() + m.obj_max.activate() + _solve_and_check(self, self.opt, m, {'objective': 5.0, 'vars': [(m.x, 5.0)]}) + + def test_two_active_objectives_at_set_instance_raises(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, 5)) + m.obj1 = pyo.Objective(expr=m.x, sense=pyo.minimize) + m.obj2 = pyo.Objective(expr=-m.x, sense=pyo.minimize) + with self.assertRaises(IncompatibleModelError): + self.opt.solve(m) + + def test_two_active_objectives_at_update_raises(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, 5)) + m.obj1 = pyo.Objective(expr=m.x, sense=pyo.minimize) + self.opt.solve(m) + + m.obj2 = pyo.Objective(expr=-m.x, sense=pyo.minimize) + with self.assertRaises(IncompatibleModelError): + self.opt.solve(m) + + def test_recovery_after_two_objectives_raises(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, 5)) + m.obj1 = pyo.Objective(expr=m.x, sense=pyo.minimize) + self.opt.solve(m) + + m.obj2 = pyo.Objective(expr=-m.x, sense=pyo.minimize) + with self.assertRaises(IncompatibleModelError): + self.opt.solve(m) + + m.obj2.deactivate() + self.opt.set_instance(m) + _solve_and_check(self, self.opt, m, {'objective': 0.0, 'vars': [(m.x, 0.0)]}) + + +@unittest.pytest.mark.solver('xpress_persistent') +class TestXpressPersistentLifecycle(unittest.TestCase): + def setUp(self): + self.opt = XpressPersistent() + + def test_eager_invalidation_on_mutation(self): + m = _simple_lp() + res = _solve_and_check( + self, self.opt, m, {'objective': -8.0, 'vars': [(m.x, 0.0), (m.y, 4.0)]} + ) + res.solution_loader.get_vars() + m.c3 = pyo.Constraint(expr=m.x + m.y >= 1) + self.opt.add_constraints([m.c3]) + with self.assertRaises(NoSolutionError): + res.solution_loader.get_vars() + + def test_eager_invalidation_on_param_change(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, 10)) + m.p = pyo.Param(mutable=True, initialize=5.0) + m.c = pyo.Constraint(expr=m.x <= m.p) + m.obj = pyo.Objective(expr=-m.x) + res = _solve_and_check( + self, self.opt, m, {'objective': -5.0, 'vars': [(m.x, 5.0)]} + ) + res.solution_loader.get_vars() + m.p.value = 7.0 + self.opt.update_parameters([m.p]) + with self.assertRaises(NoSolutionError): + res.solution_loader.get_vars() + + def test_symbolic_solver_labels_persistent(self): + m = pyo.ConcreteModel() + m.distinctive_var = pyo.Var(domain=pyo.NonNegativeReals) + m.distinctive_con = pyo.Constraint(expr=m.distinctive_var <= 5) + m.obj = pyo.Objective(expr=m.distinctive_var) + + _solve_and_check( + self, + self.opt, + m, + {'objective': 0.0, 'vars': [(m.distinctive_var, 0.0)]}, + symbolic_solver_labels=True, + ) + with tempfile.TemporaryDirectory() as tmp: + base = os.path.join(tmp, 'm') + self.opt.write(base, flags='l') + with open(base + '.lp', 'r') as f: + content = f.read() + self.assertIn('distinctive_var', content) + self.assertIn('distinctive_con', content) + + def test_auto_updates_disable_parameter_tracking(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, 10)) + m.p = pyo.Param(mutable=True, initialize=5.0) + m.c = pyo.Constraint(expr=m.x <= m.p) + m.obj = pyo.Objective(expr=-m.x) + + _solve_and_check(self, self.opt, m, {'objective': -5.0, 'vars': [(m.x, 5.0)]}) + + m.p.value = 7.0 + _solve_and_check( + self, + self.opt, + m, + {'objective': -5.0, 'vars': [(m.x, 5.0)]}, + auto_updates={'update_parameters': False}, + ) + + _solve_and_check(self, self.opt, m, {'objective': -7.0, 'vars': [(m.x, 7.0)]}) + + def test_write_mps_and_lp(self): + m = _simple_lp() + _solve_and_check( + self, self.opt, m, {'objective': -8.0, 'vars': [(m.x, 0.0), (m.y, 4.0)]} + ) + with tempfile.TemporaryDirectory() as tmp: + mps_base = os.path.join(tmp, 'mps_model') + self.opt.write(mps_base) + self.assertTrue(os.path.exists(mps_base + '.mps')) + self.assertGreater(os.path.getsize(mps_base + '.mps'), 0) + + lp_base = os.path.join(tmp, 'lp_model') + self.opt.write(lp_base, flags='l') + self.assertTrue(os.path.exists(lp_base + '.lp')) + self.assertGreater(os.path.getsize(lp_base + '.lp'), 0) + + def test_warmstart_disabled(self): + m = _simple_mip() + m.x.set_value(100) + m.y.set_value(100) + _solve_and_check( + self, + self.opt, + m, + {'objective': -8.0, 'vars': [(m.x, 0.0), (m.y, 4.0)]}, + warmstart=False, + ) + + +@unittest.pytest.mark.solver('xpress_persistent') +class TestXpressPersistentSOS(unittest.TestCase): + def setUp(self): + self.opt = XpressPersistent() + + def test_sos1_initial_and_remove(self): + m = pyo.ConcreteModel() + m.x = pyo.Var([1, 2, 3], domain=pyo.NonNegativeReals, bounds=(0, 1)) + m.sos1 = pyo.SOSConstraint(var=m.x, sos=1, weights={1: 1.0, 2: 2.0, 3: 3.0}) + m.obj = pyo.Objective(expr=m.x[1] + 2 * m.x[2] + 3 * m.x[3], sense=pyo.maximize) + + _solve_and_check( + self, + self.opt, + m, + {'objective': 3.0, 'vars': [(m.x[3], 1.0), (m.x[1], 0.0), (m.x[2], 0.0)]}, + ) + + del m.sos1 + _solve_and_check( + self, + self.opt, + m, + {'objective': 6.0, 'vars': [(m.x[1], 1.0), (m.x[2], 1.0), (m.x[3], 1.0)]}, + ) + + # Public persistent API (explicit call paths) + + def test_add_variables_public_api(self): + m = _trivial_model() + self.opt.set_instance(m) + ncols_before = self.opt._xp_prob.attributes.cols + m.y = pyo.Var(bounds=(0, 1)) + self.opt.add_variables([m.y]) + self.assertGreater(self.opt._xp_prob.attributes.cols, ncols_before) + self.assertIn(id(m.y), self.opt._maps.vars) + + def test_remove_variables_public_api(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, 5)) + m.y = pyo.Var(bounds=(0, 5)) + m.obj = pyo.Objective(expr=m.x + m.y) + self.opt.set_instance(m) + ncols_before = self.opt._xp_prob.attributes.cols + self.opt.remove_variables([m.y]) + self.assertLess(self.opt._xp_prob.attributes.cols, ncols_before) + self.assertNotIn(id(m.y), self.opt._maps.vars) + + def test_update_variables_public_api(self): + m = _trivial_model() + m.c = pyo.Constraint(expr=m.x >= 0.5) + self.opt.set_instance(m) + m.x.setub(3.0) + self.opt.update_variables([m.x]) + _solve_and_check(self, self.opt, m, {'objective': 0.5, 'vars': [(m.x, 0.5)]}) + + def test_set_objective_public_api(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(1, 5)) + m.obj = pyo.Objective(expr=m.x) + self.opt.set_instance(m) + _solve_and_check(self, self.opt, m, {'objective': 1.0, 'vars': [(m.x, 1.0)]}) + m.obj.deactivate() + m.obj2 = pyo.Objective(expr=-m.x) + self.opt.set_objective(m.obj2) + _solve_and_check(self, self.opt, m, {'objective': -5.0, 'vars': [(m.x, 5.0)]}) + + def test_add_remove_sos_constraints_public_api(self): + m = pyo.ConcreteModel() + m.x = pyo.Var([1, 2, 3], domain=pyo.NonNegativeReals, bounds=(0, 1)) + m.obj = pyo.Objective(expr=m.x[1] + 2 * m.x[2] + 3 * m.x[3], sense=pyo.maximize) + self.opt.set_instance(m) + _solve_and_check( + self, + self.opt, + m, + {'objective': 6.0, 'vars': [(m.x[1], 1.0), (m.x[2], 1.0), (m.x[3], 1.0)]}, + ) + m.sos1 = pyo.SOSConstraint(var=m.x, sos=1, weights={1: 1.0, 2: 2.0, 3: 3.0}) + self.opt.add_sos_constraints(list(m.sos1.values())) + _solve_and_check( + self, + self.opt, + m, + {'objective': 3.0, 'vars': [(m.x[3], 1.0), (m.x[1], 0.0), (m.x[2], 0.0)]}, + ) + self.opt.remove_sos_constraints(list(m.sos1.values())) + m.sos1.deactivate() + _solve_and_check( + self, + self.opt, + m, + {'objective': 6.0, 'vars': [(m.x[1], 1.0), (m.x[2], 1.0), (m.x[3], 1.0)]}, + ) + + def test_add_remove_block_public_api(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, 10)) + m.obj = pyo.Objective(expr=m.x) + self.opt.set_instance(m) + _solve_and_check(self, self.opt, m, {'objective': 0.0, 'vars': [(m.x, 0.0)]}) + m.b = pyo.Block() + m.b.c = pyo.Constraint(expr=m.x >= 5) + self.opt.add_block(m.b) + _solve_and_check(self, self.opt, m, {'objective': 5.0, 'vars': [(m.x, 5.0)]}) + self.opt.remove_block(m.b) + m.b.deactivate() + _solve_and_check(self, self.opt, m, {'objective': 0.0, 'vars': [(m.x, 0.0)]}) + + def test_xpress_control_and_attribute(self): + m = _trivial_model() + self.opt.set_instance(m) + self.opt.set_xpress_control('threads', 1) + self.assertEqual(self.opt.get_xpress_control('threads'), 1) + rows = self.opt.get_xpress_attribute('rows') + self.assertGreaterEqual(rows, 0) + + def test_get_xpress_problem_returns_problem(self): + m = _trivial_model() + m.c = pyo.Constraint(expr=m.x >= 0.5) + self.opt.set_instance(m) + prob = self.opt.get_xpress_problem() + self.assertIsNotNone(prob) + xp_con = self.opt.get_xpress_constraint(m.c) + _solve_and_check(self, self.opt, m, {'objective': 0.5, 'vars': [(m.x, 0.5)]}) + slack = prob.getSlacks(xp_con) + self.assertAlmostEqual(slack, 0.0, places=6) + + def test_update_before_set_instance_raises(self): + with self.assertRaises(RuntimeError): + XpressPersistent().update() + + def test_get_xpress_var_returns_handle(self): + m = _trivial_model() + self.opt.set_instance(m) + handle = self.opt.get_xpress_var(m.x) + self.assertIsNotNone(handle) + self.assertGreaterEqual(handle.index, 0) + + def test_get_xpress_constraint_returns_handle(self): + m = _trivial_model() + m.c = pyo.Constraint(expr=m.x >= 0.5) + self.opt.set_instance(m) + handle = self.opt.get_xpress_constraint(m.c) + self.assertIsNotNone(handle) + self.assertGreaterEqual(handle.index, 0) + + def test_get_xpress_sos_returns_handle(self): + m = pyo.ConcreteModel() + m.x = pyo.Var([1, 2], domain=pyo.NonNegativeReals, bounds=(0, 1)) + m.sos = pyo.SOSConstraint(var=m.x, sos=1, weights={1: 1.0, 2: 2.0}) + m.obj = pyo.Objective(expr=m.x[1] + m.x[2]) + self.opt.set_instance(m) + handle = self.opt.get_xpress_sos(list(m.sos.values())[0]) + self.assertIsNotNone(handle) + + def test_release_clears_state(self): + m = _trivial_model() + self.opt.set_instance(m) + self.assertIsNotNone(self.opt._xp_prob) + self.opt.release() + self.assertIsNone(self.opt._xp_prob) + self.assertIsNone(self.opt._maps) + self.assertIsNone(self.opt._change_detector) + self.assertIsNone(self.opt._pyomo_model) + self.assertIsNone(self.opt._vars) + self.assertEqual(self.opt._mutable_helpers, {}) + + def test_reset_clears_state(self): + m = _trivial_model() + self.opt.set_instance(m) + self.assertIsNotNone(self.opt._xp_prob) + self.opt.reset() + self.assertIsNone(self.opt._xp_prob) + self.assertIsNone(self.opt._maps) + self.assertIsNone(self.opt._change_detector) + self.assertIsNone(self.opt._pyomo_model) + self.assertIsNone(self.opt._vars) + self.assertEqual(self.opt._mutable_helpers, {}) + + +@unittest.pytest.mark.solver('xpress_persistent') +class TestXpressPersistentQuadratic(unittest.TestCase): + def setUp(self): + self.opt = XpressPersistent() + + def test_qp_objective_persistent(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(domain=pyo.NonNegativeReals) + m.y = pyo.Var(domain=pyo.NonNegativeReals) + m.c = pyo.Constraint(expr=m.x + m.y >= 1) + m.obj = pyo.Objective(expr=m.x**2 + m.y**2) + _solve_and_check( + self, self.opt, m, {'objective': 0.5, 'vars': [(m.x, 0.5), (m.y, 0.5)]} + ) + _solve_and_check( + self, self.opt, m, {'objective': 0.5, 'vars': [(m.x, 0.5), (m.y, 0.5)]} + ) + + def test_qcp_add_remove_persistent(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, 2)) + m.y = pyo.Var(bounds=(0, 2)) + m.obj = pyo.Objective(expr=-(m.x + m.y)) + self.opt.set_instance(m) + _solve_and_check( + self, self.opt, m, {'objective': -4.0, 'vars': [(m.x, 2.0), (m.y, 2.0)]} + ) + + m.qc = pyo.Constraint(expr=m.x**2 + m.y**2 <= 1) + self.opt.add_constraints([m.qc]) + _solve_and_check( + self, + self.opt, + m, + { + 'objective': -math.sqrt(2), + 'vars': [(m.x, math.sqrt(2) / 2), (m.y, math.sqrt(2) / 2)], + 'obj_places': 5, + 'var_places': 5, + }, + ) + + self.opt.remove_constraints([m.qc]) + m.qc.deactivate() + _solve_and_check( + self, self.opt, m, {'objective': -4.0, 'vars': [(m.x, 2.0), (m.y, 2.0)]} + ) + + def test_mutable_param_in_quadratic_obj(self): + m = pyo.ConcreteModel() + m.p = pyo.Param(mutable=True, initialize=1.0) + m.x = pyo.Var(bounds=(1, None)) + m.obj = pyo.Objective(expr=m.p * m.x**2) + _solve_check_mutate_check( + self, + self.opt, + m, + {'objective': 1.0, 'vars': [(m.x, 1.0)]}, + m.p, + 4.0, + {'objective': 4.0, 'vars': [(m.x, 1.0)], 'obj_places': 4, 'var_places': 4}, + ) + + def test_mutable_param_in_quadratic_constraint(self): + m = pyo.ConcreteModel() + m.p = pyo.Param(mutable=True, initialize=1.0) + m.x = pyo.Var(bounds=(0, None)) + m.y = pyo.Var(bounds=(0, None)) + m.qc = pyo.Constraint(expr=m.x**2 + m.p * m.y**2 <= 1) + m.obj = pyo.Objective(expr=-(m.x + m.y)) + _solve_and_check( + self, + self.opt, + m, + { + 'objective': -math.sqrt(2), + 'vars': [(m.x, math.sqrt(2) / 2), (m.y, math.sqrt(2) / 2)], + 'obj_places': 5, + 'var_places': 5, + }, + ) + m.p.set_value(4.0) + _solve_and_check( + self, + self.opt, + m, + { + 'objective': -(2 / math.sqrt(5) + 1 / (2 * math.sqrt(5))), + 'vars': [(m.x, 2 / math.sqrt(5)), (m.y, 1 / (2 * math.sqrt(5)))], + 'obj_places': 5, + 'var_places': 5, + }, + ) + + def test_mutable_param_in_quadratic_constraint_monomial_form(self): + m = pyo.ConcreteModel() + m.p = pyo.Param(mutable=True, initialize=1.0) + m.x = pyo.Var(bounds=(0, None)) + m.y = pyo.Var(bounds=(0, None)) + m.qc = pyo.Constraint(expr=m.x**2 + (m.p * m.y) * m.y <= 1) + m.obj = pyo.Objective(expr=-(m.x + m.y)) + _solve_and_check( + self, + self.opt, + m, + { + 'objective': -math.sqrt(2), + 'vars': [(m.x, math.sqrt(2) / 2), (m.y, math.sqrt(2) / 2)], + 'obj_places': 5, + 'var_places': 5, + }, + ) + m.p.set_value(4.0) + _solve_and_check( + self, + self.opt, + m, + { + 'objective': -(2 / math.sqrt(5) + 1 / (2 * math.sqrt(5))), + 'vars': [(m.x, 2 / math.sqrt(5)), (m.y, 1 / (2 * math.sqrt(5)))], + 'obj_places': 5, + 'var_places': 5, + }, + ) + + def test_mutable_quadratic_coef_plus_mutable_linear_coef_objective(self): + m = pyo.ConcreteModel() + m.p1 = pyo.Param(mutable=True, initialize=1.0) + m.p2 = pyo.Param(mutable=True, initialize=1.0) + m.p3 = pyo.Param(mutable=True, initialize=4.0) + m.x = pyo.Var() + m.y = pyo.Var() + m.obj = pyo.Objective( + expr=m.p1 * (m.x - 1) ** 2 + m.p2 * (m.y - 6) ** 2 - m.p3 * m.y + ) + m.c = pyo.Constraint(expr=m.x >= m.y) + _solve_and_check( + self, self.opt, m, {'objective': -3.5, 'vars': [(m.x, 4.5), (m.y, 4.5)]} + ) + m.p2.set_value(2.0) + _solve_and_check( + self, self.opt, m, {'objective': -2.0, 'vars': [(m.x, 5.0), (m.y, 5.0)]} + ) + + def test_mutable_quadratic_coef_persistent_analytic(self): + m = pyo.ConcreteModel() + m.p = pyo.Param(mutable=True, initialize=1.0) + m.x = pyo.Var(bounds=(0, None)) + m.y = pyo.Var(bounds=(0, None)) + m.qc = pyo.Constraint(expr=m.x**2 + m.p * m.y**2 <= 1) + m.obj = pyo.Objective(expr=-(m.x + m.y)) + _solve_and_check( + self, + self.opt, + m, + { + 'objective': -math.sqrt(2), + 'vars': [(m.x, math.sqrt(2) / 2), (m.y, math.sqrt(2) / 2)], + 'obj_places': 5, + 'var_places': 5, + }, + ) + + m.p.set_value(4.0) + x_analytic = 2.0 / math.sqrt(5) + y_analytic = 1.0 / (2.0 * math.sqrt(5)) + _solve_and_check( + self, + self.opt, + m, + { + 'objective': -(x_analytic + y_analytic), + 'vars': [(m.x, x_analytic), (m.y, y_analytic)], + }, + ) + + def test_nl_cubic_constraint_persistent(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, 10)) + m.c = pyo.Constraint(expr=m.x**3 >= 1) + m.obj = pyo.Objective(expr=m.x) + _solve_and_check(self, self.opt, m, {'objective': 1.0, 'vars': [(m.x, 1.0)]}) + + def test_nl_cubic_objective_persistent(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, 1)) + m.obj = pyo.Objective(expr=m.x**3) + _solve_and_check(self, self.opt, m, {'objective': 0.0, 'vars': [(m.x, 0.0)]}) + + +@unittest.pytest.mark.solver('xpress_persistent') +class TestXpressPersistentMisc(unittest.TestCase): + + def setUp(self): + self.opt = XpressPersistent() + + def test_mutable_param_in_objective_coefficient(self): + m = pyo.ConcreteModel() + m.p = pyo.Param(mutable=True, initialize=1.0) + m.x = pyo.Var(bounds=(0, 10)) + m.obj = pyo.Objective(expr=m.p * m.x, sense=pyo.maximize) + _solve_check_mutate_check( + self, + self.opt, + m, + {'objective': 10.0, 'vars': [(m.x, 10.0)]}, + m.p, + -1.0, + {'objective': 0.0, 'vars': [(m.x, 0.0)]}, + ) + + def test_mutable_param_as_variable_bound(self): + m = pyo.ConcreteModel() + m.p = pyo.Param(mutable=True, initialize=5.0) + m.x = pyo.Var(bounds=(0, m.p)) + m.obj = pyo.Objective(expr=-m.x) + _solve_check_mutate_check( + self, + self.opt, + m, + {'objective': -5.0, 'vars': [(m.x, 5.0)]}, + m.p, + 3.0, + {'objective': -3.0, 'vars': [(m.x, 3.0)]}, + ) + + def test_has_instance(self): + self.assertFalse(self.opt.has_instance()) + m = _trivial_model() + self.opt.set_instance(m) + self.assertTrue(self.opt.has_instance()) + self.opt.release() + self.assertFalse(self.opt.has_instance()) + + def test_add_variables_empty_list(self): + m = _trivial_model() + self.opt.set_instance(m) + ncols_before = self.opt._xp_prob.attributes.cols + self.opt.add_variables([]) + self.assertEqual(self.opt._xp_prob.attributes.cols, ncols_before) + + def test_add_constraints_empty_list(self): + m = _trivial_model() + self.opt.set_instance(m) + self.opt.add_constraints([]) + self.opt._add_constraints([]) + + def test_add_block_sos_only(self): + m = pyo.ConcreteModel() + m.x = pyo.Var([1, 2, 3], domain=pyo.NonNegativeReals, bounds=(0, 1)) + m.obj = pyo.Objective(expr=m.x[1] + 2 * m.x[2] + 3 * m.x[3], sense=pyo.maximize) + self.opt.set_instance(m) + m.b = pyo.Block() + m.b.sos = pyo.SOSConstraint(var=m.x, sos=1, weights={1: 1.0, 2: 2.0, 3: 3.0}) + self.opt.add_block(m.b) + _solve_and_check( + self, + self.opt, + m, + {'objective': 3.0, 'vars': [(m.x[1], 0.0), (m.x[2], 0.0), (m.x[3], 1.0)]}, + ) + + def test_remove_block_sos_only(self): + m = pyo.ConcreteModel() + m.x = pyo.Var([1, 2, 3], domain=pyo.NonNegativeReals, bounds=(0, 1)) + m.obj = pyo.Objective(expr=m.x[1] + 2 * m.x[2] + 3 * m.x[3], sense=pyo.maximize) + m.b = pyo.Block() + m.b.sos = pyo.SOSConstraint(var=m.x, sos=1, weights={1: 1.0, 2: 2.0, 3: 3.0}) + self.opt.set_instance(m) + _solve_and_check( + self, + self.opt, + m, + {'objective': 3.0, 'vars': [(m.x[1], 0.0), (m.x[2], 0.0), (m.x[3], 1.0)]}, + ) + self.opt.remove_block(m.b) + m.b.deactivate() + _solve_and_check( + self, + self.opt, + m, + {'objective': 6.0, 'vars': [(m.x[1], 1.0), (m.x[2], 1.0), (m.x[3], 1.0)]}, + ) + + def test_objective_sense_change_only(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, 5)) + m.obj = pyo.Objective(expr=m.x, sense=pyo.minimize) + _solve_and_check(self, self.opt, m, {'objective': 0.0, 'vars': [(m.x, 0.0)]}) + m.obj.sense = pyo.maximize + _solve_and_check(self, self.opt, m, {'objective': 5.0, 'vars': [(m.x, 5.0)]}) + + def test_constant_objective_persistent(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, 10)) + m.c = pyo.Constraint(expr=m.x >= 1) + m.obj = pyo.Objective(expr=7.0) + _solve_and_check(self, self.opt, m, {'objective': 7.0, 'vars': [(m.x, 1.0)]}) + + def test_range_constraint_persistent(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(domain=pyo.NonNegativeReals) + m.y = pyo.Var(domain=pyo.NonNegativeReals) + m.c = pyo.Constraint(expr=pyo.inequality(1, m.x + m.y, 3)) + m.obj = pyo.Objective(expr=-2 * m.x - m.y) + _solve_and_check( + self, self.opt, m, {'objective': -6.0, 'vars': [(m.x, 3.0), (m.y, 0.0)]} + ) + + def test_mutable_param_in_range_constraint(self): + m = pyo.ConcreteModel() + m.p = pyo.Param(mutable=True, initialize=1.0) + m.x = pyo.Var(domain=pyo.NonNegativeReals) + m.c = pyo.Constraint(expr=pyo.inequality(m.p, m.x, 5)) + m.obj = pyo.Objective(expr=m.x) + _solve_check_mutate_check( + self, + self.opt, + m, + {'objective': 1.0, 'vars': [(m.x, 1.0)]}, + m.p, + 3.0, + {'objective': 3.0, 'vars': [(m.x, 3.0)]}, + ) + + def test_mutable_param_in_range_constraint_ub(self): + m = pyo.ConcreteModel() + m.p = pyo.Param(mutable=True, initialize=5.0) + m.x = pyo.Var(domain=pyo.NonNegativeReals) + m.c = pyo.Constraint(expr=pyo.inequality(1, m.x, m.p)) + m.obj = pyo.Objective(expr=-m.x) + _solve_check_mutate_check( + self, + self.opt, + m, + {'objective': -5.0, 'vars': [(m.x, 5.0)]}, + m.p, + 3.0, + {'objective': -3.0, 'vars': [(m.x, 3.0)]}, + ) + + def test_range_constraint_lower_bound_direction(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(domain=pyo.NonNegativeReals) + m.y = pyo.Var(domain=pyo.NonNegativeReals) + m.c = pyo.Constraint(expr=pyo.inequality(1, m.x + m.y, 3)) + m.obj = pyo.Objective(expr=m.x + 2 * m.y) + _solve_and_check( + self, self.opt, m, {'objective': 1.0, 'vars': [(m.x, 1.0), (m.y, 0.0)]} + ) + + def test_mutable_param_changes_constraint_coefficient(self): + m = pyo.ConcreteModel() + m.p = pyo.Param(mutable=True, initialize=1.0) + m.x = pyo.Var(bounds=(0, 10)) + m.c = pyo.Constraint(expr=m.p * m.x <= 5) + m.obj = pyo.Objective(expr=-m.x) + _solve_check_mutate_check( + self, + self.opt, + m, + {'objective': -5.0, 'vars': [(m.x, 5.0)]}, + m.p, + 2.0, + {'objective': -2.5, 'vars': [(m.x, 2.5)]}, + ) + + def test_fix_unfix_variable_via_bounds(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, 10)) + m.y = pyo.Var(bounds=(0, 10)) + m.c = pyo.Constraint(expr=m.x + m.y <= 8) + m.obj = pyo.Objective(expr=-2 * m.x - m.y) + _solve_and_check( + self, self.opt, m, {'objective': -16.0, 'vars': [(m.x, 8.0), (m.y, 0.0)]} + ) + m.y.fix(2.0) + _solve_and_check( + self, self.opt, m, {'objective': -14.0, 'vars': [(m.x, 6.0), (m.y, 2.0)]} + ) + m.y.unfix() + _solve_and_check( + self, self.opt, m, {'objective': -16.0, 'vars': [(m.x, 8.0), (m.y, 0.0)]} + ) + + def test_remove_constraint_drops_mutable_helper(self): + m = pyo.ConcreteModel() + m.p = pyo.Param(mutable=True, initialize=1.0) + m.x = pyo.Var(bounds=(0, 10)) + m.c = pyo.Constraint(expr=m.p * m.x <= 5) + m.obj = pyo.Objective(expr=-m.x) + self.opt.set_instance(m) + _solve_and_check(self, self.opt, m, {'objective': -5.0, 'vars': [(m.x, 5.0)]}) + self.assertIn(m.c, self.opt._mutable_helpers) + self.opt.remove_constraints([m.c]) + m.c.deactivate() + self.assertNotIn(m.c, self.opt._mutable_helpers) + m.p.set_value(2.0) + _solve_and_check(self, self.opt, m, {'objective': -10.0, 'vars': [(m.x, 10.0)]}) + + def test_warmstart_column_indices_match_after_variable_removal(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(within=pyo.Binary) + m.y = pyo.Var(within=pyo.Binary) + m.z = pyo.Var(within=pyo.Binary) + m.obj = pyo.Objective(expr=m.x + m.y + m.z) + m.c = pyo.Constraint(expr=m.x + m.y + m.z >= 1) + self.opt.set_instance(m) + self.opt.remove_variables([m.y]) + del m.y + for j, var in enumerate(self.opt._vars): + xp_idx = self.opt._maps.vars[id(var)].index + self.assertEqual( + j, + xp_idx, + f"After variable removal: Python list position {j} != " + f"Xpress column index {xp_idx} for {var.name}", + ) + m.x.set_value(1) + m.z.set_value(0) + self.opt.remove_constraints([m.c]) + m.c.deactivate() + _solve_and_check( + self, self.opt, m, {'objective': 0.0, 'vars': [(m.x, 0.0), (m.z, 0.0)]} + ) + + +@unittest.pytest.mark.solver('xpress_persistent') +class TestXpressPersistentNLP(unittest.TestCase): + """NLP integration tests for the persistent interface.""" + + def setUp(self): + self.opt = XpressPersistent() + + def _check_optimal(self, res): + self.assertEqual( + res.termination_condition, TerminationCondition.convergenceCriteriaSatisfied + ) + + def test_nl_add_constraint_registers_nl_rebuild(self): + m = pyo.ConcreteModel() + m.p = pyo.Param(mutable=True, initialize=2.0) + m.x = pyo.Var(bounds=(0, 10)) + m.c = pyo.Constraint(expr=m.p * pyo.sin(m.x) <= 5) + m.obj = pyo.Objective(expr=m.x) + self.opt.set_instance(m) + self.opt.solve(m) + self.assertIn(m.c, self.opt._mutable_helpers) + self.assertIsNotNone(self.opt._mutable_helpers[m.c]._nl_expr) + + def test_nl_add_constraint_always_registered(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, math.pi)) + m.c = pyo.Constraint(expr=pyo.sin(m.x) <= 0.5) + m.obj = pyo.Objective(expr=m.x) + self.opt.set_instance(m) + self.opt.solve(m) + self.assertIn(m.c, self.opt._mutable_helpers) + helper = self.opt._mutable_helpers[m.c] + self.assertIsNotNone(helper._nl_expr) + self.assertEqual(len(helper._lin_coefs), 0) + self.assertEqual(len(helper._quad_coefs), 0) + + def test_nl_remove_constraint_cleans_up(self): + m = pyo.ConcreteModel() + m.p = pyo.Param(mutable=True, initialize=2.0) + m.x = pyo.Var(bounds=(0, 10)) + m.c = pyo.Constraint(expr=m.p * pyo.sin(m.x) <= 5) + m.obj = pyo.Objective(expr=m.x) + self.opt.set_instance(m) + self.opt.solve(m) + self.assertIn(m.c, self.opt._mutable_helpers) + self.opt.remove_constraints([m.c]) + m.c.deactivate() + self.assertNotIn(m.c, self.opt._mutable_helpers) + + def test_nl_mutable_linear_coef_in_nl_constraint(self): + m = pyo.ConcreteModel() + m.p = pyo.Param(mutable=True, initialize=1.0) + m.x = pyo.Var(bounds=(0, math.pi)) + m.y = pyo.Var(bounds=(0, 10)) + m.c = pyo.Constraint(expr=pyo.sin(m.x) + m.p * m.y <= 5) + m.obj = pyo.Objective(expr=-m.y) + self.opt.set_instance(m) + _solve_and_check( + self, self.opt, m, {'objective': -5.0, 'vars': [(m.x, 0.0), (m.y, 5.0)]} + ) + y1 = pyo.value(m.y) + m.p.set_value(2.0) + _solve_and_check( + self, self.opt, m, {'objective': -2.5, 'vars': [(m.x, 0.0), (m.y, 2.5)]} + ) + y2 = pyo.value(m.y) + self.assertLess(y2, y1) + + def test_nl_mutable_nl_coef_full_rebuild(self): + m = pyo.ConcreteModel() + m.p = pyo.Param(mutable=True, initialize=1.0) + m.x = pyo.Var(bounds=(0, math.pi / 2)) + m.c = pyo.Constraint(expr=m.p * pyo.sin(m.x) <= 0.5) + m.obj = pyo.Objective(expr=-m.x) + self.opt.set_instance(m) + _solve_and_check( + self, + self.opt, + m, + {'objective': -math.asin(0.5), 'vars': [(m.x, math.asin(0.5))]}, + ) + x1 = pyo.value(m.x) + m.p.set_value(2.0) + _solve_and_check( + self, + self.opt, + m, + {'objective': -math.asin(0.25), 'vars': [(m.x, math.asin(0.25))]}, + ) + x2 = pyo.value(m.x) + self.assertLess(x2, x1 - 0.1) + + def test_nl_mutable_bound(self): + m = pyo.ConcreteModel() + m.p = pyo.Param(mutable=True, initialize=0.5) + m.x = pyo.Var(bounds=(0, math.pi / 2)) + m.c = pyo.Constraint(expr=pyo.sin(m.x) >= m.p) + m.obj = pyo.Objective(expr=m.x) + self.opt.set_instance(m) + _solve_and_check( + self, + self.opt, + m, + {'objective': math.asin(0.5), 'vars': [(m.x, math.asin(0.5))]}, + ) + x1 = pyo.value(m.x) + + m.p.set_value(0.9) + _solve_and_check( + self, + self.opt, + m, + {'objective': math.asin(0.9), 'vars': [(m.x, math.asin(0.9))]}, + ) + x2 = pyo.value(m.x) + self.assertGreater(x2, x1 + 0.4) + + def test_nl_solve_modify_resolve(self): + m = pyo.ConcreteModel() + m.p = pyo.Param(mutable=True, initialize=1.0) + m.x = pyo.Var(bounds=(0, 3)) + m.c = pyo.Constraint(expr=pyo.exp(m.x) >= m.p) + m.obj = pyo.Objective(expr=m.x) + self.opt.set_instance(m) + _solve_check_mutate_check( + self, + self.opt, + m, + {'objective': 0.0, 'vars': [(m.x, 0.0)]}, + m.p, + math.e, + {'objective': 1.0, 'vars': [(m.x, 1.0)]}, + ) + + def test_fix_variable_no_nl_constraint_rebuild(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(0, math.pi)) + m.y = pyo.Var(bounds=(0, 10)) + m.c = pyo.Constraint(expr=pyo.sin(m.x) + m.y <= 5) + m.obj = pyo.Objective(expr=m.y) + self.opt.set_instance(m) + _solve_and_check( + self, self.opt, m, {'objective': 0.0, 'vars': [(m.x, 0.0), (m.y, 0.0)]} + ) + xp_con_before = self.opt._mutable_helpers[m.c]._xp_con + nrows_before = self.opt._xp_prob.attributes.rows + m.x.fix(math.pi / 6) + _solve_and_check( + self, + self.opt, + m, + {'objective': 0.0, 'vars': [(m.x, math.pi / 6), (m.y, 0.0)]}, + ) + nrows_after = self.opt._xp_prob.attributes.rows + self.assertEqual(nrows_before, nrows_after) + self.assertIs(self.opt._mutable_helpers[m.c]._xp_con, xp_con_before) + + def _run_nl_linear_shared_param_test(self, nl_first: bool): + """Test same param in NL (rebuild) and linear (chgMCoef) constraints. + + nl_first=True: NL at row 0, linear at row 1 (delConstraint shifts linear). + nl_first=False: linear at row 0, NL at row 1 (no shift on linear). + """ + m = pyo.ConcreteModel() + m.p = pyo.Param(mutable=True, initialize=1.0) + m.x = pyo.Var(bounds=(0, math.pi / 2)) + m.y = pyo.Var(bounds=(0, 10)) + if nl_first: + m.c_nl = pyo.Constraint(expr=m.p * pyo.sin(m.x) <= 0.5) + m.c_lin = pyo.Constraint(expr=m.p * m.y <= 4) + else: + m.c_lin = pyo.Constraint(expr=m.p * m.y <= 4) + m.c_nl = pyo.Constraint(expr=m.p * pyo.sin(m.x) <= 0.5) + m.obj = pyo.Objective(expr=m.x + m.y, sense=pyo.maximize) + self.opt.set_instance(m) + _solve_and_check( + self, + self.opt, + m, + { + 'objective': 4.0 + math.asin(0.5), + 'vars': [(m.y, 4.0), (m.x, math.asin(0.5))], + }, + ) + + m.p.set_value(2.0) + _solve_and_check( + self, + self.opt, + m, + { + 'objective': 2.0 + math.asin(0.25), + 'vars': [(m.y, 2.0), (m.x, math.asin(0.25))], + }, + ) + + def test_param_shared_nl_before_linear(self): + self._run_nl_linear_shared_param_test(nl_first=True) + + def test_param_shared_linear_before_nl(self): + self._run_nl_linear_shared_param_test(nl_first=False) + + def test_param_shared_multiple_nl_and_linear(self): + m = pyo.ConcreteModel() + m.p = pyo.Param(mutable=True, initialize=1.0) + m.x = pyo.Var(bounds=(0, math.pi)) + m.y = pyo.Var(bounds=(0, 10)) + m.z = pyo.Var(bounds=(0, 10)) + m.c_nl1 = pyo.Constraint(expr=m.p * pyo.sin(m.x) <= 5) + m.c_lin1 = pyo.Constraint(expr=m.p * m.y <= 4) + m.c_nl2 = pyo.Constraint(expr=m.p * pyo.cos(m.x) >= -1) + m.c_lin2 = pyo.Constraint(expr=m.p * m.z <= 3) + m.obj = pyo.Objective(expr=m.y + m.z, sense=pyo.maximize) + self.opt.set_instance(m) + _solve_and_check( + self, + self.opt, + m, + {'objective': 7.0, 'vars': [(m.x, math.pi / 2), (m.y, 4.0), (m.z, 3.0)]}, + ) + self.assertAlmostEqual(pyo.value(m.y) + pyo.value(m.z), 7.0, places=6) + + m.p.set_value(2.0) + _solve_and_check( + self, + self.opt, + m, + {'objective': 3.5, 'vars': [(m.x, math.pi / 3), (m.y, 2.0), (m.z, 1.5)]}, + ) + self.assertAlmostEqual(pyo.value(m.y) + pyo.value(m.z), 3.5, places=6) + + def test_param_quadratic_rebuild_before_linear_collect(self): + m = pyo.ConcreteModel() + m.p = pyo.Param(mutable=True, initialize=1.0) + m.x = pyo.Var(bounds=(0, 5)) + m.y = pyo.Var(bounds=(0, 10)) + m.c_quad = pyo.Constraint(expr=m.p * m.x**2 <= 5) + + m.c_lin = pyo.Constraint(expr=m.p * m.y <= 4) + m.obj = pyo.Objective(expr=m.y, sense=pyo.maximize) + self.opt.set_instance(m) + _solve_and_check( + self, self.opt, m, {'objective': 4.0, 'vars': [(m.x, 0.0), (m.y, 4.0)]} + ) + + m.p.set_value(2.0) + _solve_and_check( + self, self.opt, m, {'objective': 2.0, 'vars': [(m.x, 0.0), (m.y, 2.0)]} + ) + + def test_nl_formula_mutable_and_affine_mutable_in_same_constraint(self): + m = pyo.ConcreteModel() + m.p = pyo.Param(mutable=True, initialize=1.0) + m.q = pyo.Param(mutable=True, initialize=1.0) + m.x = pyo.Var(bounds=(0, math.pi)) + m.y = pyo.Var(bounds=(0, 10)) + m.c = pyo.Constraint(expr=pyo.sin(m.p * m.x) + m.q * m.y <= 5) + m.obj = pyo.Objective(expr=-m.y) + self.opt.set_instance(m) + _solve_and_check( + self, self.opt, m, {'objective': -5.0, 'vars': [(m.x, 0.0), (m.y, 5.0)]} + ) + self.assertIn(m.c, self.opt._mutable_helpers) + self.assertIsNotNone(self.opt._mutable_helpers[m.c]._nl_expr) + pyo.value(m.y) + + m.q.set_value(2.0) + _solve_and_check( + self, self.opt, m, {'objective': -2.5, 'vars': [(m.x, 0.0), (m.y, 2.5)]} + ) + pyo.value(m.y) + + m.p.set_value(0.0) + _solve_and_check( + self, self.opt, m, {'objective': -2.5, 'vars': [(m.x, 0.0), (m.y, 2.5)]} + ) + pyo.value(m.y) + + def test_mutable_objective_does_not_interfere_with_linear_update(self): + m = pyo.ConcreteModel() + m.p = pyo.Param(mutable=True, initialize=1.0) + m.x = pyo.Var(bounds=(0, 3)) + m.y = pyo.Var(bounds=(0, 10)) + m.c_nl = pyo.Constraint(expr=m.p * pyo.sin(m.x) <= 2) + m.c_lin = pyo.Constraint(expr=m.p * m.y <= 4) + m.obj = pyo.Objective(expr=m.p * m.y, sense=pyo.maximize) + self.opt.set_instance(m) + _solve_and_check( + self, self.opt, m, {'objective': 4.0, 'vars': [(m.x, 1.5), (m.y, 4.0)]} + ) + + m.p.set_value(2.0) + _solve_and_check( + self, self.opt, m, {'objective': 4.0, 'vars': [(m.x, 1.5), (m.y, 2.0)]} + ) + + def test_nl_mutable_objective_persistent(self): + m = pyo.ConcreteModel() + m.p = pyo.Param(mutable=True, initialize=1.0) + m.x = pyo.Var(bounds=(0, math.pi / 2)) + m.obj = pyo.Objective(expr=m.p * pyo.sin(m.x), sense=pyo.maximize) + _solve_and_check( + self, self.opt, m, {'objective': 1.0, 'vars': [(m.x, math.pi / 2)]} + ) + pyo.value(m.x) + + m.p.set_value(-1.0) + _solve_and_check(self, self.opt, m, {'objective': 0.0, 'vars': [(m.x, 0.0)]}) + pyo.value(m.x) + + def test_nl_objective_stable_xp_not_mutated_by_constant_update(self): + import math + + m = pyo.ConcreteModel() + m.p = pyo.Param(mutable=True, initialize=1.0) + m.x = pyo.Var(bounds=(0, 1)) + m.z = pyo.Var(bounds=(0, math.pi / 2)) + m.obj = pyo.Objective(expr=2 * m.x + m.p + pyo.sin(m.z)) + _solve_and_check( + self, self.opt, m, {'objective': 1.0, 'vars': [(m.x, 0.0), (m.z, 0.0)]} + ) + + m.p.set_value(3.0) + _solve_and_check( + self, self.opt, m, {'objective': 3.0, 'vars': [(m.x, 0.0), (m.z, 0.0)]} + ) + + m.p.set_value(5.0) + _solve_and_check( + self, self.opt, m, {'objective': 5.0, 'vars': [(m.x, 0.0), (m.z, 0.0)]} + ) + + def test_nl_constraint_stable_quadratic_term(self): + m = pyo.ConcreteModel() + m.p = pyo.Param(mutable=True, initialize=1.0) + m.x = pyo.Var(bounds=(0, 3)) + m.y = pyo.Var(bounds=(0, 1)) + m.c = pyo.Constraint(expr=2 * m.x**2 + m.p * pyo.exp(m.y) <= 5) + m.obj = pyo.Objective(expr=-m.x) + _solve_and_check( + self, + self.opt, + m, + {'objective': -math.sqrt(2), 'vars': [(m.x, math.sqrt(2)), (m.y, 0.0)]}, + ) + x1 = pyo.value(m.x) + + m.p.set_value(2.0) + _solve_and_check( + self, + self.opt, + m, + {'objective': -math.sqrt(1.5), 'vars': [(m.x, math.sqrt(1.5)), (m.y, 0.0)]}, + ) + x2 = pyo.value(m.x) + self.assertLess(x2, x1 - 0.1) + + def test_nl_objective_stable_lin_quad_terms(self): + m = pyo.ConcreteModel() + m.p = pyo.Param(mutable=True, initialize=1.0) + m.x = pyo.Var(bounds=(0, 1)) + m.y = pyo.Var(bounds=(0, 1)) + m.z = pyo.Var(bounds=(0, math.pi / 2)) + m.obj = pyo.Objective(expr=2 * m.x + m.y**2 + m.p * pyo.sin(m.z)) + _solve_and_check( + self, + self.opt, + m, + {'objective': 0.0, 'vars': [(m.x, 0.0), (m.y, 0.0), (m.z, 0.0)]}, + ) + + m.p.set_value(-1.0) + _solve_and_check( + self, + self.opt, + m, + {'objective': -1.0, 'vars': [(m.x, 0.0), (m.y, 0.0), (m.z, math.pi / 2)]}, + ) + + def test_nl_cubic_constraint_mutable_param_persistent(self): + m = pyo.ConcreteModel() + m.p = pyo.Param(mutable=True, initialize=1.0) + m.x = pyo.Var(bounds=(0, 10)) + m.c = pyo.Constraint(expr=m.p * m.x**3 >= 1) + m.obj = pyo.Objective(expr=m.x) + _solve_check_mutate_check( + self, + self.opt, + m, + {'objective': 1.0, 'vars': [(m.x, 1.0)]}, + m.p, + 8.0, + {'objective': 0.5, 'vars': [(m.x, 0.5)]}, + ) + + def test_add_remove_readd_changes_row_ordering(self): + m = pyo.ConcreteModel() + m.p = pyo.Param(mutable=True, initialize=1.0) + m.x = pyo.Var(bounds=(0, math.pi)) + m.y = pyo.Var(bounds=(0, 10)) + m.c_lin = pyo.Constraint(expr=m.p * m.y <= 4) + m.c_nl = pyo.Constraint(expr=m.p * pyo.sin(m.x) <= 5) + m.obj = pyo.Objective(expr=m.y, sense=pyo.maximize) + self.opt.set_instance(m) + _solve_and_check( + self, + self.opt, + m, + {'objective': 4.0, 'vars': [(m.x, math.pi / 2), (m.y, 4.0)]}, + ) + + m.c_lin.deactivate() + self.opt.remove_constraints([m.c_lin]) + m.c_lin.activate() + self.opt.add_constraints([m.c_lin]) + self.assertIn(m.c_lin, self.opt._mutable_helpers) + + m.p.set_value(2.0) + _solve_and_check( + self, + self.opt, + m, + {'objective': 2.0, 'vars': [(m.x, math.pi / 2), (m.y, 2.0)]}, + ) + + +@unittest.pytest.mark.solver('xpress_persistent') +class TestXpressPersistentIIS(unittest.TestCase): + + def setUp(self): + self.opt = XpressPersistent() + + def _infeasible_model(self): + m = pyo.ConcreteModel() + m.x = pyo.Var(within=pyo.Binary) + m.y = pyo.Var(within=pyo.NonNegativeReals) + m.c1 = pyo.Constraint(expr=m.y <= 100.0 * m.x) + m.c2 = pyo.Constraint(expr=m.y <= -100.0 * m.x) + m.c3 = pyo.Constraint(expr=m.x >= 0.5) + m.obj = pyo.Objective(expr=-m.y) + return m + + def test_write_iis_produces_file(self): + import os, tempfile + + m = self._infeasible_model() + self.opt.solve( + m, + raise_exception_on_nonoptimal_result=False, + load_solutions=False, + symbolic_solver_labels=True, + ) + with tempfile.TemporaryDirectory() as tmp: + base = os.path.join(tmp, 'iis') + result = self.opt.write_iis(base) + self.assertEqual(result, base) + lp_file = base + '.lp' + self.assertTrue(os.path.exists(lp_file)) + with open(lp_file) as f: + content = f.read() + self.assertIn('c2', content) + self.assertIn('c3', content) + + def test_get_iis_returns_pyomo_objects(self): + m = self._infeasible_model() + self.opt.solve( + m, raise_exception_on_nonoptimal_result=False, load_solutions=False + ) + iis = self.opt.get_iis() + self.assertIn('constraints', iis) + self.assertIn('variables', iis) + con_names = {c.name for c in iis['constraints']} + self.assertIn('c2', con_names) + self.assertIn('c3', con_names) + var_names = {v.name for v in iis['variables']} + self.assertIn('y', var_names) + + def test_get_iis_objects_are_model_constraints(self): + m = self._infeasible_model() + self.opt.solve( + m, raise_exception_on_nonoptimal_result=False, load_solutions=False + ) + iis = self.opt.get_iis() + model_cons = list(m.component_data_objects(pyo.Constraint, active=True)) + model_vars = list(m.component_data_objects(pyo.Var)) + for con in iis['constraints']: + self.assertTrue( + any(con is c for c in model_cons), + f"{con.name} is not a model constraint object", + ) + for var in iis['variables']: + self.assertTrue( + any(var is v for v in model_vars), + f"{var.name} is not a model variable object", + ) + + +if __name__ == '__main__': + unittest.main() diff --git a/pyomo/contrib/solver/tests/solvers/test_xpress_pool.py b/pyomo/contrib/solver/tests/solvers/test_xpress_pool.py new file mode 100644 index 00000000000..c354ee5a410 --- /dev/null +++ b/pyomo/contrib/solver/tests/solvers/test_xpress_pool.py @@ -0,0 +1,179 @@ +# ____________________________________________________________________________________ +# +# 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.contrib.solver.common.util import ( + NoDualsError, + NoReducedCostsError, + NoSolutionError, +) +from pyomo.contrib.solver.common.results import SolutionStatus +from pyomo.contrib.solver.solvers.xpress import XpressDirect, XpressPersistent +from pyomo.contrib.solver.tests.solvers._xpress_test_utils import _simple_lp + +if not XpressDirect().available(): + raise unittest.SkipTest('Xpress not available') + + +def _make_mip_with_many_solutions(): + """Small binary MIP with many feasible integer solutions. + + max x[0] + x[1] + x[2] + x[3] + x[4] + s.t. x[0] + x[1] + x[2] + x[3] + x[4] <= 3 + x in {0,1}^5 + + Optimal value = 3 (choose any 3 of 5 items). There are C(5,3) = 10 + feasible solutions with obj=3, plus all sub-optimal assignments. + """ + m = pyo.ConcreteModel() + m.I = pyo.RangeSet(0, 4) + m.x = pyo.Var(m.I, within=pyo.Binary) + m.c = pyo.Constraint(expr=sum(m.x[i] for i in m.I) <= 3) + m.obj = pyo.Objective(expr=sum(m.x[i] for i in m.I), sense=pyo.maximize) + return m + + +@unittest.pytest.mark.solver('xpress_direct') +class TestXpressSolutionPool(unittest.TestCase): + """Integration tests for the solution pool (pool_solutions config).""" + + def setUp(self): + self.opt = XpressDirect() + + def test_pool_disabled_by_default(self): + """Default config collects no pool: exactly 1 solution accessible.""" + m = _make_mip_with_many_solutions() + res = self.opt.solve(m) + loader = res.solution_loader + self.assertEqual(loader.get_number_of_solutions(), 1) + self.assertEqual(loader.get_solution_ids(), [0]) + self.assertEqual(res.solution_status, SolutionStatus.optimal) + + def test_pool_collects_solutions(self): + """pool_solutions=5 collects multiple feasible solutions. + + Verifies: + - get_number_of_solutions() > 1 + - solution 0 is the optimal incumbent + - solution 1 is a distinct valid assignment + """ + m = _make_mip_with_many_solutions() + res = self.opt.solve(m, pool_solutions=5, load_solutions=False) + loader = res.solution_loader + + n = loader.get_number_of_solutions() + self.assertGreater(n, 1) + self.assertEqual(loader.get_solution_ids(), list(range(n))) + + # Load and check solution 0 (incumbent -- optimal). + loader.solution(0).load_vars() + obj0 = pyo.value(m.obj) + self.assertAlmostEqual(obj0, res.incumbent_objective, places=6) + self.assertEqual(obj0, 3.0) + + # Load solution 1 and verify it is a valid feasible assignment. + loader.solution(1).load_vars() + for i in m.I: + self.assertIn(round(pyo.value(m.x[i])), (0, 1)) + total = sum(pyo.value(m.x[i]) for i in m.I) + self.assertLessEqual(total, 3.0 + 1e-6) + + def test_pool_context_manager(self): + """solution(k) context manager restores incumbent on exit.""" + m = _make_mip_with_many_solutions() + res = self.opt.solve(m, pool_solutions=5, load_solutions=False) + loader = res.solution_loader + + if loader.get_number_of_solutions() < 2: + self.skipTest('Solver found fewer than 2 solutions -- cannot test pool.') + + # Load incumbent into model variables. + loader.solution(0).load_vars() + incumbent_vals = {i: pyo.value(m.x[i]) for i in m.I} + + # Context manager temporarily activates solution 1. + with loader.solution(1): + loader.load_vars() + + # After exiting the context, active id is restored to 0 (incumbent). + # Reload to confirm values match the original incumbent. + loader.load_vars() + for i in m.I: + self.assertAlmostEqual(pyo.value(m.x[i]), incumbent_vals[i], places=6) + + def test_pool_solution_raises_out_of_range(self): + """Requesting solution(1) on a default (no-pool) solve must raise NoSolutionError. + Without the fix, the silent fallthrough would load incumbent values instead.""" + m = _make_mip_with_many_solutions() + res = self.opt.solve(m, load_solutions=False) + loader = res.solution_loader + self.assertEqual(loader.get_number_of_solutions(), 1) + with self.assertRaises(NoSolutionError): + loader.solution(1).load_vars() + + def test_pool_duals_raise_for_nonzero_id(self): + """get_duals() and get_reduced_costs() inside solution(1) context must raise.""" + m = _make_mip_with_many_solutions() + # Use LP model for duals (MIP has no duals anyway); LP has an incumbent only. + lp = _simple_lp() + self.opt.solve(lp, pool_solutions=0, load_solutions=False) + # Pool is empty; solution(1) will raise NoSolutionError, which already + # confirms the guard fires. For the duals/RC guard test we need pool_solutions>0 + # so a solution(1) context can be entered. Use the MIP model for that. + mip_res = self.opt.solve(m, pool_solutions=5, load_solutions=False) + loader = mip_res.solution_loader + if loader.get_number_of_solutions() < 2: + self.skipTest( + 'Solver found fewer than 2 solutions -- cannot test duals guard.' + ) + with loader.solution(1): + with self.assertRaises(NoDualsError): + loader.get_duals() + with self.assertRaises(NoReducedCostsError): + loader.get_reduced_costs() + + def test_pool_persistent(self): + """pool_solutions works through XpressPersistent and callbacks do not + accumulate across consecutive solves.""" + opt = XpressPersistent() + m = _make_mip_with_many_solutions() + + res1 = opt.solve(m, pool_solutions=5) + n1 = res1.solution_loader.get_number_of_solutions() + self.assertGreaterEqual(n1, 1) + + # Second solve: pool must be freshly collected, not accumulated from solve 1. + res2 = opt.solve(m, pool_solutions=5) + n2 = res2.solution_loader.get_number_of_solutions() + self.assertGreaterEqual(n2, 1) + # The pool from the first solve is in res1's loader; res2 has its own pool. + self.assertIsNot(res1.solution_loader, res2.solution_loader) + + def test_pool_rolling_window(self): + """pool_solutions=N keeps a rolling window of the last N solutions found. + + The pool size must not exceed N even when more than N solutions are found + during B&B. Uses a model with many feasible integer solutions. + """ + m = _make_mip_with_many_solutions() + window = 2 + res = self.opt.solve(m, pool_solutions=window, load_solutions=False) + loader = res.solution_loader + n = loader.get_number_of_solutions() + # Pool entries are at most window + 1 (incumbent + window collected) + # but the solver might have found fewer than window non-incumbent solutions. + self.assertGreaterEqual(n, 1) + self.assertLessEqual(n, window + 1) + + +if __name__ == '__main__': + unittest.main() diff --git a/pyomo/contrib/solver/tests/solvers/test_xpress_walker.py b/pyomo/contrib/solver/tests/solvers/test_xpress_walker.py new file mode 100644 index 00000000000..bda1b6049d9 --- /dev/null +++ b/pyomo/contrib/solver/tests/solvers/test_xpress_walker.py @@ -0,0 +1,668 @@ +# ____________________________________________________________________________________ +# +# 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. +# ____________________________________________________________________________________ + +"""Walker unit tests: XpressExpressionWalker xp expression building.""" + +import pyomo.common.unittest as unittest +import pyomo.environ as pyo + +from parameterized import parameterized + +from pyomo.common.dependencies import attempt_import +from pyomo.contrib.solver.common.util import IncompatibleModelError +from pyomo.core.expr import sin, ceil +from pyomo.core.expr.numeric_expr import LinearExpression +from pyomo.contrib.solver.solvers.xpress.xpress_base import _before_linear +from pyomo.core.expr import MinExpression, MaxExpression +from pyomo.core.expr.numvalue import NumericConstant +from pyomo.environ import Expr_if +from pyomo.contrib.solver.solvers.xpress.xpress_base import _EXIT_HANDLERS + +xp, xpress_available = attempt_import('xpress', catch_exceptions=(Exception,)) + +try: + from pyomo.contrib.solver.solvers.xpress.xpress_base import XpressExpressionWalker +except Exception: + pass + +if not xpress_available: + raise unittest.SkipTest('Xpress not available') + + +def _make_walker(): + """Return (model, walker, xv, yv, zv) with x/y/z vars registered.""" + m = pyo.ConcreteModel() + m.x = pyo.Var() + m.y = pyo.Var() + m.z = pyo.Var() + prob = xp.problem() + xv = prob.addVariable(name='x') + yv = prob.addVariable(name='y') + zv = prob.addVariable(name='z') + var_map = {id(m.x): xv, id(m.y): yv, id(m.z): zv} + walker = XpressExpressionWalker(var_map, prob) + return m, walker, xv, yv, zv + + +def _walk_and_assert( + test_case, + walker, + expr, + expected_type, + expected_lin=[], + expected_quad=[], + expected_nlp_vars=[], + expected_nlp_scalars=[], + expected_nlp_operators=[], +): + """Walk a Pyomo expression and assert structural properties of the result.""" + result = walker.walk_expression(expr) + + tc = test_case + var_map = walker.var_map + + tc.assertIsInstance(result, expected_type) + + lin_vars, lin_coefs = result.extractLinear() + lin_by_id = {id(v): c for v, c in zip(lin_vars, lin_coefs)} + tc.assertEqual(len(lin_by_id), len(expected_lin)) + for pyo_var, expected_coef in expected_lin: + xp_var = var_map[id(pyo_var)] + tc.assertIn(id(xp_var), lin_by_id, f"{pyo_var.name} missing from linear terms") + tc.assertAlmostEqual(lin_by_id[id(xp_var)], expected_coef, places=12) + + def qkey(v1, v2): + id1, id2 = id(v1), id(v2) + return (id1, id2) if id1 < id2 else (id2, id1) + + q_v1s, q_v2s, q_coefs = result.extractQuadratic() + quad_by_ids = {qkey(v1, v2): c for v1, v2, c in zip(q_v1s, q_v2s, q_coefs)} + tc.assertEqual(len(q_coefs), len(expected_quad)) + for pyo_v1, pyo_v2, expected_coef in expected_quad: + key = qkey(var_map[id(pyo_v1)], var_map[id(pyo_v2)]) + tc.assertIn( + key, + quad_by_ids, + f"({pyo_v1.name}, {pyo_v2.name}) missing from quadratic terms", + ) + tc.assertAlmostEqual(quad_by_ids[key], expected_coef, places=12) + + test_con = xp.constraint(body=result, lb=-1e20, ub=1e20) + n_before = walker.prob.attributes.rows + walker.prob.addConstraint(test_con) + tc.assertEqual(walker.prob.attributes.rows, n_before + 1) + tc.assertEqual(test_con, walker.prob.getConstraint(n_before)) + + types, values = walker.prob.nlpGetFormula(test_con, 0) + + xc = xp.constants + _IFUN = { + 'log10': xc.IFUN_LOG10, + 'ln': xc.IFUN_LN, + 'exp': xc.IFUN_EXP, + 'abs': xc.IFUN_ABS, + 'sqrt': xc.IFUN_SQRT, + 'sin': xc.IFUN_SIN, + 'cos': xc.IFUN_COS, + 'tan': xc.IFUN_TAN, + 'asin': xc.IFUN_ARCSIN, + 'arcsin': xc.IFUN_ARCSIN, + 'acos': xc.IFUN_ARCCOS, + 'arccos': xc.IFUN_ARCCOS, + 'atan': xc.IFUN_ARCTAN, + 'arctan': xc.IFUN_ARCTAN, + 'min': xc.IFUN_MIN, + 'max': xc.IFUN_MAX, + } + _OP = { + 'uminus': xc.OP_UMINUS, + 'pow': xc.OP_EXPONENT, + 'mul': xc.OP_MULTIPLY, + 'div': xc.OP_DIVIDE, + 'plus': xc.OP_PLUS, + 'minus': xc.OP_MINUS, + } + + col_indices = {int(v) for t, v in zip(types, values) if int(t) == xc.TOK_COL} + scalars = [v for t, v in zip(types, values) if int(t) == xc.TOK_CON] + ifun_codes = {int(v) for t, v in zip(types, values) if int(t) == xc.TOK_IFUN} + op_codes = {int(v) for t, v in zip(types, values) if int(t) == xc.TOK_OP} + + for pyo_var in expected_nlp_vars: + xp_var = var_map[id(pyo_var)] + tc.assertIn( + xp_var.index, + col_indices, + f"Variable {pyo_var.name} not found in NLP formula", + ) + + tc.assertEqual(len(scalars), len(expected_nlp_scalars)) + for s in expected_nlp_scalars: + tc.assertTrue( + any(abs(v - s) < 1e-10 for v in scalars), + f"Scalar {s} not found in NLP formula scalars {scalars}", + ) + + for op_name in expected_nlp_operators: + key = op_name.lower() + if key in _IFUN: + tc.assertIn( + _IFUN[key], + ifun_codes, + f"NLP function '{op_name}' (code {_IFUN[key]}) not found in formula", + ) + elif key in _OP: + tc.assertIn( + _OP[key], + op_codes, + f"Operator '{op_name}' (code {_OP[key]}) not found in formula", + ) + else: + raise ValueError(f"Unknown NLP operator name '{op_name}'") + + walker.prob.delConstraint(test_con) + tc.assertEqual(walker.prob.attributes.rows, n_before) + return result + + +@unittest.pytest.mark.solver('xpress_direct') +class TestXpressWalkerLinear(unittest.TestCase): + + def test_linear_float_coef(self): + m, w, _, _, _ = _make_walker() + _walk_and_assert( + self, + w, + 3.0 * m.x + 2.0 * m.y, + expected_type=xp.expression, + expected_lin=[(m.x, 3.0), (m.y, 2.0)], + ) + + def test_linear_mutable_coef(self): + m, w, _, _, _ = _make_walker() + m.p = pyo.Param(mutable=True, initialize=5.0) + _walk_and_assert( + self, + w, + m.p * m.x + m.y, + expected_type=xp.expression, + expected_lin=[(m.x, 5.0), (m.y, 1.0)], + ) + + def test_linear_zero_coef(self): + m, w, _, _, _ = _make_walker() + _walk_and_assert( + self, + w, + 0 * m.x + m.y, + expected_type=xp.expression, + expected_lin=[(m.y, 1.0)], + ) + + def test_fixed_var_kept_as_column(self): + m, w, _, _, _ = _make_walker() + m.x.fix(2.0) + _walk_and_assert( + self, + w, + 3.0 * m.x + m.y, + expected_type=xp.expression, + expected_lin=[(m.x, 3.0), (m.y, 1.0)], + ) + + def test_sum_with_constant(self): + m, w, _, _, _ = _make_walker() + _walk_and_assert( + self, + w, + 3.0 * m.x + 5.0, + expected_type=xp.expression, + expected_lin=[(m.x, 3.0)], + ) + + def test_zero_coef_monomial_is_constant(self): + m, w, _, _, _ = _make_walker() + _walk_and_assert( + self, + w, + sin(m.x) + 0 * m.y, + expected_type=xp.nonlin, + expected_lin=[], + expected_quad=[], + ) + + def test_mutable_constant_body_const(self): + m, w, _, _, _ = _make_walker() + m.p = pyo.Param(mutable=True, initialize=3.0) + _walk_and_assert( + self, w, m.x + m.p, expected_type=xp.expression, expected_lin=[(m.x, 1.0)] + ) + + def test_negation(self): + m, w, _, _, _ = _make_walker() + _walk_and_assert( + self, w, -m.x, expected_type=xp.linterm, expected_lin=[(m.x, -1.0)] + ) + + def test_division_by_constant(self): + m, w, _, _, _ = _make_walker() + _walk_and_assert( + self, w, m.x / 2.0, expected_type=xp.linterm, expected_lin=[(m.x, 0.5)] + ) + + def test_linear_two_vars_no_coef(self): + m, w, _, _, _ = _make_walker() + _walk_and_assert( + self, + w, + m.x + m.y, + expected_type=xp.expression, + expected_lin=[(m.x, 1.0), (m.y, 1.0)], + ) + + def test_all_constant_expr(self): + m, w, _, _, _ = _make_walker() + result = w.walk_expression(3.0 + 2.0) + self.assertAlmostEqual(result, 5.0) + + def test_linear_fast_path(self): + m, w, _, _, _ = _make_walker() + expr = m.x + m.y + self.assertIsInstance(expr, LinearExpression) + _, result = _before_linear(w, expr) + self.assertIsNotNone(result) + + def test_before_linear_dispatcher_contract(self): + m, w, _, _, _ = _make_walker() + expr = 3.0 * m.x + 2.0 * m.y + self.assertIsInstance(expr, LinearExpression) + should_descend, result = _before_linear(w, expr) + self.assertFalse(should_descend) + self.assertIsNotNone(result) + + +@unittest.pytest.mark.solver('xpress_direct') +class TestXpressWalkerQuadratic(unittest.TestCase): + + def test_power_squared(self): + m, w, _, _, _ = _make_walker() + _walk_and_assert( + self, + w, + m.x**2, + expected_type=xp.quadterm, + expected_lin=[], + expected_quad=[(m.x, m.x, 1.0)], + ) + + def test_product_two_vars(self): + m, w, _, _, _ = _make_walker() + _walk_and_assert( + self, + w, + m.x * m.y, + expected_type=xp.quadterm, + expected_lin=[], + expected_quad=[(m.x, m.y, 1.0)], + ) + + def test_mutable_quad_coef(self): + m, w, _, _, _ = _make_walker() + m.p = pyo.Param(mutable=True, initialize=3.0) + _walk_and_assert( + self, + w, + m.p * m.x * m.y, + expected_type=xp.quadterm, + expected_lin=[], + expected_quad=[(m.x, m.y, 3.0)], + ) + + def test_bilinear_p_times_paren_xy(self): + m, w, _, _, _ = _make_walker() + m.p = pyo.Param(mutable=True, initialize=3.0) + _walk_and_assert( + self, + w, + m.p * (m.x * m.y), + expected_type=xp.quadterm, + expected_lin=[], + expected_quad=[(m.x, m.y, 3.0)], + ) + + def test_monomial_times_var(self): + m, w, _, _, _ = _make_walker() + m.p = pyo.Param(mutable=True, initialize=2.0) + _walk_and_assert( + self, + w, + (m.p * m.x) * m.y, + expected_type=xp.quadterm, + expected_lin=[], + expected_quad=[(m.x, m.y, 2.0)], + ) + + def test_monomial_times_same_var(self): + m, w, _, _, _ = _make_walker() + m.p = pyo.Param(mutable=True, initialize=4.0) + _walk_and_assert( + self, + w, + (m.p * m.x) * m.x, + expected_type=xp.quadterm, + expected_lin=[], + expected_quad=[(m.x, m.x, 4.0)], + ) + + def test_mutable_coef_var_squared(self): + m, w, _, _, _ = _make_walker() + m.p = pyo.Param(mutable=True, initialize=3.0) + _walk_and_assert( + self, + w, + m.p * m.x**2, + expected_type=xp.quadterm, + expected_lin=[], + expected_quad=[(m.x, m.x, 3.0)], + ) + + def test_non_mutable_coef_var_squared(self): + m, w, _, _, _ = _make_walker() + _walk_and_assert( + self, + w, + 3 * m.x**2, + expected_type=xp.quadterm, + expected_lin=[], + expected_quad=[(m.x, m.x, 3.0)], + ) + + +@unittest.pytest.mark.solver('xpress_direct') +class TestXpressWalkerNonlinear(unittest.TestCase): + + def test_division_by_variable(self): + m, w, _, _, _ = _make_walker() + _walk_and_assert( + self, + w, + m.x / m.y, + expected_type=xp.nonlin, + expected_nlp_vars=[m.x, m.y], + expected_nlp_operators=['div'], + ) + + def test_power_cubic_is_nl(self): + m, w, _, _, _ = _make_walker() + _walk_and_assert(self, w, m.x * m.y * m.z, expected_type=xp.nonlin) + + @parameterized.expand( + [ + ('sin', pyo.sin, [], ['sin']), + ('cos', pyo.cos, [], ['cos']), + ('sinh', pyo.sinh, [0.5], ['exp']), + ('cosh', pyo.cosh, [0.5], ['exp']), + ('tanh', pyo.tanh, [], ['exp', 'div']), + ('asinh', pyo.asinh, [1.0], ['ln', 'sqrt']), + ('acosh', pyo.acosh, [-1.0], ['ln', 'sqrt']), + ('atanh', pyo.atanh, [1.0, 1.0, 0.5], ['ln']), + ] + ) + def test_nl_trig_hyperbolic_plus_linear(self, fn_name, fn, nlp_scalars, nlp_ops): + m, w, _, _, _ = _make_walker() + _walk_and_assert( + self, + w, + fn(m.x) + m.y, + expected_type=xp.nonlin, + expected_lin=[(m.y, 1.0)], + expected_quad=[], + expected_nlp_vars=[m.x], + expected_nlp_scalars=nlp_scalars, + expected_nlp_operators=nlp_ops, + ) + + def test_nl_mutable_outside_nl(self): + m, w, _, _, _ = _make_walker() + m.p = pyo.Param(mutable=True, initialize=2.0) + _walk_and_assert( + self, + w, + m.p * sin(m.x), + expected_type=xp.nonlin, + expected_nlp_vars=[m.x], + expected_nlp_scalars=[2.0], + expected_nlp_operators=['sin'], + ) + + def test_nl_mutable_linear_coexist(self): + m, w, _, _, _ = _make_walker() + m.p = pyo.Param(mutable=True, initialize=2.0) + _walk_and_assert( + self, + w, + sin(m.x) + m.p * m.y, + expected_type=xp.nonlin, + expected_lin=[(m.y, 2.0)], + expected_nlp_vars=[m.x], + expected_nlp_scalars=[], + expected_nlp_operators=['sin'], + ) + + def test_nl_mutable_sin_arg(self): + m, w, _, _, _ = _make_walker() + m.p = pyo.Param(mutable=True, initialize=2.0) + _walk_and_assert( + self, + w, + sin(m.p * m.x), + expected_type=xp.nonlin, + expected_nlp_vars=[m.x], + expected_nlp_scalars=[2.0], + expected_nlp_operators=['sin'], + ) + + def test_min_expression(self): + m, w, _, _, _ = _make_walker() + _walk_and_assert( + self, + w, + MinExpression([m.x, m.y, m.z]), + expected_type=xp.nonlin, + expected_nlp_vars=[m.x, m.y, m.z], + expected_nlp_scalars=[], + expected_nlp_operators=['min'], + ) + + def test_abs_expression(self): + m, w, _, _, _ = _make_walker() + _walk_and_assert( + self, + w, + abs(m.x + m.y), + expected_type=xp.nonlin, + expected_nlp_vars=[m.x, m.y], + expected_nlp_scalars=[], + expected_nlp_operators=['abs'], + ) + + def test_max_expression(self): + m, w, _, _, _ = _make_walker() + _walk_and_assert( + self, + w, + MaxExpression([m.x, m.y, m.z]), + expected_type=xp.nonlin, + expected_nlp_vars=[m.x, m.y, m.z], + expected_nlp_scalars=[], + expected_nlp_operators=['max'], + ) + + def test_max_all_constants(self): + _, w, _, _, _ = _make_walker() + _walk_and_assert( + self, + w, + MaxExpression([NumericConstant(1.0), NumericConstant(2.0)]), + expected_type=xp.nonlin, + expected_nlp_scalars=[1.0, 2.0], + expected_nlp_operators=['max'], + ) + + def test_mutable_param_in_max(self): + m, w, _, _, _ = _make_walker() + m.p = pyo.Param(mutable=True, initialize=3.0) + _walk_and_assert( + self, + w, + MaxExpression([m.x, m.p]), + expected_type=xp.nonlin, + expected_nlp_vars=[m.x], + expected_nlp_scalars=[3.0], + expected_nlp_operators=['max'], + ) + + def test_sum_divided_by_constant(self): + m, w, _, _, _ = _make_walker() + _walk_and_assert( + self, + w, + (m.x + m.y) / 2.0, + expected_type=xp.expression, + expected_lin=[(m.x, 0.5), (m.y, 0.5)], + ) + + def test_nl_times_nl(self): + m, w, _, _, _ = _make_walker() + _walk_and_assert( + self, + w, + sin(m.x) * sin(m.y), + expected_type=xp.nonlin, + expected_nlp_vars=[m.x, m.y], + expected_nlp_operators=['sin'], + ) + + def test_mutable_param_in_nl_product(self): + m, w, _, _, _ = _make_walker() + m.p = pyo.Param(mutable=True, initialize=2.0) + _walk_and_assert( + self, + w, + pyo.sin(m.p * m.x) * pyo.sin(m.y), + expected_type=xp.nonlin, + expected_nlp_vars=[m.x, m.y], + expected_nlp_scalars=[2.0], + expected_nlp_operators=['sin'], + ) + + def test_before_npv_evaluation_error_arithmetic(self): + m, w, _, _, _ = _make_walker() + m.q = pyo.Param(mutable=True, initialize=1.0) + m.r = pyo.Param(mutable=True, initialize=0.0) + with self.assertRaises(ZeroDivisionError): + w.walk_expression(sin(m.x) + m.q / m.r) + + def test_exit_unary_const_domain_error(self): + m, w, _, _, _ = _make_walker() + m.s = pyo.Param(mutable=True, initialize=1.0) + expr = sin(m.x) + pyo.sqrt(m.s) + m.s.set_value(-1.0) + with self.assertRaises(ValueError): + w.walk_expression(expr) + + def test_npv_evaluation_error(self): + m, w, _, _, _ = _make_walker() + m.p = pyo.Param(mutable=True, initialize=-1.0) + expr = pyo.log(m.p) + m.x + with self.assertRaises(ValueError): + w.walk_expression(expr) + + +@unittest.pytest.mark.solver('xpress_direct') +class TestXpressWalkerCache(unittest.TestCase): + + def test_named_expr_cache_populated(self): + m, w, _, _, _ = _make_walker() + m.e = pyo.Expression(expr=m.x + m.y) + self.assertEqual(len(w.subexpression_cache), 0) + _walk_and_assert( + self, + w, + m.e + 1.0, + expected_type=xp.expression, + expected_lin=[(m.x, 1.0), (m.y, 1.0)], + ) + self.assertEqual(len(w.subexpression_cache), 1) + + def test_named_expr_cache_hit(self): + m, w, _, _, _ = _make_walker() + m.e = pyo.Expression(expr=m.x + m.y) + self.assertEqual(len(w.subexpression_cache), 0) + _walk_and_assert( + self, + w, + m.e + 1.0, + expected_type=xp.expression, + expected_lin=[(m.x, 1.0), (m.y, 1.0)], + ) + self.assertEqual(len(w.subexpression_cache), 1) + _walk_and_assert( + self, + w, + m.e + 2.0, + expected_type=xp.expression, + expected_lin=[(m.x, 1.0), (m.y, 1.0)], + ) + self.assertEqual(len(w.subexpression_cache), 1) + + def test_named_expr_produces_valid_xp_expr(self): + m, w, _, _, _ = _make_walker() + m.p = pyo.Param(mutable=True, initialize=2.0) + m.e = pyo.Expression(expr=m.p * m.x) + _walk_and_assert( + self, + w, + m.e + m.y, + expected_type=xp.expression, + expected_lin=[(m.x, 2.0), (m.y, 1.0)], + ) + + +@unittest.pytest.mark.solver('xpress_direct') +class TestXpressWalkerErrors(unittest.TestCase): + + def test_unsupported_function_ceil_raises(self): + m, w, _, _, _ = _make_walker() + with self.assertRaises(IncompatibleModelError) as ctx: + w.walk_expression(ceil(m.x)) + self.assertIn('ceil', str(ctx.exception)) + + def test_unsupported_function_floor_raises(self): + m, w, _, _, _ = _make_walker() + with self.assertRaises(IncompatibleModelError) as ctx: + w.walk_expression(pyo.floor(m.x)) + self.assertIn('floor', str(ctx.exception)) + + def test_expr_if_raises(self): + m, w, _, _, _ = _make_walker() + expr = Expr_if(IF=m.x > 0, THEN=m.x, ELSE=-m.x) + with self.assertRaises(IncompatibleModelError): + w.walk_expression(expr) + + def test_unregistered_expression_type_raises(self): + class _CustomUnregisteredExpr: + pass + + with self.assertRaises(IncompatibleModelError): + _EXIT_HANDLERS[_CustomUnregisteredExpr] + + +if __name__ == '__main__': + unittest.main()