From d0d5a736e10f01cc82122ad2d3d98bb511f0c0ec Mon Sep 17 00:00:00 2001 From: Francesco Cavaliere Date: Tue, 30 Jun 2026 12:29:43 +0200 Subject: [PATCH 1/2] Txp 8774 nlp (#3) * TXP-8770: add new Xpress direct solver (LP/MIP core) for new Pyomo solver interface * TXP-8770: optimize and clean up Xpress direct solver post-review * TXP-8770: add Xpress-specific test suite covering controls, labels and error paths * TXP-8770: add warmstart, solver stats and parity/beyond-Gurobi improvements * TXP-8770: tighten test assertions and add solver marker * TXP-8771: add XpressPersistent solver and refactor connector shared code New xpress_base.py (replaces xpress_direct_base.py): - EntityMaps dataclass: groups id(VarData)->xp.var, ConstraintData->list[xp.constraint], SOSConstraintData->xp.sos stable handle maps; avoids integer index rebuilds after deletions - XpressSolutionLoaderBase: unified loader with _fetch_var_values helper, get_duals walking handle spans with max-abs for 2-row range constraints - build_loadmip_arrays, extract_sos_arrays, get_rhs_and_sense: shared compile utilities - build_xpress_problem: full compile pipeline; deactivates SOSConstraint components before LinearStandardFormCompiler (which rejects them) and reactivates in a finally block - XpressSolverMixin: solve scaffold and _create_xpress_model dispatch extracted from XpressDirect, now shared with XpressPersistent New xpress_persistent.py: - XpressPersistent: reuses xp.problem() across solves via ModelChangeDetector observer - Incremental updates: _update_variables (add/remove/modify, fix/unfix reprocesses rows), _update_constraints (remove+re-add on expr change), _update_sos_constraints, _update_objectives (added/removed/sense), _update_parameters - _initializing flag suppresses callbacks during ModelChangeDetector initial snapshot - Fixed variable handling: compiler folds fixed vars into constants; extra zero-coefficient columns registered so unfix triggers correct row reprocessing - Bulk API: chgBounds with interleaved L/U layout via np.repeat/np.tile, chgColType with np.int8 array, addRows/addSets with pre-allocated np.int64 start arrays - Public persistent API delegates to change detector to keep observer callbacks in sync xpress_direct.py: reduced to ~30 lines delegating entirely to build_xpress_problem xpress_direct_base.py: deleted, superseded by xpress_base.py test_xpress_direct.py: - TestXpressPersistentObjective: objective removal between solves (Reason.removed path) - TestXpressPersistentSOS: SOS1 via set_instance + incremental remove via del - load_vars/get_vars/get_reduced_costs subset tests for direct solver loader test_solvers.py: XpressPersistent added to all_solvers and mip_solvers * TXP-8771: Phase A -- eager invalidation, symbolic_solver_labels propagation, write(), multi-active-obj hardening * TXP-8771: add XpressPersistentConfig with auto_updates and wire through ModelChangeDetector * TXP-8771: mutable-param helpers -- remove loadMIP from set_instance, incremental build, chgMCoef batching, minimal helper state * TXP-8771: type annotations, update_parameters signature fix, and minor cleanups across connector * TXP-8771: handle accessors, release/reset, fix public API bugs, fill test coverage to 95% * TXP-8771: fix real LSP errors; suppress false positives from Pyomo untyped API * TXP-8771: fix TC.solverFailure bug (use TC.error), proper cast() for value() return type, str() for Path arg * TXP-8771: add get_xpress_problem() for full Xpress API access * TXP-8772: QP/QCP support and connector refactor Quadratic support: - _quad_arrays() encodes both Xpress scaling conventions (diag_scale=1 for addQMatrix constraints, diag_scale=2 for chgMQObj 0.5*x'Hx objective) - _add_constr_impl calls addQMatrix for QCP constraints after addRows - _set_objective_impl calls chgMQObj for QP objectives - _MutableQuadConstraint: pre-computed handles + 0.5 off-diagonal scaling; rebuild() does delQMatrix+addQMatrix on param change - _MutableObjective extended with quadratic state; diagonal coefs pre-scaled by 2 at construction so update() calls value() unconditionally Compilation: LinearStandardFormCompiler replaced by per-constraint generate_standard_repn(quadratic=True, compute_values=False). Fixed variables become explicit lb=ub columns. EntityMaps.cons holds a single xp.constraint handle per row; range constraints use Xpress 'R'-type rows via get_rhs_and_sense(). Shared _impl statics on XpressSolverMixin: _add_vars_impl, _add_constr_impl, _add_sos_impl, _set_objective_impl, _compile_constraint_repns, _compile_objective_repn (unified NL guard), _xp_obj_sense. Mutable helpers consume repn directly with no intermediate lists; _xp_con passed via constructor. _update_parameters merged to single pass over affected_cons. _populate_results uses lazy-built status lookup dicts. _update_objectives uses get_objective() and skips the active-obj scan when no objective was added. Tests: 18 new tests across QP/QCP, edge cases (empty model, constant objective, range constraints, mutable coef path, fix/unfix rebuild, SOS-only block). Coverage: xpress_direct 100%, xpress_persistent 99%, xpress_base 95%. * TXP-8774: NLP support, unified walker, type annotations, and correctness fixes Core: replace generate_standard_repn-based LP/QP-only path with XpressExpressionWalker (StreamBasedExpressionVisitor) that handles linear, quadratic, and nonlinear expressions uniformly. Variables now registered lazily via _register_variable (addVariable with lb/ub/type in one call), eliminating the pre-scan and duplicate-column bug in SOS handling. Pyomo-to-Xpress maps (_OBJ_SENSE_MAP, _SOL_STATUS_MAP, _STOP_TYPE_MAP, _XP_FUNCTION_MAP, _VAR_XP_TYPE_MAP) consolidated into _init_xpress() called once at module load. NL functions and their decompositions (sinh/cosh/...) merged into a single _XP_FUNCTION_MAP; AbsExpression unified with other unary functions. Exit handlers refactored with _make_binary_const, _make_binary_general, _make_variadic_general, _flag_mutable, _nl_unary to eliminate repetition. Type annotations added throughout (ExpressionBase, UnaryFunctionExpression, _NodeResult alias, visitor: XpressExpressionWalker). Bug fixes: mutable param sign error in NegationExpression contexts (fixed via NegationExpression in _DEPTH_NODES); mutable body constants not tracked for RHS updates (fixed is_mutable detection and _has_mutable_nonlinear flag); _update_parameters single-pass ordering could apply chgMCoef to a stale integer row index after delConstraint shifted rows (fixed with two-pass: NL rebuilds first, linear collects after). Tests: 208 pass. New tests cover NL/quadratic/linear param interactions in all row orderings (including add/remove/re-add cycles), SOS duplicate-column and KeyError bugs, and walker unit coverage for NLP expressions. * TXP-8774: mutable tracking refactor -- _coef_acc/chgRHS/body-const, walker perf, bug fixes, new tests * TXP-8775: register Xpress in NLP/QCP test suite, add has_instance(), add warmstart config flag * TXP-8775: fix mutable quadratic objective with mixed linear coefs; register Xpress in QCP/NLP test suites * TXP-8775: add xpress_persistent to nlp_solvers -- NLP Persistent not available in Gurobi * TXP-8775: coverage 93->97%; new tests for edge cases, has_instance, NLP Persistent registration * TXP-8774: track_mutable=False for XpressDirect -- skip mutable tracking overhead * TXP-8774/8771: quadratic coef tracking, affine expansion, optimizations, coverage, and docs * TXP-8774: walker simplification + ExternalFunction support; mutable tracking moved to persistent xpress_base.py -- removed ~800 lines: - _NodeResult (ExprType compound tuple), WalkerResult dataclass, _coef_acc/_coef_stack mechanism, track_mutable flag, walk_linear, finalizeResult all deleted - all ExprType-keyed exit handlers (_exit_negation/sum/product_*/pow_*/division_*/ unary_const, _define_exit_handlers) deleted - _expand_affine_product, _before_pow, _before_product_bilinear, _extract_affine_form, _flag_mutable, _track_mutable_lin_coef, _track_mutable_const and helpers deleted - ExitNodeDispatcher, ExprType, initialize_exit_node_dispatcher dropped from imports - walker is now a pure xp expression builder; exit dispatch replaced by flat _EXIT_HANDLERS (_ExitHandlerMap dict with MRO fallback for unknown subclasses) - _before_leaf unifies _before_param/_before_npv/_before_native_*; domain errors propagate as ValueError -- no silent nan from try-except in _before_npv - _exit_external_function: ExternalFunctionExpression -> xp.user; supports _fcn/ _grad/_fgh callbacks; unique __name__ per PythonCallbackFunction avoids name clash - SolStatus.FEASIBLE added to _SOL_STATUS_MAP for SLP local-optimum (COMPLETED+FEASIBLE) xpress_persistent.py: mutable tracking moved here, based on generate_standard_repn; _UpdateBatch gains nl_old/new_cons for batched NL del+add; _MutableConstraint NL path xpress_direct.py: API alignment (use_names, plain result object, no track_mutable) tests: 5 new ExternalFunction tests; walker and direct tests updated for new interfaces * TXP-8775: solution pool, test coverage overhaul, hyperbolic decomposition, tolerance tightening to places=6 * TXP-8775: test suite overhaul for the Xpress connector - Add _solve_and_check helper (TypedDict expected, mandatory obj+vars for optimal, nvariables() completeness check) and split shared utilities into _xpress_test_utils.py; split test_xpress_direct.py into separate direct/persistent files - Tighten all assertAlmostEqual tolerances to places=6 (matching Xpress internal precision); redesign degenerate tests with unique optimal solutions - Parametrize 8 structurally identical trig/hyperbolic walker tests into one @parameterized.expand method; add _trivial_model, _solve_lp_no_load, and _solve_check_mutate_check helpers to reduce per-test boilerplate - Reduce .solve( calls from 163 to 34 (79%) by replacing solve+assert sequences with _solve_and_check throughout direct and persistent suites * TXP-8775: add Sphinx documentation for XpressDirect and XpressPersistent connectors * TXP-8775: walker NLP formula assertions, IIS support, pool test split, bug comment fixes - Add NLP formula checks to _walk_and_assert (column indices, scalar CON tokens, operator IFUN/OP tokens via nlpGetFormula); replace product-of-NL*linear tests with product-of-NL*NL to avoid the SLP coefficient-formula storage path where nlpGetFormula returns EOF by design - Add write_iis() and get_iis() to XpressPersistent; write_iis writes the IIS in LP format, get_iis returns conflicting Pyomo ConstraintData and VarData objects mapped via _maps; 3 tests added to TestXpressPersistentIIS - Split TestXpressSolutionPool into test_xpress_pool.py (shared direct/persistent feature, no longer tied to the direct-only test file) - Fix stale comment on test_sos1_vars_not_in_objective (bug fixed by guard in _create_xpress_model; comment was describing the old broken behavior) * TXP-8775: _assemble_xp_expr refactor, pool simplification, stable_xp mutation fix - Extract _assemble_xp_expr helper shared by _make_xp_con, _MutableObjective.update, and _build_xp_from_repn; uses xp.Sum([...]) to avoid aliasing stable_xp - Fix stable_xp mutation bug: xp.expression += float mutates in place for multi-term expressions; _MutableObjective now builds expression without aliasing stable_xp - pool_solutions: drop negative-N rolling-window variant; N>=0 always keeps last N solutions (NonNegativeInt domain); serializepreintsol=1 for determinism - Update config descriptions and test expectations to match new semantics - Fix test_time_limit_zero to assert timelimit instead of removed worklimit branch - Update hyperbolic scalar expectations after /2 -> *0.5 substitution * TXP-8775: update pool_solutions docs to reflect rolling-window-only API * TXP-8775: add symbolic_labels support for SOS set names via addNames Co-authored-by: Claude --- .../reference/topical/solvers/index.rst | 1 + .../reference/topical/solvers/xpress.rst | 234 +++ pyomo/contrib/solver/plugins.py | 12 + .../contrib/solver/solvers/xpress/__init__.py | 11 + .../solver/solvers/xpress/xpress_base.py | 876 ++++++++++ .../solver/solvers/xpress/xpress_direct.py | 122 ++ .../solvers/xpress/xpress_persistent.py | 1006 +++++++++++ .../tests/solvers/_xpress_test_utils.py | 141 ++ .../solver/tests/solvers/test_solvers.py | 13 +- .../tests/solvers/test_xpress_direct.py | 796 +++++++++ .../tests/solvers/test_xpress_persistent.py | 1537 +++++++++++++++++ .../solver/tests/solvers/test_xpress_pool.py | 179 ++ .../tests/solvers/test_xpress_walker.py | 743 ++++++++ 13 files changed, 5669 insertions(+), 2 deletions(-) create mode 100644 doc/OnlineDocs/reference/topical/solvers/xpress.rst create mode 100644 pyomo/contrib/solver/solvers/xpress/__init__.py create mode 100644 pyomo/contrib/solver/solvers/xpress/xpress_base.py create mode 100644 pyomo/contrib/solver/solvers/xpress/xpress_direct.py create mode 100644 pyomo/contrib/solver/solvers/xpress/xpress_persistent.py create mode 100644 pyomo/contrib/solver/tests/solvers/_xpress_test_utils.py create mode 100644 pyomo/contrib/solver/tests/solvers/test_xpress_direct.py create mode 100644 pyomo/contrib/solver/tests/solvers/test_xpress_persistent.py create mode 100644 pyomo/contrib/solver/tests/solvers/test_xpress_pool.py create mode 100644 pyomo/contrib/solver/tests/solvers/test_xpress_walker.py 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..3dfe25f1648 --- /dev/null +++ b/doc/OnlineDocs/reference/topical/solvers/xpress.rst @@ -0,0 +1,234 @@ +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 connectors 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 connector 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 +------------- + +Both connectors accept a common set of configuration options passed as +keyword arguments to :meth:`~XpressDirect.solve`. Options labelled +*framework* are defined in the Pyomo solver framework base class and are +expected to be supported by every compliant connector. Options labelled +*connector* are specific to this Xpress implementation. + +.. list-table:: + :header-rows: 1 + :widths: 25 12 63 + + * - Option + - Scope + - Description + * - ``time_limit`` + - framework + - Wall-clock time limit in seconds. + * - ``threads`` + - framework + - Number of solver threads. + * - ``rel_gap`` + - framework + - Relative MIP optimality gap tolerance. + * - ``abs_gap`` + - framework + - Absolute MIP optimality gap tolerance. + * - ``symbolic_solver_labels`` + - framework + - Use Pyomo component names in the Xpress problem (aids debugging). + * - ``solver_options`` + - framework + - Dict of raw solver control names forwarded directly to the solver. + * - ``warmstart`` + - connector + - Pass variable values as a MIP start hint (default ``True``). + * - ``pool_solutions`` + - connector + - Collect multiple MIP solutions during B&B (0 = disabled). ``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 + +:attr:`~pyomo.contrib.solver.common.results.Results.termination_condition` +reports why the solver stopped; +:attr:`~pyomo.contrib.solver.common.results.Results.solution_status` +reports what was returned. For NLP problems solved via Xpress SLP, +``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 connector ships with a test suite covering LP, MIP, QP, QCP, NLP, +MINLP, SOS, mutable parameter tracking, incremental structural updates, +and the solution pool. + +.. note:: + + Development of this connector was aided by + `Claude Code `_ (Anthropic). + +API Reference +------------- + +.. autosummary:: + :toctree: generated/ + + XpressDirect + XpressPersistent diff --git a/pyomo/contrib/solver/plugins.py b/pyomo/contrib/solver/plugins.py index c16c369abd3..ec962a02fe2 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 @@ -48,3 +50,13 @@ def load(): legacy_name="knitro_direct", doc="Direct interface to KNITRO solver", )(KnitroDirectSolver) + SolverFactory.register( + name="xpress_direct", + legacy_name="xpress_direct_v2", + doc="Direct (scipy-based) 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..5f3f79f6771 --- /dev/null +++ b/pyomo/contrib/solver/solvers/xpress/xpress_base.py @@ -0,0 +1,876 @@ +# ____________________________________________________________________________________ +# +# 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 infrastructure for the Xpress connector: expression walker, solution loaders, +and the constraint/variable build helpers used by both XpressDirect and XpressPersistent. + +Architecture overview +--------------------- +The connector converts a Pyomo ConcreteModel into an xp.problem by walking every +constraint body with XpressExpressionWalker, which turns each Pyomo expression tree +into an equivalent Xpress expression object (xp.Sum, xp.constraint, etc.). + +This per-constraint walk approach handles LP, MIP, QP, QCP, and NLP in a single code +path without special-casing. The linear walk fast path, PauseGC, and single-pass +to_bounded_expression keep per-constraint overhead low enough that the flexibility +comes at no practical cost. + +Mutable parameter tracking (XpressPersistent) +---------------------------------------------- +XpressPersistent uses generate_standard_repn (in xpress_persistent.py) to identify +mutable coefficients at set_instance time. The walker here is a pure expression +builder: it produces xp expressions from Pyomo expression trees but does not track +which params affect which matrix entries. All mutable tracking logic lives in +xpress_persistent.py. +""" + +import datetime +import io +import math +import os +import time +from dataclasses import dataclass +from typing import Any, Iterator, Mapping, Optional, Sequence, cast + +import numpy as np + +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.expr import native_numeric_types +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.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, +) + +# ---- Pure-Python module-level constants (no xp dependency) ------------------ + +_BOUND_TYPE_CODES = np.array([76, 85], dtype=np.int8) # 'L', 'U' + +_VAR_TYPE_CODES: dict[tuple[bool, bool], int] = { + (True, True): 66, # 'B' -- binary (implies integer) + (False, True): 73, # 'I' -- integer + (False, False): 67, # 'C' -- continuous +} + +# ---- xp-dependent maps (empty until _init_xp_maps() is called) -------------- +# Mutable dicts so that `from .xpress_base import _OBJ_SENSE_MAP` in sibling +# modules receives the same object and sees updates via .update(). + +_VAR_XP_TYPE_MAP: dict = {} +_OBJ_SENSE_MAP: dict = {} +_SOL_STATUS_MAP: dict = {} +_STOP_TYPE_MAP: dict = {} +_XP_FUNCTION_MAP: dict = {} +_CON_TYPE_MAP: dict = {} # (is_range, is_equality, has_ub) -> xp constraint type + + +class _ExitHandlerMap(dict): + """Dict with MRO fallback: unknown subclasses resolve via their base class. + + Pyomo has many concrete expression subclasses (e.g. ScalarExpression) that + are not registered explicitly. Walking the Python MRO finds the registered + base class handler (e.g. NamedExpressionComponent for all named expressions) + and caches the result so subsequent lookups are direct dict hits. + Only the minimal set of base classes needs to be registered. + """ + + 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 nodes via xp.user. + + `data` is (xp_arg1, ..., xp_argN, fcn_id_int). The last element is the + evaluated _PythonCallbackFunctionID (a plain Python int, not an xp expression). + Only PythonCallbackFunction is supported; AMPLExternalFunction has no + callable Python callbacks. + """ + pyo_fcn = node._fcn + if not isinstance(pyo_fcn, PythonCallbackFunction): + raise IncompatibleModelError( + f"ExternalFunction of type '{type(pyo_fcn).__name__}' uses AMPL callbacks " + "which are not supported by the Xpress connector. " + "Use PythonCallbackFunction instead." + ) + xp_args = data[:-1] # strip fcn_id (last element is always a plain int) + has_grad = pyo_fcn._grad is not None or pyo_fcn._fgh is not None + + if pyo_fcn._fgh is not None: + _fgh = pyo_fcn._fgh + + def xp_cb(*vals): + f, g, _ = _fgh(list(vals), 1, None) + return (f, *g) if g is not None else f + + elif pyo_fcn._grad is not None: + _fcn, _grad = pyo_fcn._fcn, pyo_fcn._grad + + def xp_cb(*vals): + args = list(vals) + return (_fcn(*args), *_grad(args, None)) + + else: + _fcn = pyo_fcn._fcn + + def xp_cb(*vals): + return _fcn(*vals) + + # Xpress identifies user functions by their __name__; each distinct Pyomo + # ExternalFunction needs a unique name to avoid a "duplicate user function" error. + 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() -> tuple: + """Populate all Pyomo<->Xpress entity maps. No-op if already done or if + xpress did not import successfully.""" + + xp, xpress_available = attempt_import('xpress', catch_exceptions=(Exception,)) + + if not xpress_available: + return xp, xpress_available + + _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), + # SLP (xp.user) reports FEASIBLE on completion: local optimum found. + 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, + } + ) + + # 1:1 mapping from pyomo nodes to Xpress operators + _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, + } + ) + + # Pyomo defines many unary operators with a single class + a name attribute. + # We need a 2 levels dispatch with this second dict. + _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, + # Decomposed via xpress primitives -- same call signature as direct fns + '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 + } + ) + + return xp, xpress_available + + +xp, xpress_available = _init_xpress() + + +def _register_pool_collector(prob, pool_limit: int) -> 'list[list[float]]': + """Register an IntsolCallback to collect MIP solutions during solve. + + Returns the list that will be populated by the callback. Solution 0 + (the incumbent) is always available via prob.getSolution(); pool[k] holds + the (k+1)-th solution kept in the rolling window. + + pool_limit > 0: keep a rolling window of the last N solutions found. + When the window is full the oldest entry is evicted on each new solution, + so the pool always contains the N most recently found feasible solutions. + + Zero overhead for LP/QP: the callback never fires for continuous problems. + """ + pool: list[list[float]] = [] + + def _intsol_cb(cbprob, cbdata): + if len(pool) == pool_limit: + pool.pop(0) + pool.append(cbprob.getCallbackSolution()) + + # We set SERIALIZEPREINTSOL to ensure determinism in how solutions are found + prob.controls.serializepreintsol = 1 + prob.addIntsolCallback(_intsol_cb) + return pool + + +@dataclass +class EntityMaps: + """Stable handle maps from Pyomo entities to Xpress entity objects. + + Xpress entity handles (xp.var, xp.constraint, xp.sos) remain valid after + other entities are deleted, unlike integer indices, which Xpress renumbers + after every deletion. Storing handles here eliminates all index rebuilds. + + vars: keyed by id(VarData), not VarData itself, because VarData objects are + not hashable in all Pyomo versions and id() is stable for the lifetime of + the solve session. + + cons: each constraint maps to a single xp.constraint handle. Range constraints + are stored as Xpress 'R'-type rows, one row per constraint. + + sos: values are single xp.sos handles (SOS sets are never split). + """ + + vars: dict[int, Any] # id(VarData) -> xp.var + cons: dict[ConstraintData, Any] # ConstraintData -> xp.constraint + sos: dict[SOSConstraintData, Any] # SOSConstraintData -> xp.sos + + +class XpressSolutionLoaderBase(SolutionLoader): + """Base solution loader shared by direct and persistent Xpress solvers.""" + + 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 + # list[list[float]]: pool[k] = solution (k+1) in B&B discovery order. + # Populated by _register_pool_collector during prob.optimize(). + self._pool: list = pool_solutions if pool_solutions is not None else [] + # Active solution id: 0 = incumbent (default). k > 0 = pool[k-1]. + # _set_solution_id(None) is treated the same as 0 for framework compat. + self._active_id: int = 0 + + def _query_vars( + self, variables: Optional[Sequence[VarData]], fn, exc_type: type[Exception] + ) -> Iterator[tuple[VarData, float]]: + """Query a variable-valued attribute from the Xpress problem, returning + (VarData, value) pairs. Raises exc_type on xp.ModelError.""" + 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: + # LP/QP: pool stays empty, so returns 1 (the incumbent). + # MIP with pool_solutions > 0: incumbent + len(pool) collected solutions. + 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 currently active solution. + + Active solution 0 or None -> incumbent via prob.getSolution(). + Active solution k > 0 -> pool[k-1] (B&B discovery order). + + For pool entries, values are indexed by xp.var.index (Xpress column + index), which matches the column registration order in self._vars. + """ + 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: + # Note: Xpress cons order might differ from maps.cons dict order + # thus, prob.getDuals() with no args is not viable here + 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 logic for XpressDirect and XpressPersistent.""" + + _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: + 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: + """Bulk-set column types (C/I/B) for pyo_vars.""" + 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: + """Bulk-set variable bounds for pyo_vars, respecting fixed-var pinning.""" + n = len(pyo_vars) + cbounds = [b for var in pyo_vars for b in self._var_bounds(var)] + prob.chgBounds(np.repeat(xp_vars, 2), np.tile(_BOUND_TYPE_CODES, n), cbounds) + + def _add_vars_impl(self, prob, pyo_vars: list[VarData], symbolic_labels: bool): + """Add columns, set types and bounds. Returns the xp.var array. + + Uses addVariables(n, name='') to get sequential C1/C2/... auto-names + that never repeat across incremental calls. Without name='', Xpress + generates x(0)/x(1)/... which conflict on a second batch. + """ + 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 pyo_cons and build xp.constraint objects in a single tight loop. + + Performance choices: + - PauseGC wraps the loop: each walk creates many short-lived objects; + deferring GC eliminates unpredictable mid-loop pauses. + - to_bounded_expression called once per constraint to obtain (lb, body, ub) + together, avoiding three separate property accesses. + - _before_linear fast path: LinearExpression bodies (the common case for LP/MIP) + bypass the full StreamBasedExpressionVisitor dispatch, saving + initializeWalker + enterNode overhead. + """ + 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. Returns the xp.sos handle array.""" + n = len(pyo_sos) + if n == 0: + return [] + settype = np.empty(n, dtype=np.int8) + setstart = np.empty(n + 1, dtype=np.int64) + setstart[0] = 0 + setind: list[int] = [] + refval: list[float] = [] + for i, con in enumerate(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[i] = ord('1' if con.level == 1 else '2') + setstart[i + 1] = 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: + n = len(entind) + ws_vals = np.empty(n, dtype=np.float64) + ws_cols = np.empty(n, dtype=np.int32) + count = 0 + for j in entind: + var = vars[j] + if var.value is not None: + ws_vals[count] = var.value + ws_cols[count] = j + count += 1 + if count > 0: + prob.addMipSol(ws_vals[:count], ws_cols[:count]) + + 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): + """Return the xp.var for pyo_var, registering it in prob/var_map on first use.""" + 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): + # child.args bypasses LinearExpression._build_cache. + # xp.Sum(list) is a single C call; Python + accumulation is O(n) calls. + 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. + + Pure expression builder: no ExprType tracking, no mutable param tracking. + Each node type maps directly to its Xpress operation in _EXIT_HANDLERS. + Mutable tracking is handled by generate_standard_repn in xpress_persistent.py. + + _before_linear provides a fast path for LinearExpression bodies (the dominant + case in LP/MIP models) via direct call from _add_cons_impl. + """ + + 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..169b0baa7c1 --- /dev/null +++ b/pyomo/contrib/solver/solvers/xpress/xpress_direct.py @@ -0,0 +1,122 @@ +# ____________________________________________________________________________________ +# +# 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.config import ConfigValue, NonNegativeInt +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, + _OBJ_SENSE_MAP, + _register_variable, + _register_pool_collector, + xp, +) + + +class XpressDirectSolutionLoader(XpressSolutionLoaderBase): + """Solution loader for the non-persistent XpressDirect solver.""" + + +class XpressDirectConfig(BranchAndBoundConfig): + """Configuration for XpressDirect.""" + + 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.' + ), + ), + ) + + +class XpressDirect(XpressSolverMixin): + CONFIG = XpressDirectConfig() + + 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) + # SOS vars may not appear in constraints/objective -- register if missing. + 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) + if 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..34cefb42ab4 --- /dev/null +++ b/pyomo/contrib/solver/solvers/xpress/xpress_persistent.py @@ -0,0 +1,1006 @@ +# ____________________________________________________________________________________ +# +# 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.common.config import ConfigValue, NonNegativeInt +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, + _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 for a single flush. + + LP/QP affine updates are applied first (no structural changes, indices stable). + NL rebuilds follow as a single batched delConstraint + addConstraint pair. + The new xp.constraint objects are built during collect() so flush() needs + only prob, maps, and mutable_helpers. + """ + + def __init__(self): + # parallel lists for chgMCoef(rows, cols, vals) + self.coef_rows = [] + self.coef_cols = [] + self.coef_vals = [] + # parallel lists for chgRHS(rows, vals) + self.rhs_rows = [] + self.rhs_vals = [] + # parallel lists for chgRHSRange(rows, vals) + self.rng_rows = [] + self.rng_vals = [] + # list of (row, col1, col2, val) tuples for chgQRowCoeff + self.quad_updates = [] + # NL row replacement: old handles for delConstraint, new for addConstraint + 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: + """Handles mutable-param updates for a constraint row. + + LP/QP path (nl_expr is None): targeted chgMCoef/chgRHS/chgQRowCoeff updates. + NL path (nl_expr set): full row rebuild -- re-walks nl_expr only; reconstructs + the rest from a precomputed stable xp expression and parallel coef lists. + + generate_standard_repn guarantees unique variable entries so parallel lists + suffice -- no dict accumulation needed. + + _rhs_expr: bound - body_constant (Pyomo expr). For NL, always set (even when + non-mutable) since it is needed for the rebuild's xp.constraint call. + _rng_expr: ub - lb for range constraints (Pyomo expr). Same rule as _rhs_expr. + _type: xp constraint type (leq/geq/eq/rng), stored for NL rebuild only. + _stable_xp: precomputed xp expression for the non-mutable lin+quad body terms. + None for LP/QP (not needed for targeted updates). + """ + + __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: + # ConstraintData; needed to update maps.cons after NL rebuild + self._con = con + # xp.constraint handle; updated in-place after each NL rebuild + self._xp_con = xp_con + # mutable linear xp.var handles (parallel with _lin_coefs) + self._lin_vars = lin_vars + # mutable linear Pyomo coef expressions + self._lin_coefs = lin_coefs + # mutable quad var handles, parallel triplet (_quad_v1s, _quad_v2s, _quad_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 + # repn.nonlinear_expr; None for LP/QP, set for NL constraints + 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 a single Xpress expression from precomputed parts. + + Coefs may be Pyomo expressions or pre-evaluated floats; value() is called on + each one. All vars must already be xp.var handles. + stable_xp is a precomputed xp expression (or None if absent). + constant is a Pyomo expression, float, or None; value() is called if not None. + nl_expr is a walked xp expression (or None). + + Uses xp.Sum([...]) so stable_xp is never aliased or mutated. + """ + 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: + """Handles mutable-param updates for the objective function. + + Mirrors _MutableConstraint: both LP/QP and NL paths share parallel coef lists. + NL additionally carries _stable_xp and _nl_expr (re-walked on every update). + + LP/QP: chgObj(keys+[-1], vals+[const]) + chgMQObj (Hessian-scaled at update time). + NL: setObjective(stable_xp + mutable_terms + const + nl_part). + """ + + __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: + # mutable linear xp.var handles and Pyomo coef expressions (parallel) + self._lin_vars = lin_vars + self._lin_coefs = lin_coefs + # mutable quad var handles and coefs, parallel triplet + self._quad_v1s = quad_v1s + self._quad_v2s = quad_v2s + self._quad_coefs = quad_coefs + # LP/QP: mutable Pyomo expr or None; NL: repn.constant always (may be 0) + self._constant = constant + # non-mutable lin+quad body as a precomputed xp expr; None for LP/QP + self._stable_xp = stable_xp + # repn.nonlinear_expr; None for LP/QP objectives + self._nl_expr = nl_expr + + def update(self, prob, walker: 'XpressExpressionWalker') -> None: + if self._nl_expr is None: + if self._constant is not None: + 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 stores the QP objective as (1/2)*x'Qx, so chgMQObj + # expects Hessian-scaled values. Off-diagonal entries appear twice + # in the symmetric Hessian (which compensates), diagonal entries + # must be doubled to match the user-visible coefficient. + 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 for XpressPersistent -- invalidated before each 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(BranchAndBoundConfig): + """Configuration for XpressPersistent.""" + + 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.auto_updates = self.declare('auto_updates', AutoUpdateConfig()) + self.warmstart: bool = self.declare( + 'warmstart', + ConfigValue( + default=True, + domain=bool, + description='Pass current integer variable values as a MIP warm start.', + ), + ) + 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.' + ), + ), + ) + + +class XpressPersistent(XpressSolverMixin, PersistentSolverBase, Observer): + """Persistent Xpress solver: reuses the xp.problem() across re-solves.""" + + CONFIG = XpressPersistentConfig() + + def __init__(self, **kwds): + super().__init__(**kwds) + # active xp.problem(); None until set_instance + self._xp_prob = None + # model currently loaded; identity-checked on re-solve + self._pyomo_model: Optional[BlockData] = None + # ordered VarData list; used for batch unfix before repn calls + self._vars: Optional[list[VarData]] = None + # Pyomo->Xpress handle maps (vars, cons, sos) + self._maps: Optional[EntityMaps] = None + # active Pyomo objective; None if model has no objective + self._objective: Optional[ObjectiveData] = None + # invalidated on every model change; prevents stale solution reads + self._last_solution_loader: Optional[XpressPersistentSolutionLoader] = None + # watches model for changes and notifies this observer + self._change_detector: Optional[ModelChangeDetector] = None + # per-constraint update helpers (LP/QP targeted or NL rebuild) + self._mutable_helpers: dict[ConstraintData, _MutableConstraint] = {} + # None when objective has no mutable params + self._mutable_objective: Optional[_MutableObjective] = None + # drives symbolic names in Xpress (symbolic_solver_labels config) + self._use_names: bool = False + # lazily initialised on first NL walk; reused across all walks + 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 an Xpress expression from a StandardRepn (linear + quadratic parts).""" + 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, descend_into=True, active=True) + ) + old_sos = list( + block.component_data_objects(SOSConstraint, descend_into=True, 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: + """Return True if set_instance has been called and a model is loaded.""" + return self._pyomo_model is not None + + def get_xpress_problem(self) -> Any: + """Return the underlying xp.problem object 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 the 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 the 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 the xp.sos handle for a Pyomo SOS constraint.""" + assert self._maps is not None + return self._maps.sos[con] + + def release(self) -> None: + """Drop the xp.problem and release all solver resources.""" + self._clear() + + def reset(self) -> None: + """Clear all model data from the 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 the loaded Xpress problem to a file.""" + self._xp_prob.writeProb(filename, flags) + + def write_iis(self, filename: str) -> str: + """Compute the IIS and write it to filename in LP format. + + Must be called after an infeasible solve. Raises if no model is + loaded or if Xpress cannot compute an IIS. + + Returns the filename written. + """ + 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 the IIS and return the conflicting Pyomo objects. + + Must be called after an infeasible solve. Raises if no model is + loaded or if Xpress cannot compute an IIS. + + Returns a dict with keys: + 'constraints': list[ConstraintData] -- Pyomo constraints in the IIS + 'variables': list[VarData] -- Pyomo variables whose bounds + contribute to the IIS + """ + 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..f52ec46940c --- /dev/null +++ b/pyomo/contrib/solver/tests/solvers/_xpress_test_utils.py @@ -0,0 +1,141 @@ +# ____________________________________________________________________________________ +# +# 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 connector test suite. + +Imported by test_xpress_direct.py and test_xpress_persistent.py. +Not a test module itself (underscore prefix prevents pytest collection). +""" + +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, enforce baseline health checks, return the result. + + Always asserts termination_condition and solution_status. + When status resolves to optimal, 'objective' and 'vars' are mandatory + in expected -- a test that solves optimally without checking the solution + is not a useful test. + + **solve_kwargs: forwarded verbatim to opt.solve(). + """ + 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(): + """Single bounded variable, no constraints, minimise x. + + Used by persistent API surface tests that need a live xp.problem but do not + care about the specific solution (handles, controls, state checks, etc.). + """ + 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): + """Return (model, result) for _simple_lp() solved with load_solutions=False. + + Avoids the two-line preamble that is repeated in every test that exercises + the solution-loader interface (get_vars, get_duals, get_reduced_costs, etc.). + """ + 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, +): + """Two-step mutable param test: solve + check, mutate param, solve + check again. + + Covers the pervasive pattern in persistent tests where one param change is + applied between two consecutive solves and both solutions are verified. + Returns (res_before, res_after). + + Only use when there are no assertions that compare values across the two + solves (e.g. assertLess(x2, x1)). Those tests must keep explicit solve calls. + """ + 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 686c45adb74..d5e63897763 100644 --- a/pyomo/contrib/solver/tests/solvers/test_solvers.py +++ b/pyomo/contrib/solver/tests/solvers/test_solvers.py @@ -42,6 +42,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 @@ -76,6 +77,8 @@ def param_as_standalone_func(cls, p, func, name): ('highs', Highs), ('gams', GAMS), ('knitro_direct', KnitroDirectSolver), + ('xpress_direct', XpressDirect), + ('xpress_persistent', XpressPersistent), ] mip_solvers = [ ('gurobi_persistent', GurobiPersistent), @@ -83,17 +86,23 @@ def param_as_standalone_func(cls, p, func, name): ('gurobi_direct_minlp', GurobiDirectMINLP), ('highs', Highs), ('knitro_direct', KnitroDirectSolver), + ('xpress_direct', XpressDirect), + ('xpress_persistent', XpressPersistent), ] nlp_solvers = [ ('gurobi_direct_minlp', GurobiDirectMINLP), ('ipopt', Ipopt), ('knitro_direct', KnitroDirectSolver), + ('xpress_direct', XpressDirect), + ('xpress_persistent', XpressPersistent), ] qcp_solvers = [ ('gurobi_persistent', GurobiPersistent), ('gurobi_direct_minlp', GurobiDirectMINLP), ('ipopt', Ipopt), ('knitro_direct', KnitroDirectSolver), + ('xpress_direct', XpressDirect), + ('xpress_persistent', XpressPersistent), ] qp_solvers = qcp_solvers + [("highs", Highs)] miqcqp_solvers = [ @@ -1688,8 +1697,8 @@ def test_log(self, name: str, opt_class: Type[SolverBase], use_presolve: bool): m.obj = pyo.Objective(expr=m.x**2 + m.y**2) m.c1 = pyo.Constraint(expr=m.y <= pyo.log(m.x)) res = opt.solve(m) - self.assertAlmostEqual(m.x.value, 0.6529186341994245) - self.assertAlmostEqual(m.y.value, -0.42630274815985264) + self.assertAlmostEqual(m.x.value, 0.6529186341994245, places=5) + self.assertAlmostEqual(m.y.value, -0.42630274815985264, places=5) @mark_parameterized.expand(input=_load_tests(all_solvers)) def test_with_numpy( 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..813731a2cd6 --- /dev/null +++ b/pyomo/contrib/solver/tests/solvers/test_xpress_direct.py @@ -0,0 +1,796 @@ +# ____________________________________________________________________________________ +# +# 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 +from unittest.mock import MagicMock + +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 +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, + ) + # Verify that Pyomo names appear in the LP file -- symbolic_solver_labels + # must not be silently dropped. + 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}, + ) + # Invalid control names are forwarded to Xpress and raise -- proves + # solver_options are not silently ignored. + with self.assertRaises(Exception): + self.opt.solve(m, solver_options={'_invalid_control_xyz': 1}) + + def test_rel_gap(self): + # Verify rel_gap is accepted and the solve completes without error. + # A trivial MIP solves to optimality regardless of gap tolerance. + 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() + # Provide a valid feasible hint: x=0, y=4 satisfies all constraints + # and is the optimal solution. Verify that the hint is consumed (at least + # one MIP integer solution found, which includes the warm-start point). + 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) + # LP always runs at least one iteration + self.assertGreaterEqual(res.extra_info.simplex_iterations, 1) + self.assertGreaterEqual(res.extra_info.barrier_iterations, 0) + # LP has no B&B nodes + self.assertEqual(res.extra_info.node_count, 0) + # LP has no MIP solutions + 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) # untouched + 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) + # two variables -- both orders must produce correct mapping + 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): + # y is in-basis at the LP optimum, so its reduced cost is 0. + 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): + # get_vars() with no argument exercises the full-variable path in + # _get_solution_vals (vars_to_load=None -> return all self._vars). + 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) + # SOS1 optimal: only x[3] can be nonzero (highest weight/coefficient), x[3]=1 -> obj=3.0 + _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): + # SOS vars that do NOT appear in any constraint or objective. + 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) # x vars intentionally absent + _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): + # SOS vars that also appear in the objective must NOT get duplicate + # Xpress columns -- the problem should have exactly 3 columns (x[1..3]). + 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): + # Xpress may return a scalar (not a list) when queried for a single + # constraint dual -- the scalar-to-list normalization must be exercised. + 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): + # _simple_lp() optimal: x=0 (non-basic at lb), y=4 (in-basis). + # Analytically: RC[x] = c_x - u1 * a_{x,c1} = -1 - (-2)*1 = 1.0 exactly. + # RC[y] = 0 (y is basic in the LP). + 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): + # c1 (x+y<=4) is binding at the optimum. + # Shadow price: relaxing c1 by 1 allows y to increase by 1, improving obj by -2. + # Xpress convention: dual[c1] = -2.0 for this minimization problem. + # c2 (2x+y<=6) has slack 2 at the optimum -> dual = 0. + 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): + # min x^2 + y^2 s.t. x + y >= 1 -> optimal x=y=0.5, obj=0.5 + 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): + # max x + y s.t. x^2 + y^2 <= 1, x >= 0, y >= 0 -> optimal at x=y=1/sqrt(2) + 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)) # maximize x+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): + # x**3 >= 1, min x: optimal solution is x=1 via Xpress NLP solver. + 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): + """Edge cases and rarely-hit branches for the direct connector.""" + + def setUp(self): + self.opt = XpressDirect() + + def test_nl_cubic_objective_direct(self): + # min x**3 for x in [0, 1]: optimal is x=0 (boundary minimum). + 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): + # Verifies abs_gap config is forwarded to mipabsstop. Trivial MIP solves + # to optimality regardless of gap tolerance. + 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): + # working_dir must chdir into the directory before optimize and restore + # the original cwd after, even if optimize raises. + 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): + # No constraints at all: only bounds and objective. Optimal at lower bound. + 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): + # No objective: Xpress treats it as a feasibility problem. The reported + # incumbent_objective should be None (has_obj=False path). + 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): + # Objective with no variable terms (constant only). Exercises the + # len(xp_vars) == 0 branch in _set_objective_impl. + 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) + # x value is solver-determined (feasibility only); Xpress returns x=1.0 (lb of constraint) + _solve_and_check(self, self.opt, m, {'objective': 5.0, 'vars': [(m.x, 1.0)]}) + + def test_range_constraint_lp(self): + # Range constraint: 1 <= x + y <= 3. Exercises the 'R' rowtype path + # in get_rhs_and_sense (rng_arr branch). + 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)) + # obj=-2x-y: unique optimal at upper bound is x=3, y=0 (x has larger coefficient) + 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): + # Default cons_to_load=None path exercises maps.cons.values() ordering. + 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): + # var.fix() with no argument sets fixed=True but leaves value=None. + # _var_bounds must raise a descriptive error rather than the opaque + # TypeError: float() argument must be a real number, not NoneType. + m = pyo.ConcreteModel() + m.x = pyo.Var() + m.x.fix() # fixed=True, value stays None + 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): + # Unit test for _apply_solver_controls: verify controls are actually written + # to xp.problem. A bug that silently drops the call would leave all four + # solver-config tests green while this one fails. + + 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): + # Model with lb > ub (inverted bounds) is detected as infeasible by Xpress. + # The connector must return provenInfeasible without raising. + m = pyo.ConcreteModel() + m.x = pyo.Var(bounds=(5, 3)) # inverted: lb=5 > ub=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): + # time_limit=0 passes 0.0 directly to timelimit (Xpress interprets 0 as no limit). + 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): + # The finally block in solve() restores cwd even when an exception propagates + # (e.g., from an invalid solver option -- distinct from the InfeasibleConstraintException + # catch which is handled separately). + 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 # expected: invalid control raises inside the try block + self.assertEqual(os.getcwd(), original_cwd) + + +@unittest.pytest.mark.solver('xpress_persistent') +class TestXpressDirectNLP(unittest.TestCase): + """NLP integration tests for the direct connector.""" + + def setUp(self): + self.opt = XpressDirect() + + def _check_optimal(self, res): + self.assertEqual( + res.termination_condition, TerminationCondition.convergenceCriteriaSatisfied + ) + + def test_nl_exp_objective_linear_constraints(self): + # min exp(x) s.t. x >= 1, x in [0,3] + 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): + # min x+y s.t. sin(x) + y <= 1, x in [0, pi/2], y in [0,1] + 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): + # min sin(x) s.t. exp(x) <= 2, x in [0, 2] + 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)) + # exp(x) <= 2 -> x <= ln(2); minimizing sin(x) on [0, ln(2)] -> x=0 + _solve_and_check(self, self.opt, m, {'objective': 0.0, 'vars': [(m.x, 0.0)]}) + + def test_nl_range_constraint(self): + # 0.5 <= sin(x) <= 1, min x for x in [0, pi] + 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): + # sin(x_fixed) + y <= 5, min y; with x fixed at pi/2, sin(x)=1. + 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): + # min abs(x - 2), x in [0, 5]: optimal x=2, obj=0 + 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): + """Integration tests for ExternalFunction support via xp.user (SLP).""" + + def setUp(self): + self.opt = XpressDirect() + + def _check_solved(self, res): + # xp.user (SLP) reports SolutionStatus.feasible on local-optimum convergence, + # not SS.optimal (which requires global-optimality proof). + self.assertIn( + res.solution_status, (SolutionStatus.optimal, SolutionStatus.feasible) + ) + self.assertIsNotNone(res.incumbent_objective) + + def test_external_function_no_gradient(self): + """ExternalFunction without gradient: model builds and solves correctly. + + min f(x,y) = x^2 + y s.t. x in [0,5], y in [1,5]. + Optimum: x=0, y=1, obj=1. + """ + 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): + """ExternalFunction with gradient: derivatives='always' path; solves correctly. + + min f(x,y) = x^2 + y s.t. x in [0,5], y in [1,5]. + Optimum: x=0, y=1, obj=1. + """ + + 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_external_function_ampl_raises(self): + """Non-PythonCallbackFunction must raise IncompatibleModelError. + + Calls _exit_external_function directly with a mock to avoid needing a + working AMPLExternalFunction (which requires a shared library on disk). + """ + node = MagicMock() + node._fcn = MagicMock() # not a PythonCallbackFunction instance + with self.assertRaises(IncompatibleModelError): + _exit_external_function(None, node) + + def test_external_function_in_constraint(self): + """External function in a constraint: constraint is respected at optimum. + + g(y) = y (identity), constraint g(y) >= 1, objective min y, y in [0,5]. + Optimum: y=1, objective=1. The external function identity forces y>=1. + """ + + 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): + """Two external functions in the same model. + + Objective: f1(x) + f2(y) = x^2 + (y+1), x in [0,5], y in [0,5]. + Optimum: x=0, y=0, obj=1. + """ + + 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)) + _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): + """ExternalFunction with fgh= callback: exercises the _fgh code path in + _exit_external_function (xpress_base.py lines 176-180). + + fgh(args, fgh_flag, fixed) returns (f, g, None): + f = x^2 + y (function value) + g = [2x, 1.0] (gradient, computed when fgh_flag != 0) + + min x^2 + y s.t. x in [0,5], y in [1,5] -> optimal x=0, y=1, obj=1. + """ + + 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..3e11cdde804 --- /dev/null +++ b/pyomo/contrib/solver/tests/solvers/test_xpress_persistent.py @@ -0,0 +1,1537 @@ +# ____________________________________________________________________________________ +# +# 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): + # Exercises the Reason.removed path in _update_objectives, which calls + # _set_objective(None). Not covered by any parameterized test because + # those use set_instance (new model), not incremental update. + 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): + # Pyomo idiom: define multiple Objective components as alternatives, + # toggle which one is active via activate()/deactivate(). The change + # detector treats deactivate as Reason.removed and activate as + # Reason.added; our handler clears the old and sets the new. + 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): + # Solve once with one active obj, then activate a second + # without deactivating the first. The detector fires _update_objectives + # with only the new obj added; the multi-active hardening must catch + # the model state and raise. + 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): + # After a two-active-obj raise, the solver must be recoverable via + # set_instance. The IncompatibleModelError leaves the Xpress problem in + # a partially-updated state (obj2 was written before the raise). Calling + # set_instance rebuilds from scratch and resolves cleanly. + 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) + + # Recovery: deactivate the conflicting objective and reset via set_instance. + 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): + # After solve, the loader works. After any mutation, the loader must + # be eagerly invalidated -- subsequent reads raise NoSolutionError. + m = _simple_lp() + res = _solve_and_check( + self, self.opt, m, {'objective': -8.0, 'vars': [(m.x, 0.0), (m.y, 4.0)]} + ) + # baseline: loader works right after solve + res.solution_loader.get_vars() + # mutate: add a constraint + 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): + # With symbolic_solver_labels=True the user's Pyomo names should + # appear in the LP-format output; without it they would not. + 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): + # auto_updates.update_parameters=False tells the detector to skip + # parameter change tracking. Changing a mutable param between solves + # should then have no effect on the loaded problem. + 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 + # Parameter change was ignored: optimum stays at 5.0, not 7.0 + _solve_and_check( + self, + self.opt, + m, + {'objective': -5.0, 'vars': [(m.x, 5.0)]}, + auto_updates={'update_parameters': False}, + ) + + # Positive control: with default (full) auto_updates the param change IS + # picked up -- distinguishes "tracking correctly disabled" from + # "_update_parameters always broken". + _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): + # warmstart=False: the solve must complete correctly even with stale + # variable values set on the model (the hint is not passed to the solver). + # This test verifies the config branch does not crash or corrupt the solve; + # it cannot verify that addMipSol was skipped without a mock. + 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): + # SOS1: at most one variable nonzero. Maximize x1+2*x2+3*x3 s.t. xi in [0,1]. + # With SOS1: optimal is x3=1, others=0 (obj=3). + # After removing SOS1: all vars reach upper bound (obj=6). + 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)]}) + # Deactivate old obj in Pyomo model before adding new one -- the + # _update_objectives active-objective scan would raise on two active objs. + 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())) + # Deactivate in Pyomo model so auto-update does not re-add it on next solve. + 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) + # Deactivate block in Pyomo model so auto-update does not re-add it on next solve. + 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) + # Full Xpress API accessible: use prob + entity handles to read slack + 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): + # min x^2 + y^2 s.t. x + y >= 1 -> optimal x=y=0.5, obj=0.5 + 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)]} + ) + # re-solve (tests persistent re-use) + _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): + # Without QC: min -(x+y) s.t. x,y in [0,2] -> optimal at x=y=2. + # Add QC x^2+y^2<=1 -> optimal at x=y=1/sqrt(2). + # Remove QC -> optimal returns to x=y=2. + 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, + }, + ) + + # Remove QC: optimal returns to x=y=2 (no longer constrained by unit circle). + 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): + # min p*x^2 s.t. x >= 1. Change p, verify chgMQObj is called. + 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) + # chgMQObj updates the quadratic obj coefficient in-place; accuracy ~1e-5 + _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): + # min -(x+y) s.t. x^2 + p*y^2 <= 1. Change p, verify rebuild. + 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, + { + # QCP solver for quadratic constraints achieves ~6e-7 on the circle + 'objective': -math.sqrt(2), + 'vars': [(m.x, math.sqrt(2) / 2), (m.y, math.sqrt(2) / 2)], + 'obj_places': 5, + 'var_places': 5, + }, + ) + # p=4: ellipse x^2 + 4*y^2 <= 1. Lagrange: x=2/sqrt(5), y=1/(2*sqrt(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): + # Same math as test_mutable_param_in_quadratic_constraint but the + # expression is written as (m.p * m.y) * m.y (MonomialTerm * Var path) + # to exercise Case 2 of _before_product_bilinear. + 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, + }, + ) + # p=4: ellipse x^2 + 4*y^2 <= 1. Lagrange: x=2/sqrt(5), y=1/(2*sqrt(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): + # Regression: p1*(x-1)^2 + p2*(y-6)^2 - p3*y has mutable quadratic coefs + # (p1, p2) AND a mutable linear coef (p3). The linear term populated + # _mutable_lin_vars which prevented the fallback from setting + # _has_mutable_nl_formula, so the objective was wrongly tracked as + # _MutableObjective (chgObj only) instead of requiring a full rebuild. + 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)]} + ) + # Change p2: new unconstrained min at (1, 8), constrained at x=y=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): + # min -(x+y) s.t. x^2 + p*y^2 <= 1. + # p=4: ellipse x^2+4y^2<=1; by Lagrange multipliers: + # max x+y on ellipse -> x=2/sqrt(5) ~ 0.8944, y=1/(2*sqrt(5)) ~ 0.2236. + 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, + { + # p=1: circle, QCP ~6e-7 + '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) + # Analytic: x=2/sqrt(5), y=1/(2*sqrt(5)) + 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): + # x**3 >= 1, min x: optimal solution is x=1 via Xpress NLP solver. + 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): + # min x**3 for x in [0, 1]: optimal is x=0 (boundary minimum). + 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): + """Edge cases and rarely-hit branches for the persistent connector.""" + + 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): + # has_instance() returns False before set_instance and True after (L565). + 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): + # _add_vars_impl early-return when list is empty (line 531 in xpress_base.py). + 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): + # _add_constraints early-return when list is empty. + m = _trivial_model() + self.opt.set_instance(m) + self.opt.add_constraints([]) # public API path (via change detector) + self.opt._add_constraints( + [] + ) # direct path -- exercises the len==0 early return + + def test_add_block_sos_only(self): + # add_block path with only SOS constraints (no linear constraints). + # Exercises xpress_persistent.py:_add_block SOS branch. + 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): + # remove_block path with only SOS constraints. + 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): + # Toggling sense without touching expr exercises the Reason.sense-only + # branch in _update_objectives. + 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): + # Persistent path with constant-only objective. + 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) + # x value is solver-determined (feasibility only); Xpress returns x=1.0 (lb of constraint) + _solve_and_check(self, self.opt, m, {'objective': 7.0, 'vars': [(m.x, 1.0)]}) + + def test_range_constraint_persistent(self): + # Range constraint on persistent solver. + 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)) + # obj=-2x-y: unique optimal at upper bound is x=3, y=0 (x has larger coefficient) + 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): + # Range constraint with mutable lb. Exercises the lb_h branch in + # _build_constraint_helper for rowtype == 'R'. + 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): + # Symmetric to test_mutable_param_in_range_constraint (which tests mutable lb). + # This tests mutable ub: 1 <= x <= p. Exercises the ub_h branch in + # _build_constraint_helper for rowtype == 'R'. + 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) # maximize 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): + # The existing test_range_constraint_lp maximizes (upper bound is binding). + # This test minimizes so the LOWER bound is binding. + 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)) + # obj=x+2y: unique optimal at lower bound is x=1, y=0 (y has larger coeff, min -> y=0) + 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): + # Coefficient updated via _MutableConstraint.collect path (chgMCoef). + 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): + # Fixing a variable sets its Xpress column bounds to [val, val] via + # chgBounds; Xpress handles the fixed variable natively without a + # constraint rebuild. + 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) + # obj=-2x-y: unique optimal at x+y=8 is x=8, y=0 (x has larger coefficient) + 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) + # With y=2 fixed: x<=6, objective=-2*6-2=-14. Unique. + _solve_and_check( + self, self.opt, m, {'objective': -14.0, 'vars': [(m.x, 6.0), (m.y, 2.0)]} + ) + m.y.unfix() + # After unfix: back to unique x=8, y=0. + _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): + # After remove_constraints, the mutable helper for that row must be + # dropped from _mutable_helpers (otherwise stale handle is dereferenced + # on the next param update). + 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) + # Should not raise even though the only mutable helper was removed. + _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): + # _warmstart stores Python list-position j as the Xpress column index. + # After _remove_variables, Xpress renumbers columns. The invariant + # self._vars[j] <-> Xpress column j must hold; otherwise addMipSol + # receives warm-start values in the wrong columns. + 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) + # Remove the middle variable -- this compacts self._vars and renumbers + # Xpress columns. Afterwards self._vars = [x, z]. + self.opt.remove_variables([m.y]) + # Remove m.y from the Pyomo model so model.nvariables() reflects only the + # two remaining variables (x and z). The _solve_and_check assertion + # requires vars to cover all active Pyomo variables. + del m.y + # Verify that Python list position j equals the Xpress column index for + # every remaining variable. If they diverge, _warmstart passes wrong columns. + 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}", + ) + # Set feasible warm-start hint values on the remaining variables and re-solve. + # If column indices are wrong, addMipSol passes values to incorrect columns + # which could produce a wrong solution or a solver error. + m.x.set_value(1) + m.z.set_value(0) + # Remove the constraint that referenced m.y so the Pyomo model is consistent + # with the reduced Xpress problem (y column deleted, constraint cleaned up). + 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 connector.""" + + 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): + # NL constraint -> registered in _mutable_helpers with nl_expr set. + 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): + # NL constraint with no mutable params: still registered in _mutable_helpers + # (nl_expr set) because the NL row rebuild path is always available; + # no mutable coef dicts populated since there are no mutable params. + 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): + # sin(x) + p*y <= 5: the whole constraint is NL (full row rebuild on any + # param change); p is stored in _lin_coefs and re-evaluated at rebuild time. + 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) + # With p=1: sin(x) + y <= 5, maximize y -> y = 5 - sin(x) near 5 + 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) + # With p=2: 2y must leave room for sin(x), so y <= (5 - sin(x))/2 < y1 + self.assertLess(y2, y1) + + def test_nl_mutable_nl_coef_full_rebuild(self): + # p*sin(x) <= 0.5: mutable param inside NL -> full row rebuild. + # Bound x to [0, pi/2] so sin is strictly increasing and the constraint + # x <= arcsin(0.5/p) is always the binding upper limit. + # Objective maximize x (min -x) so the constraint is always active. + # p=1: x <= arcsin(0.5) = pi/6 ~ 0.5236. p=2: x <= arcsin(0.25) ~ 0.2527. + 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) + # Larger p tightens the constraint: x2 should be noticeably smaller than x1 + self.assertLess(x2, x1 - 0.1) + + def test_nl_mutable_bound(self): + # sin(x) >= p: mutable lower bound (NL row rebuild path when p changes). + # x in [0, pi/2] so sin is strictly increasing; constraint is always binding. + # p=0.5: x >= arcsin(0.5) = pi/6 ~ 0.5236. + # p=0.9: x >= arcsin(0.9) ~ 1.1198. + 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) + # Tighter lower bound means x must be larger + self.assertGreater(x2, x1 + 0.4) + + def test_nl_solve_modify_resolve(self): + # Solve NLP, change mutable bound param, re-solve. + # exp(x) >= p is a pure-NL constraint; p change triggers a full NL row rebuild. + # exp(x) >= p, min x: optimal is x = ln(p). + 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) + # exp(x) >= 1 -> x >= 0, minimize -> x=0 + _solve_check_mutate_check( + self, + self.opt, + m, + {'objective': 0.0, 'vars': [(m.x, 0.0)]}, + m.p, + math.e, # exp(x) >= e -> x >= 1, minimize -> x=1 + {'objective': 1.0, 'vars': [(m.x, 1.0)]}, + ) + + def test_fix_variable_no_nl_constraint_rebuild(self): + # Fixing a variable in an NL constraint must NOT rebuild the row -- + # Xpress handles fixed vars via bounds natively. Row count equality is + # a weak proxy (del+add restores count). Use xp_con handle identity as + # the definitive proof: if the handle survives, no del+add occurred. + 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) + # First solve: get internal state before fixing x. + # sin(x)+y<=5, min y -> y=0. x free, so x can be anything in [0,pi]. + _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) # fix x at pi/6 -> sin(x)=0.5, constraint: 0.5+y<=5 + # Second solve: sin(pi/6)=0.5, constraint becomes 0.5+y<=5, min y -> y=0. + _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) + # Handle identity: same xp.constraint object means no del+add occurred + self.assertIs(self.opt._mutable_helpers[m.c]._xp_con, xp_con_before) + + def _run_nl_linear_shared_param_test(self, nl_first: bool): + """Helper: same param in NL (rebuild) and linear (chgMCoef) constraints. + nl_first=True -> NL at row 0, linear at row 1 (delConstraint(NL) shifts linear). + nl_first=False -> linear at row 0, NL at row 1 (no shift on linear). + Both orderings must produce the correct solution after param change. + + The NL constraint IS binding at the optimum (p*sin(x) <= 0.5 limits x). + This is necessary so that a silently dropped NL rebuild would be detectable. + Objective max x+y so both the NL and linear constraints contribute to the result. + p=1: sin(x)<=0.5->x<=pi/6~0.524, y<=4 -> x+y~4.524. + p=2: 2*sin(x)<=0.5->x<=arcsin(0.25)~0.253, y<=2 -> x+y~2.253.""" + 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): + # NL at row 0, linear at row 1: delConstraint(NL) shifts linear's row. + # The single-pass bug would corrupt the linear constraint update here. + self._run_nl_linear_shared_param_test(nl_first=True) + + def test_param_shared_linear_before_nl(self): + # Linear at row 0, NL at row 1: no row shift on linear after NL deletion. + self._run_nl_linear_shared_param_test(nl_first=False) + + def test_param_shared_multiple_nl_and_linear(self): + # Multiple NL constraints interleaved with linear ones. + # Each NL rebuild (delConstraint) can shift remaining row indices. + 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)) + # Declare in interleaved order: NL, linear, NL, linear + 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): + # Mutable quadratic coef (has_mutable_quadratic -> _nl_rebuild_cons) at + # a lower row than a linear mutable constraint. Same row-shift risk as NL. + 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 + ) # mutable_quadratic -> NL rebuild + m.c_lin = pyo.Constraint(expr=m.p * m.y <= 4) # mutable_lin -> chgMCoef + 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): + # sin(p*x) + q*y <= 5: p is inside the NL formula (-> full rebuild on + # p change), q is a mutable linear coef stored in _mutable_helpers[con]._lin_coefs + # (re-evaluated at rebuild time via _make_xp_con). Constraint is in _mutable_helpers. + 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)]} + ) + # NL rebuild registered; mutable lin coef (q) stored in same helper. + self.assertIn(m.c, self.opt._mutable_helpers) + self.assertIsNotNone(self.opt._mutable_helpers[m.c]._nl_expr) + pyo.value(m.y) + + # Change affine param q only: constraint updated, solution changes + m.q.set_value(2.0) + # q=2, p=1: sin(1*x) + 2*y <= 5. At x~0, y <= 5/2 = 2.5. + _solve_and_check( + self, self.opt, m, {'objective': -2.5, 'vars': [(m.x, 0.0), (m.y, 2.5)]} + ) + pyo.value(m.y) + + # Change NL formula param p: full rebuild, solution changes + m.p.set_value(0.0) # sin(0)=0, so constraint becomes q*y <= 5, y <= 2.5 + _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): + # Mutable param in objective AND in a linear constraint. + # setObjective does not shift constraint row indices, so the linear + # update must be applied correctly regardless of objective rebuild. + 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) # NL rebuild + m.c_lin = pyo.Constraint(expr=m.p * m.y <= 4) # chgMCoef + m.obj = pyo.Objective(expr=m.p * m.y, sense=pyo.maximize) # mutable obj + 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) + # p=2: y<=2, obj = 2*y = 4. + _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): + # Coverage gap: xpress_persistent.py lines 287-295 (_MutableObjective NL rebuild). + # Objective = p*sin(x), mutable p inside NL formula -> full NL objective rebuild. + 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) + # p=1, max sin(x) on [0, pi/2] -> x = pi/2; sense=maximize -> obj = sin(pi/2) = 1.0 + _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) # now minimizing sin(x) -> x = 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): + # Regression test for a mutation bug in _MutableObjective.update(). + # + # Trigger conditions: + # - NL objective + # - stable (non-mutable) linear terms -> stable_xp is xp.expression + # - no mutable linear/quad terms -> mut_terms stays empty + # - mutable constant in repn -> _constant is set + # + # Model: obj = 2*x + p + sin(z) + # 2*x : stable linear term (non-mutable coef) -> goes into stable_xp + # p : mutable param, no variable -> becomes repn.constant -> _constant + # sin(z): NL part -> nl_expr, mut_terms is empty + 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)) + # min 2*x + p + sin(z): with x,z in bounds, optimum is x=0, z=0. + # obj_value = 2*0 + p + sin(0) = p + 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)]}, + ) + + # After the first update p=3: obj = 2*0 + 3 + sin(0) = 3. + # With the bug, stable_xp has been mutated to (2*x + 1), so the next + # update would produce (2*x + 1 + 3 + sin(z)) = 4 instead of 3. + 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)]}, + ) + + # A third solve with p=5 would give 7 with the bug (accumulated 1+3+5-2 offset) + # but should give 5. + 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): + # Coverage gap: xpress_persistent.py line 619 (stable quadratic in NL constraint). + # 2*x^2 is a non-mutable quadratic term; it goes into stable_xp at registration. + # p*exp(y) is the NL mutable term that triggers row rebuilds. + # Objective max x (min -x) so the constraint is always the binding limit. + # p=1: 2*x^2 + exp(0) <= 5 -> x <= sqrt(2) ~ 1.414. + # p=2: 2*x^2 + 2*exp(0) <= 5 -> x <= sqrt(1.5) ~ 1.225. + 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) + + # Change mutable param: stable quad term must survive the NL rebuild. + 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): + # Coverage gap: xpress_persistent.py lines 715/726 (stable lin/quad in NL objective). + # 2*x is stable linear, y^2 is stable quadratic, p*sin(z) is NL with mutable p. + 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)) + # objective = 2*x + y^2 + p*sin(z): stable lin (2x), stable quad (y^2), NL (p*sin) + m.obj = pyo.Objective(expr=2 * m.x + m.y**2 + m.p * pyo.sin(m.z)) + # Minimizing: x=0 (coef 2>0), y=0 (convex quad), z=0 (p=1>0, sin>=0) -> obj=0 + _solve_and_check( + self, + self.opt, + m, + {'objective': 0.0, 'vars': [(m.x, 0.0), (m.y, 0.0), (m.z, 0.0)]}, + ) + + # Change mutable param (triggers NL objective rebuild) + m.p.set_value(-1.0) # now -sin(z) is negative -> want max z + # Minimizing 2x + y^2 - sin(z): x=0, y=0, z=pi/2 -> obj = -1 + _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): + # p*x^3 >= 1, min x, x in [0,10]. + # p=1: x^3 >= 1 -> x >= 1 -> optimal x=1.0 + # p=8: 8*x^3 >= 1 -> x^3 >= 1/8 -> x >= 0.5 -> optimal x=0.5 + 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): + # Start with linear at row 0, NL at row 1 (safe ordering -- no shift). + # Remove linear: NL shifts to row 0. + # Re-add linear: NL=0, linear=1 (now the dangerous ordering). + # A param change must still correctly update the linear constraint. + 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) # row 0 initially + m.c_nl = pyo.Constraint(expr=m.p * pyo.sin(m.x) <= 5) # row 1 initially + 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)]}, + ) + + # Remove c_lin (row 0): c_nl shifts to row 0. + # Correct Pyomo API ordering: deactivate first, then remove from solver. + m.c_lin.deactivate() + self.opt.remove_constraints([m.c_lin]) + # Re-add c_lin: now c_nl=row 0, c_lin=row 1 (dangerous ordering). + m.c_lin.activate() + self.opt.add_constraints([m.c_lin]) + # c_lin must be re-registered in _mutable_helpers after re-add. + self.assertIn(m.c_lin, self.opt._mutable_helpers) + + m.p.set_value(2.0) + # p=2: c_lin becomes p*y<=4 -> y<=2 + _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): + """Tests for IIS computation via write_iis() and get_iis().""" + + def setUp(self): + self.opt = XpressPersistent() + + def _infeasible_model(self): + # x in {0,1}, y >= 0 + # c1: y <= 100*x c2: y <= -100*x c3: x >= 0.5 + # c2 and c3 together with x binary force infeasibility. + 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() + # The two infeasible constraints must appear in the IIS file. + 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']} + # c2 and c3 are the infeasible pair; c1 may or may not be included. + self.assertIn('c2', con_names) + self.assertIn('c3', con_names) + # y's lower bound (0) contributes to the IIS. + var_names = {v.name for v in iis['variables']} + self.assertIn('y', var_names) + + def test_get_iis_objects_are_model_constraints(self): + # Verify the returned ConstraintData and VarData objects are the + # actual Pyomo components, not copies or proxies (identity check). + 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..3fdb9f02c08 --- /dev/null +++ b/pyomo/contrib/solver/tests/solvers/test_xpress_walker.py @@ -0,0 +1,743 @@ +# ____________________________________________________________________________________ +# +# 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: verify that XpressExpressionWalker produces correct xp expressions. + +Mutable parameter tracking is now done by generate_standard_repn in +xpress_persistent.py, not by the walker. These tests verify only xp_expr +building correctness across LP/MIP/QP/NLP expression patterns. +""" + +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. + + expected_type: xp.expression (linear/mixed), xp.linterm (single scaled var), + xp.quadterm (pure quad), xp.nonlin (nonlinear). + expected_lin: [] = assert no linear terms; [(pyo_var, coef), ...] = exact. + expected_quad: [] = assert no quadratic terms; [(v1, v2, coef), ...] = exact. + Quadratic lookup is symmetric (x*y == y*x). + + NLP formula checks: + expected_nlp_vars: list of Pyomo VarData -- all must appear as COL tokens. + expected_nlp_scalars: list of float -- all must appear as CON tokens. + expected_nlp_operators: list of str -- each must appear at least once as an + IFUN or OP token. Supported names (case-insensitive): + 'sin','cos','tan','asin','acos','atan', + 'exp','ln','log10','sqrt','abs','min','max', + 'uminus','mul','div','plus','minus','pow'. + + Returns the xp expression produced by the walker. + """ + result = walker.walk_expression(expr) + + tc = test_case + var_map = walker.var_map + + tc.assertIsInstance(result, expected_type) + + # Test linear part + 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) + + # Test quadratic part + 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 constraint is a valid Xpress constraint + 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)) + + # NLP formula checks -- only when at least one expected_nlp_* is provided. + tok_types, tok_values = walker.prob.nlpGetFormula(test_con, 0) + + # Token type constants (from xprs.h) + _TOK_CON = 1 # numeric constant + _TOK_COL = 10 # variable column index + _TOK_IFUN = 12 # intrinsic function + _TOK_OP = 31 # arithmetic operator + # Intrinsic function codes + _IFUN = { + 'log10': 14, + 'ln': 15, + 'exp': 16, + 'abs': 17, + 'sqrt': 18, + 'sin': 27, + 'cos': 28, + 'tan': 291, + 'asin': 30, + 'arcsin': 30, + 'acos': 31, + 'arccos': 31, + 'atan': 32, + 'arctan': 32, + 'min': 33, + 'max': 34, + } + # Arithmetic operator codes + _OP = {'uminus': 1, 'pow': 2, 'mul': 3, 'div': 4, 'plus': 5, 'minus': 6} + + col_indices = {int(v) for t, v in zip(tok_types, tok_values) if int(t) == _TOK_COL} + scalars = [v for t, v in zip(tok_types, tok_values) if int(t) == _TOK_CON] + ifun_codes = {int(v) for t, v in zip(tok_types, tok_values) if int(t) == _TOK_IFUN} + op_codes = {int(v) for t, v in zip(tok_types, tok_values) if int(t) == _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, xv, yv, _ = _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, xv, yv, _ = _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): + # 0*x is folded to constant 0 by _before_monomial; only y remains as a term. + m, w, xv, yv, _ = _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): + # x fixed at 2.0 must still appear as an xp.var column, not folded to 6.0. + m, w, xv, yv, _ = _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): + # The scalar 5.0 is embedded as a constant; extractLinear does not return it. + m, w, xv, yv, _ = _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): + # 0*y -> constant 0 (not a linear term); sin(x) is the only meaningful part. + m, w, xv, yv, _ = _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): + # x + p where p=3.0 (mutable scalar). The constant is embedded; only x linear. + m, w, xv, yv, _ = _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, xv, yv, _ = _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, xv, yv, _ = _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): + # _before_linear branch (a): bare VarData with implicit coefficient 1. + m, w, xv, yv, _ = _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, xv, yv, _ = _make_walker() + result = w.walk_expression(3.0 + 2.0) + self.assertAlmostEqual(result, 5.0) + + def test_linear_fast_path(self): + """_before_linear handles LinearExpression directly without full walker descent.""" + + m, w, xv, yv, _ = _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): + """_before_linear must return (False, result) -- False prevents descent into children.""" + + m, w, xv, yv, _ = _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, xv, yv, _ = _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, xv, yv, _ = _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, xv, yv, _ = _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, xv, yv, _ = _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, xv, yv, _ = _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, xv, yv, _ = _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, xv, yv, _ = _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, xv, yv, _ = _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, xv, yv, _ = _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): + # Xpress stores x*y*z as a pure NL term; col list is empty in the formula. + m, w, xv, yv, zv = _make_walker() + _walk_and_assert(self, w, m.x * m.y * m.z, expected_type=xp.nonlin) + + @parameterized.expand( + [ + # (name, fn, expected_nlp_scalars, primary_operator) + # sin/cos map directly; hyperbolic decompose via exp so scalars appear. + # sinh=0.5*(exp(x)-exp(-x)): scalar 0.5 + # cosh=0.5*(exp(x)+exp(-x)): scalar 0.5 + # tanh=(exp(x)-exp(-x))/(exp(x)+exp(-x)): no scalars + # asinh=ln(x+sqrt(x^2+1)): scalar 1.0 + # acosh=ln(x+sqrt(x^2-1)): scalar -1.0 + # atanh=0.5*ln((1+x)/(1-x)): scalars [1.0, 1.0, 0.5] + ('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): + # fn(x) + y: NL formula contains x and the primary operator. + # y is in the linear part (extractLinear), not in the NLP formula. + m, w, xv, yv, _ = _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): + # p*sin(x): p=2.0 is folded as a CON token in the NLP formula. + m, w, xv, yv, _ = _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): + # sin(x) + p*y: sin(x) is in the NLP formula; p*y is linear. + m, w, xv, yv, _ = _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): + # sin(p*x): p=2.0 is a CON token inside the argument. + m, w, xv, yv, _ = _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, xv, yv, zv = _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, xv, yv, _ = _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, xv, yv, zv = _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, xv, _, _ = _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, xv, yv, _ = _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): + # sin(x)*sin(y): product of two NL expressions. Both variables and both + # sin operators appear in the row formula. + m, w, xv, yv, _ = _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): + # sin(p*x)*sin(y): mutable param inside the NL argument of a product of two + # NL expressions. Exercises the ProductExpression handler with a mutable param. + m, w, xv, yv, _ = _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): + """NPV sub-expression m.q/m.r with r=0 raises ZeroDivisionError at walk time. + _before_npv calls value() on the sub-expression; the error propagates out.""" + m, w, xv, yv, _ = _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): + """sqrt(-1) in an NPV sub-expression raises ValueError from value() -- loud, + not silent nan. _before_npv has no try-except so domain errors propagate.""" + + m, w, xv, _, _ = _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): + """NPV sub-expression log(m.p) with p<0 raises ValueError at walk time. + _before_npv calls value() on the sub-expression; the error propagates out.""" + m, w, xv, _, _ = _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): + """_exit_named_expression must insert into subexpression_cache on first walk.""" + m, w, xv, yv, _ = _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): + """Named expression walked once; cache grows only on first walk.""" + m, w, xv, yv, _ = _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) # Hit + + def test_named_expr_produces_valid_xp_expr(self): + """Named expression with mutable coef: valid xp expression on cache miss.""" + m, w, xv, yv, _ = _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): + """ceil is not supported by Xpress -- must raise IncompatibleModelError.""" + m, w, xv, yv, _ = _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): + """floor is not supported by Xpress -- must raise IncompatibleModelError.""" + m, w, xv, yv, _ = _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, xv, yv, _ = _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): + """_ExitHandlerMap.__missing__ raises IncompatibleModelError when the + expression type's MRO has no registered exit handler.""" + + class _CustomUnregisteredExpr: + """Synthetic class that is not part of the Pyomo expression hierarchy.""" + + with self.assertRaises(IncompatibleModelError): + _EXIT_HANDLERS[_CustomUnregisteredExpr] + + +if __name__ == '__main__': + unittest.main() From 38f73d33736b16cad911d1378a21b0d38b32c639 Mon Sep 17 00:00:00 2001 From: Francesco Cavaliere Date: Tue, 30 Jun 2026 14:36:08 +0200 Subject: [PATCH 2/2] Formatting last few additions that were left off --- .../solver/solvers/xpress/xpress_persistent.py | 4 +++- .../solver/tests/solvers/_xpress_test_utils.py | 14 ++++++++------ .../solver/tests/solvers/test_xpress_persistent.py | 10 ++++------ 3 files changed, 15 insertions(+), 13 deletions(-) diff --git a/pyomo/contrib/solver/solvers/xpress/xpress_persistent.py b/pyomo/contrib/solver/solvers/xpress/xpress_persistent.py index 34cefb42ab4..a5a54f3003f 100644 --- a/pyomo/contrib/solver/solvers/xpress/xpress_persistent.py +++ b/pyomo/contrib/solver/solvers/xpress/xpress_persistent.py @@ -711,7 +711,9 @@ def _remove_constraints(self, cons: list[ConstraintData]) -> 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) + 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: diff --git a/pyomo/contrib/solver/tests/solvers/_xpress_test_utils.py b/pyomo/contrib/solver/tests/solvers/_xpress_test_utils.py index f52ec46940c..579d2e9a603 100644 --- a/pyomo/contrib/solver/tests/solvers/_xpress_test_utils.py +++ b/pyomo/contrib/solver/tests/solvers/_xpress_test_utils.py @@ -41,11 +41,11 @@ def _simple_mip(): 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 + 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): @@ -135,7 +135,9 @@ def _solve_check_mutate_check( Only use when there are no assertions that compare values across the two solves (e.g. assertLess(x2, x1)). Those tests must keep explicit solve calls. """ - res_before = _solve_and_check(test_case, opt, model, expected_before, **solve_kwargs) + 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_xpress_persistent.py b/pyomo/contrib/solver/tests/solvers/test_xpress_persistent.py index 3e11cdde804..d18616ade27 100644 --- a/pyomo/contrib/solver/tests/solvers/test_xpress_persistent.py +++ b/pyomo/contrib/solver/tests/solvers/test_xpress_persistent.py @@ -1305,6 +1305,7 @@ def test_nl_objective_stable_xp_not_mutated_by_constant_update(self): # p : mutable param, no variable -> becomes repn.constant -> _constant # sin(z): NL part -> nl_expr, mut_terms is empty import math + m = pyo.ConcreteModel() m.p = pyo.Param(mutable=True, initialize=1.0) m.x = pyo.Var(bounds=(0, 1)) @@ -1313,8 +1314,7 @@ def test_nl_objective_stable_xp_not_mutated_by_constant_update(self): # obj_value = 2*0 + p + sin(0) = p 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)]}, + self, self.opt, m, {'objective': 1.0, 'vars': [(m.x, 0.0), (m.z, 0.0)]} ) # After the first update p=3: obj = 2*0 + 3 + sin(0) = 3. @@ -1322,16 +1322,14 @@ def test_nl_objective_stable_xp_not_mutated_by_constant_update(self): # update would produce (2*x + 1 + 3 + sin(z)) = 4 instead of 3. 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)]}, + self, self.opt, m, {'objective': 3.0, 'vars': [(m.x, 0.0), (m.z, 0.0)]} ) # A third solve with p=5 would give 7 with the bug (accumulated 1+3+5-2 offset) # but should give 5. 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)]}, + self, self.opt, m, {'objective': 5.0, 'vars': [(m.x, 0.0), (m.z, 0.0)]} ) def test_nl_constraint_stable_quadratic_term(self):