diff --git a/pyomo/contrib/iis/tests/test_iis.py b/pyomo/contrib/iis/tests/test_iis.py index ef892241ece..42553ac396f 100644 --- a/pyomo/contrib/iis/tests/test_iis.py +++ b/pyomo/contrib/iis/tests/test_iis.py @@ -118,17 +118,27 @@ def _test_iis(solver_name): def _validate_ilp(file_name): - lines_found = {"c2: 100 x + y <= 0": False, "c3: x >= 0.5": False} + # Xpress may write constraint terms in different orders across versions + # (e.g., "100 x + y <= 0" vs "y + 100 x <= 0"), so check both forms for c2. + c2_forms = {"c2: 100 x + y <= 0", "c2: y + 100 x <= 0"} + c2_found = False + c3_found = False with open(file_name, "r") as f: for line in f.readlines(): - for k, v in lines_found.items(): - if (not v) and k in line: - lines_found[k] = True - - if not all(lines_found.values()): + if any(form in line for form in c2_forms): + c2_found = True + if "c3: x >= 0.5" in line: + c3_found = True + + missing = [] + if not c2_found: + missing.append("c2: 100 x + y <= 0") + if not c3_found: + missing.append("c3: x >= 0.5") + if missing: raise Exception( f"The file {file_name} is not as expected. Missing constraints:\n" - + "\n".join(k for k, v in lines_found.items() if not v) + + "\n".join(missing) ) diff --git a/pyomo/repn/plugins/standard_form.py b/pyomo/repn/plugins/standard_form.py index 70b1b74b8cc..e565133b354 100644 --- a/pyomo/repn/plugins/standard_form.py +++ b/pyomo/repn/plugins/standard_form.py @@ -82,15 +82,31 @@ class LinearStandardFormInfo: rhs : numpy.ndarray - The constraint right-hand sides. + The constraint right-hand sides. For range rows (``bound_type + == 2``, only produced when ``keep_range_constraints=True``), + this holds the adjusted *upper* bound ``ub - offset``. + + rhs_range : numpy.ndarray or None + + Range widths for range rows: ``rhs_range[i] = ub - lb`` for + ``bound_type == 2`` rows. ``None`` when + ``keep_range_constraints=False`` (the default) or when the + model contains no range constraints. When not ``None``, the + adjusted lower bound for a range row can be recovered as + ``rhs[i] - rhs_range[i]``. rows : List[Tuple[ConstraintData, int]] The list of Pyomo constraint objects corresponding to the rows in `A`. Each element in the list is a 2-tuple of - (ConstraintData, row_multiplier). The `row_multiplier` will be - +/- 1 indicating if the row was multiplied by -1 (corresponding - to a constraint lower bound) or +1 (upper bound). + (ConstraintData, bound_type). ``bound_type`` values: + + * ``+1`` – upper-bound row (``Ax ≤ rhs``); + * ``-1`` – lower-bound row (see mode-dependent sign conventions); + * ``0`` – equality row (``mixed_form`` only); + * ``+2`` – range row (``lb - offset ≤ Ax ≤ ub - offset``, + coefficients in the upper-bound sense; only produced when + ``keep_range_constraints=True``). columns : List[VarData] @@ -99,7 +115,9 @@ class LinearStandardFormInfo: objectives : List[ObjectiveData] - The list of Pyomo objective objects corresponding to the active objectives + The list of Pyomo objective objects corresponding to the active + objectives whose expressions are purely linear (and thus appear + in `c`). eliminated_vars: List[Tuple[VarData, NumericExpression]] @@ -111,17 +129,47 @@ class LinearStandardFormInfo: all variables appearing in the expression must either have appeared in the standard form, or appear *earlier* in this list. + nonlinear_constraints : List[ConstraintData] or None + + Constraints skipped because they contain nonlinear terms. ``None`` + when ``allow_nonlinear=False`` (the default). When + ``allow_nonlinear=True``, holds the list of constraints with nonlinear + terms that were omitted from the compiled matrices (may be empty). + + nonlinear_objectives : List[ObjectiveData] or None + + Objectives skipped because they contain nonlinear terms. ``None`` + when ``allow_nonlinear=False`` (the default). When + ``allow_nonlinear=True``, holds the list of objectives with nonlinear + terms that were omitted from the compiled matrices (may be empty). + """ - def __init__(self, c, c_offset, A, rhs, rows, columns, objectives, eliminated_vars): + def __init__( + self, + c, + c_offset, + A, + rhs, + rhs_range, + rows, + columns, + objectives, + eliminated_vars, + nonlinear_constraints=None, + nonlinear_objectives=None, + ): self.c = c self.c_offset = c_offset self.A = A self.rhs = rhs + self.rhs_range = rhs_range self.rows = rows self.columns = columns self.objectives = objectives self.eliminated_vars = eliminated_vars + self.nonlinear_constraints = nonlinear_constraints + self.nonlinear_objectives = nonlinear_objectives @property def x(self): @@ -180,6 +228,19 @@ class LinearStandardFormCompiler: 'mix of <=, ==, and >=)', ), ) + CONFIG.declare( + 'keep_range_constraints', + ConfigValue( + default=False, + domain=bool, + description='Emit range constraints (finite lb ≠ ub) as a single ' + 'row with bound_type=2 rather than splitting them into separate ' + 'upper- and lower-bound rows. The rhs entry for such a row is the ' + 'adjusted upper bound (ub - offset); the range width (ub - lb) is ' + 'stored in the rhs_range array of the returned ' + 'LinearStandardFormInfo. Cannot be combined with slack_form.', + ), + ) CONFIG.declare( 'set_sense', ConfigValue( @@ -188,6 +249,30 @@ class LinearStandardFormCompiler: description='If not None, map all objectives to the specified sense.', ), ) + CONFIG.declare( + 'allow_nonlinear', + ConfigValue( + default=False, + domain=bool, + description='If True, constraints and objectives containing nonlinear ' + 'terms are collected into ``LinearStandardFormInfo.nonlinear_constraints`` ' + 'and ``LinearStandardFormInfo.nonlinear_objectives`` rather than raising ' + 'an exception. The nonlinear components are omitted from the compiled ' + 'matrices.', + ), + ) + CONFIG.declare( + 'extra_valid_ctypes', + ConfigValue( + default=[], + description='Additional component types that are permitted to appear ' + 'in the model without causing an error, but that are not processed by ' + 'the compiler. Use this when the model contains component types ' + '(e.g., :class:`~pyomo.core.base.sos.SOSConstraint`) that are valid ' + 'for the calling solver but that the standard-form compiler does not ' + 'know how to handle.', + ), + ) CONFIG.declare( 'show_section_timing', ConfigValue( @@ -309,7 +394,8 @@ def write(self, model): RangeSet, Port, # TODO: Piecewise, Complementarity - }, + } + | set(self.config.extra_valid_ctypes), targets={Suffix, Objective}, ) if unknown: @@ -364,6 +450,7 @@ def write(self, model): # Process objective # set_sense = self.config.set_sense + allow_nonlinear = self.config.allow_nonlinear objectives = [] for blk in component_map[Objective]: objectives.extend( @@ -376,29 +463,54 @@ def write(self, model): obj_data = [] obj_index = [] obj_index_ptr = [0] + linear_objectives = [] + nonlinear_objectives = [] for obj in objectives: if hasattr(obj, 'template_expr'): - offset, linear_index, linear_data, lb, ub = ( - template_visitor.expand_expression(obj, obj.template_expr()) - ) + if allow_nonlinear: + try: + offset, linear_index, linear_data, lb, ub = ( + template_visitor.expand_expression(obj, obj.template_expr()) + ) + except InvalidExpressionError: + nonlinear_objectives.append(obj) + if with_debug_timing: + timer.toc( + 'Objective %s (nonlinear)', obj, level=logging.DEBUG + ) + continue + else: + offset, linear_index, linear_data, lb, ub = ( + template_visitor.expand_expression(obj, obj.template_expr()) + ) assert lb is None and ub is None N = len(linear_index) obj_index.append(linear_index) obj_data.append(linear_data) obj_offset.append(offset) + linear_objectives.append(obj) else: repn = visitor.walk_expression(obj.expr) - N = len(repn.linear) - obj_index.append(map(var_recorder.var_order.__getitem__, repn.linear)) - obj_data.append(repn.linear.values()) - obj_offset.append(repn.constant) if repn.nonlinear is not None: + if allow_nonlinear: + nonlinear_objectives.append(obj) + if with_debug_timing: + timer.toc( + 'Objective %s (nonlinear)', obj, level=logging.DEBUG + ) + continue raise InvalidExpressionError( f"Model objective ({obj.name}) contains nonlinear terms that " "cannot be compiled to standard (linear) form." ) + N = len(repn.linear) + obj_index.append(map(var_recorder.var_order.__getitem__, repn.linear)) + obj_data.append(repn.linear.values()) + obj_offset.append(repn.constant) + linear_objectives.append(obj) + obj_nnz += N if set_sense is not None and set_sense != obj.sense: obj_data[-1] = -self._to_vector(obj_data[-1], float, N) @@ -406,20 +518,28 @@ def write(self, model): obj_index_ptr.append(obj_index_ptr[-1] + N) if with_debug_timing: timer.toc('Objective %s', obj, level=logging.DEBUG) + objectives = linear_objectives # # Tabulate constraints # slack_form = self.config.slack_form mixed_form = self.config.mixed_form + keep_range_constraints = self.config.keep_range_constraints if slack_form and mixed_form: raise ValueError("cannot specify both slack_form and mixed_form") + if slack_form and keep_range_constraints: + raise ValueError( + "cannot specify both slack_form and keep_range_constraints" + ) rows = [] rhs = [] + rhs_range = [] if keep_range_constraints else None con_nnz = 0 con_data = [] con_index = [] con_index_ptr = [0] + nonlinear_constraints = [] last_parent = None for con in ordered_active_constraints(model, self.config): if with_debug_timing and con._component is not last_parent: @@ -428,9 +548,22 @@ def write(self, model): last_parent = con._component if hasattr(con, 'template_expr'): - offset, linear_index, linear_data, lb, ub = ( - template_visitor.expand_expression(con, con.template_expr()) - ) + if allow_nonlinear: + try: + offset, linear_index, linear_data, lb, ub = ( + template_visitor.expand_expression(con, con.template_expr()) + ) + except InvalidExpressionError: + nonlinear_constraints.append(con) + if with_debug_timing: + timer.toc( + 'Constraint %s (nonlinear)', con, level=logging.DEBUG + ) + continue + else: + offset, linear_index, linear_data, lb, ub = ( + template_visitor.expand_expression(con, con.template_expr()) + ) N = len(linear_data) else: # Note: lb and ub could be a number, expression, or None. @@ -442,6 +575,13 @@ def write(self, model): ub = value(ub) repn = visitor.walk_expression(body) if repn.nonlinear is not None: + if allow_nonlinear: + nonlinear_constraints.append(con) + if with_debug_timing: + timer.toc( + 'Constraint %s (nonlinear)', con, level=logging.DEBUG + ) + continue raise InvalidConstraintError( f"Model constraint ({con.name}) contains nonlinear terms that " "cannot be compiled to standard (linear) form." @@ -453,6 +593,14 @@ def write(self, model): linear_index = map(var_recorder.var_order.__getitem__, repn.linear) linear_data = repn.linear.values() + # Normalize ±inf to None: both kernel constraints and AML + # RangedExpressions can return ±inf instead of None for unbounded + # sides (e.g., `(-inf, x, 5)`). Treat them as unbounded. + if lb == -float('inf'): + lb = None + if ub == float('inf'): + ub = None + if lb is None and ub is None: # Note: you *cannot* output trivial (unbounded) # constraints in matrix format. I suppose we could add a @@ -473,6 +621,18 @@ def write(self, model): con_nnz += N rows.append(RowEntry(con, 0)) rhs.append(ub - offset) + if keep_range_constraints: + rhs_range.append(0.0) + con_data.append(linear_data) + con_index.append(linear_index) + con_index_ptr.append(con_nnz) + elif lb is not None and ub is not None and keep_range_constraints: + # Range constraint: single row, coefficients in the upper- + # bound sense (not negated), bound_type=2. + con_nnz += N + rows.append(RowEntry(con, 2)) + rhs.append(ub - offset) + rhs_range.append(ub - lb) con_data.append(linear_data) con_index.append(linear_index) con_index_ptr.append(con_nnz) @@ -483,6 +643,8 @@ def write(self, model): con_nnz += N rows.append(RowEntry(con, 1)) rhs.append(ub - offset) + if keep_range_constraints: + rhs_range.append(0.0) con_data.append(linear_data) con_index.append(linear_index) con_index_ptr.append(con_nnz) @@ -490,6 +652,8 @@ def write(self, model): con_nnz += N rows.append(RowEntry(con, -1)) rhs.append(lb - offset) + if keep_range_constraints: + rhs_range.append(0.0) con_data.append(linear_data) con_index.append(linear_index) con_index_ptr.append(con_nnz) @@ -524,22 +688,37 @@ def write(self, model): con_index.append(linear_index) con_index_ptr.append(con_nnz) else: - if ub is not None: - if lb is not None: - linear_index = list(linear_index) + is_range = lb is not None and ub is not None and lb != ub + if is_range and keep_range_constraints: + # Range constraint: single row, bound_type=2. con_nnz += N - rows.append(RowEntry(con, 1)) + rows.append(RowEntry(con, 2)) rhs.append(ub - offset) + rhs_range.append(ub - lb) con_data.append(linear_data) con_index.append(linear_index) con_index_ptr.append(con_nnz) - if lb is not None: - con_nnz += N - rows.append(RowEntry(con, -1)) - rhs.append(offset - lb) - con_data.append(-np.array(list(linear_data))) - con_index.append(linear_index) - con_index_ptr.append(con_nnz) + else: + if ub is not None: + if lb is not None: + linear_index = list(linear_index) + con_nnz += N + rows.append(RowEntry(con, 1)) + rhs.append(ub - offset) + if keep_range_constraints: + rhs_range.append(0.0) + con_data.append(linear_data) + con_index.append(linear_index) + con_index_ptr.append(con_nnz) + if lb is not None: + con_nnz += N + rows.append(RowEntry(con, -1)) + rhs.append(offset - lb) + if keep_range_constraints: + rhs_range.append(0.0) + con_data.append(-np.array(list(linear_data))) + con_index.append(linear_index) + con_index_ptr.append(con_nnz) if with_debug_timing: # report the last constraint @@ -613,7 +792,17 @@ def write(self, model): eliminated_vars = [] info = LinearStandardFormInfo( - c, obj_offset, A, rhs, rows, columns, objectives, eliminated_vars + c, + obj_offset, + A, + rhs, + np.array(rhs_range) if keep_range_constraints else None, + rows, + columns, + objectives, + eliminated_vars, + nonlinear_constraints=nonlinear_constraints if allow_nonlinear else None, + nonlinear_objectives=nonlinear_objectives if allow_nonlinear else None, ) timer.toc("Generated linear standard form representation", delta=False) return info diff --git a/pyomo/repn/tests/test_standard_form.py b/pyomo/repn/tests/test_standard_form.py index eaa515c21bc..90de0b6075c 100644 --- a/pyomo/repn/tests/test_standard_form.py +++ b/pyomo/repn/tests/test_standard_form.py @@ -351,6 +351,222 @@ def test_alternative_forms(self): self.assertTrue(np.all(repn.c == ref)) self._verify_solution(soln, repn, True) + def test_keep_range_constraints(self): + m = pyo.ConcreteModel() + m.x = pyo.Var() + m.y = pyo.Var([0, 1, 3], bounds=lambda m, i: (0, 10)) + # Pure lower-bound constraint + m.c = pyo.Constraint(expr=m.x + 2 * m.y[1] >= 3) + # Pure upper-bound constraint + m.d = pyo.Constraint(expr=m.y[1] + 4 * m.y[3] <= 5) + # Range constraint: -2 <= y[0] + 1 + 6*y[1] <= 7 → -3 <= y[0] + 6*y[1] <= 6 + m.e = pyo.Constraint(expr=pyo.inequality(-2, m.y[0] + 1 + 6 * m.y[1], 7)) + # Equality + m.f = pyo.Constraint(expr=m.x + m.y[0] == 8) + m.o = pyo.Objective(expr=5 * m.x) + + col_order = [m.x, m.y[0], m.y[1], m.y[3]] + + # --- mixed_form + keep_range_constraints --- + repn = LinearStandardFormCompiler().write( + m, mixed_form=True, keep_range_constraints=True, column_order=col_order + ) + # m.e: single range row (bound_type=2); all others are normal mixed rows + self.assertEqual(repn.rows, [(m.c, -1), (m.d, 1), (m.e, 2), (m.f, 0)]) + ref_A = np.array([[1, 0, 2, 0], [0, 0, 1, 4], [0, 1, 6, 0], [1, 1, 0, 0]]) + self.assertTrue(np.all(repn.A.toarray() == ref_A)) + # m.e: rhs = ub - offset = 7 - 1 = 6 + self.assertTrue(np.all(repn.rhs == np.array([3, 5, 6, 8]))) + # rhs_range: only m.e is a range row; range = 7 - (-2) = 9 + self.assertTrue(np.all(repn.rhs_range == np.array([0.0, 0.0, 9.0, 0.0]))) + + # --- default form + keep_range_constraints --- + repn2 = LinearStandardFormCompiler().write( + m, keep_range_constraints=True, column_order=col_order + ) + # lb-only (m.c) → negated ≤ row; ub-only (m.d) → ≤ row; + # range (m.e) → single row; equality (m.f) → two rows (ub + negated lb) + self.assertEqual( + repn2.rows, [(m.c, -1), (m.d, 1), (m.e, 2), (m.f, 1), (m.f, -1)] + ) + self.assertTrue(np.all(repn2.rhs_range == np.array([0.0, 0.0, 9.0, 0.0, 0.0]))) + + # --- without keep_range_constraints m.e still splits into two rows --- + repn3 = LinearStandardFormCompiler().write( + m, mixed_form=True, column_order=col_order + ) + e_rows = [ + (r.constraint, r.bound_type) for r in repn3.rows if r.constraint is m.e + ] + self.assertEqual(e_rows, [(m.e, 1), (m.e, -1)]) + # rhs_range is None when keep_range_constraints=False + self.assertIsNone(repn3.rhs_range) + + # --- slack_form + keep_range_constraints must raise --- + with self.assertRaises(ValueError): + LinearStandardFormCompiler().write( + m, slack_form=True, keep_range_constraints=True + ) + + def test_nonlinear_fields_none_when_not_allowed(self): + m = pyo.ConcreteModel() + m.x = pyo.Var() + m.c = pyo.Constraint(expr=m.x <= 1) + m.o = pyo.Objective(expr=m.x) + + repn = LinearStandardFormCompiler().write(m, mixed_form=True) + self.assertIsNone(repn.nonlinear_constraints) + self.assertIsNone(repn.nonlinear_objectives) + + def test_allow_nonlinear_constraints(self): + m = pyo.ConcreteModel() + m.x = pyo.Var() + m.y = pyo.Var() + m.c_lin = pyo.Constraint(expr=m.x + m.y <= 5) + m.c_nl = pyo.Constraint(expr=m.x**2 + m.y <= 3) + m.o = pyo.Objective(expr=m.x + m.y) + + # Default (allow_nonlinear=False) must raise on the nonlinear constraint. + with self.assertRaises(Exception): + LinearStandardFormCompiler().write(m, mixed_form=True) + + # allow_nonlinear=True: nonlinear constraint is collected separately; + # the linear constraint still appears in A. + repn = LinearStandardFormCompiler().write( + m, mixed_form=True, allow_nonlinear=True + ) + self.assertEqual(repn.nonlinear_constraints, [m.c_nl]) + self.assertEqual(repn.nonlinear_objectives, []) + # Only the linear constraint appears in A. + self.assertEqual(len(repn.rows), 1) + self.assertEqual(repn.rows[0].constraint, m.c_lin) + # Linear objective is still compiled into c. + self.assertEqual(repn.objectives, [m.o]) + self.assertTrue(np.all(repn.c.toarray() != 0)) + + def test_allow_nonlinear_objective(self): + m = pyo.ConcreteModel() + m.x = pyo.Var() + m.y = pyo.Var() + m.c_lin = pyo.Constraint(expr=m.x + m.y <= 5) + m.o_nl = pyo.Objective(expr=m.x**2 + m.y) + + # Default must raise on the nonlinear objective. + with self.assertRaises(Exception): + LinearStandardFormCompiler().write(m, mixed_form=True) + + repn = LinearStandardFormCompiler().write( + m, mixed_form=True, allow_nonlinear=True + ) + # Nonlinear objective is NOT compiled into c; it appears in nonlinear_objectives. + self.assertEqual(repn.nonlinear_objectives, [m.o_nl]) + self.assertEqual(repn.objectives, []) + # c is empty (no linear objectives). + self.assertEqual(repn.c.shape[0], 0) + # The linear constraint is still compiled normally. + self.assertEqual(len(repn.rows), 1) + self.assertEqual(repn.rows[0].constraint, m.c_lin) + + def test_allow_nonlinear_mixed(self): + """Linear constraints/objectives compiled; nonlinear ones passed through.""" + m = pyo.ConcreteModel() + m.x = pyo.Var() + m.y = pyo.Var() + m.c_lin = pyo.Constraint(expr=m.x + 2 * m.y >= 1) + m.c_nl = pyo.Constraint(expr=m.x * m.y <= 4) + m.o_lin = pyo.Objective(expr=m.x + m.y) + + repn = LinearStandardFormCompiler().write( + m, mixed_form=True, allow_nonlinear=True + ) + + # Exactly one linear row, one nonlinear constraint. + self.assertEqual(len(repn.rows), 1) + self.assertEqual(repn.rows[0].constraint, m.c_lin) + self.assertEqual(repn.nonlinear_constraints, [m.c_nl]) + # Linear objective compiles normally. + self.assertEqual(repn.objectives, [m.o_lin]) + self.assertEqual(repn.nonlinear_objectives, []) + # Both variables appear as columns (referenced by the linear constraint). + col_ids = {id(v) for v in repn.columns} + self.assertIn(id(m.x), col_ids) + self.assertIn(id(m.y), col_ids) + + def test_inf_bounds_normalized(self): + """Constraints returning ±inf bounds are treated as unbounded (None). + + Both pyomo.kernel and AML constraints can return ±inf from + to_bounded_expression() when the user explicitly passes float('inf'). + LinearStandardFormCompiler must normalize these so that a one-sided + constraint is not misclassified as a range constraint, and a fully + unbounded constraint is skipped rather than emitted as a range row. + """ + import pyomo.kernel as pmo + + # --- kernel --- + mk = pmo.block() + mk.x = pmo.variable() + # lb=-inf (unbounded below) → should become a pure ≤ row, not a range row + mk.c_ub = pmo.constraint(ub=2.0, body=mk.x) + # ub=+inf (unbounded above) → should become a pure ≥ row, not a range row + mk.c_lb = pmo.constraint(lb=-3.0, body=mk.x) + # Explicit finite range → should still be a range row + mk.c_rng = pmo.constraint((-1.0, mk.x, 4.0)) + + repn = LinearStandardFormCompiler().write( + mk, mixed_form=True, keep_range_constraints=True + ) + by_con = {r.constraint: r.bound_type for r in repn.rows} + self.assertEqual(by_con[mk.c_ub], 1) # ≤, not range + self.assertEqual(by_con[mk.c_lb], -1) # ≥, not range + self.assertEqual(by_con[mk.c_rng], 2) # finite range + rhs_map = {r.constraint: repn.rhs[i] for i, r in enumerate(repn.rows)} + self.assertEqual(rhs_map[mk.c_ub], 2.0) + self.assertEqual(rhs_map[mk.c_lb], -3.0) + self.assertEqual(rhs_map[mk.c_rng], 4.0) # ub of range row + rr_map = {r.constraint: repn.rhs_range[i] for i, r in enumerate(repn.rows)} + self.assertEqual(rr_map[mk.c_ub], 0.0) + self.assertEqual(rr_map[mk.c_lb], 0.0) + self.assertEqual(rr_map[mk.c_rng], 5.0) # 4 - (-1) + + # --- AML: explicit float('inf') in RangedExpression --- + m = pyo.ConcreteModel() + m.x = pyo.Var() + # One-sided: (-inf, x, 5) — lb should be treated as unbounded + m.c_ub = pyo.Constraint(rule=lambda m: (float('-inf'), m.x, 5)) + # Fully unbounded: (-inf, x, inf) — should be skipped entirely + m.c_skip = pyo.Constraint(rule=lambda m: (float('-inf'), m.x, float('inf'))) + + repn2 = LinearStandardFormCompiler().write( + m, mixed_form=True, keep_range_constraints=True + ) + by_con2 = {r.constraint: r.bound_type for r in repn2.rows} + self.assertEqual(by_con2[m.c_ub], 1) # ≤, not range + self.assertNotIn(m.c_skip, by_con2) # fully unbounded: skipped + + def test_extra_valid_ctypes(self): + """Component types in extra_valid_ctypes are permitted but not compiled.""" + m = pyo.ConcreteModel() + m.x = pyo.Var([1, 2, 3]) + m.y = pyo.Var() + m.obj = pyo.Objective(expr=m.y) + m.sos = pyo.SOSConstraint(var=m.x, sos=1) + + # Without extra_valid_ctypes, LSFC raises on the SOSConstraint. + with self.assertRaises(ValueError): + LinearStandardFormCompiler().write(m, mixed_form=True) + + # With extra_valid_ctypes, the SOSConstraint is silently skipped. + repn = LinearStandardFormCompiler().write( + m, mixed_form=True, extra_valid_ctypes=[pyo.SOSConstraint] + ) + # Only m.y appears in the objective; m.x[i] are unreferenced by + # linear constraints/objectives so not included in repn.columns. + self.assertEqual(len(repn.rows), 0) + col_ids = {id(v) for v in repn.columns} + self.assertIn(id(m.y), col_ids) + self.assertNotIn(id(m.x[1]), col_ids) + class TestTemplatedLinearStandardFormCompiler(TestLinearStandardFormCompiler): def setUp(self): diff --git a/pyomo/repn/tests/test_util.py b/pyomo/repn/tests/test_util.py index b8010d167f1..f26546a4186 100644 --- a/pyomo/repn/tests/test_util.py +++ b/pyomo/repn/tests/test_util.py @@ -878,6 +878,23 @@ def test_TemplateVarRecorder_user_varmap(self): vr.var_order, ) + def test_TemplateVarRecorder_kernel_variable(self): + """TemplateVarRecorder.add() must handle kernel variables, which do + not have a parent_component() method.""" + from pyomo.core.kernel.variable import variable + + v1 = variable() + v2 = variable() + + vm = {} + vr = TemplateVarRecorder(vm, SortComponents.deterministic) + vr.add(v1) + self.assertIn(id(v1), vm) + vr.add(v2) + self.assertIn(id(v2), vm) + self.assertEqual(len(vm), 2) + self.assertEqual({id(v1): 0, id(v2): 1}, vr.var_order) + if __name__ == "__main__": unittest.main() diff --git a/pyomo/repn/util.py b/pyomo/repn/util.py index 26d5ca20e1d..e80b162f29d 100644 --- a/pyomo/repn/util.py +++ b/pyomo/repn/util.py @@ -886,7 +886,12 @@ def add(self, var): # Note: the following is mostly a copy of # LinearBeforeChildDispatcher.record_var, but with extra # handling to update the env in the same loop - var_comp = var.parent_component() + try: + var_comp = var.parent_component() + except AttributeError: + # kernel variable objects do not have parent_component(); treat + # the variable itself as a scalar component. + var_comp = var # Double-check that the component has not already been processed # (through an individual var data) name = self.symbolmap.getSymbol(var_comp) @@ -903,9 +908,8 @@ def add(self, var): try: _iter = var_comp.items(self.sorter) except AttributeError: - # Note that this only works for the AML, as kernel does not - # provide a parent_component() - _iter = (None, var) + # kernel variables have no items(); treat as a scalar with index None + _iter = ((None, var),) if self._var_order is None: for i, (idx, v) in enumerate(_iter, start=len(vm)): vm[id(v)] = v diff --git a/pyomo/solvers/plugins/solvers/xpress_direct.py b/pyomo/solvers/plugins/solvers/xpress_direct.py index eacf67fd65f..df05efebb44 100644 --- a/pyomo/solvers/plugins/solvers/xpress_direct.py +++ b/pyomo/solvers/plugins/solvers/xpress_direct.py @@ -14,14 +14,21 @@ import time from pyomo.common.collections import ComponentSet, ComponentMap, Bunch -from pyomo.common.dependencies import attempt_import -from pyomo.common.errors import ApplicationError +from pyomo.common.dependencies import attempt_import, numpy as np +from pyomo.common.errors import ApplicationError, InvalidExpressionError from pyomo.common.tee import capture_output from pyomo.common.tempfiles import TempfileManager from pyomo.core.expr.numvalue import is_fixed from pyomo.core.expr.numvalue import value from pyomo.core.staleflag import StaleFlagManager from pyomo.repn import generate_standard_repn +from pyomo.repn.linear_template import LinearTemplateRepnVisitor +from pyomo.repn.plugins.standard_form import LinearStandardFormCompiler +from pyomo.repn.util import ( + FileDeterminism, + FileDeterminism_to_SortComponents, + TemplateVarRecorder, +) from pyomo.solvers.plugins.solvers.direct_solver import DirectSolver from pyomo.solvers.plugins.solvers.direct_or_persistent_solver import ( DirectOrPersistentSolver, @@ -33,6 +40,8 @@ from pyomo.opt.base import SolverFactory from pyomo.core.base.suffix import Suffix import pyomo.core.base.var +import pyomo.core.base.sos +import pyomo.core.base.constraint logger = logging.getLogger('pyomo.solvers') @@ -232,6 +241,15 @@ def _addConstraint( lambda self, prob, *args, **kwargs: prob.getIndexFromName(*args, **kwargs) ) XpressDirect._getObjIndex = lambda self, prob, obj: prob.getIndex(obj) + XpressDirect._addRows = lambda self, prob, rowtype, rhs, start, colind, rowcoef, rhsrange=None, names=None: prob.addrows( + rowtype, rhs, start, colind, rowcoef, range=rhsrange, names=names + ) + XpressDirect._chgObj = lambda self, prob, *args, **kwargs: prob.chgobj( + *args, **kwargs + ) + XpressDirect._chgObjSense = ( + lambda self, prob, *args, **kwargs: prob.chgobjsense(*args, **kwargs) + ) def _getLB(self, prob, *args, **kwargs): lb = [] @@ -281,6 +299,12 @@ def _getUB(self, prob, *args, **kwargs): XpressDirect._getUB = lambda self, prob, *args, **kwargs: prob.getUB( *args, **kwargs ) + XpressDirect._chgObj = lambda self, prob, *args, **kwargs: prob.chgObj( + *args, **kwargs + ) + XpressDirect._chgObjSense = ( + lambda self, prob, *args, **kwargs: prob.chgObjSense(*args, **kwargs) + ) def _addCols(self, prob, objx, mstart, mrwind, dmatval, bdl, bdu, names, types): first_col_ind = prob.attributes.cols @@ -294,6 +318,19 @@ def _addCols(self, prob, objx, mstart, mrwind, dmatval, bdl, bdu, names, types): XpressDirect._addCols = _addCols + def _addRows( + self, prob, rowtype, rhs, start, colind, rowcoef, rhsrange=None, names=None + ): + first_row_ind = prob.attributes.rows + prob.addRows( + rowtype, rhs, rng=rhsrange, start=start, colind=colind, rowcoef=rowcoef + ) + last_row_ind = prob.attributes.rows - 1 + if names is not None: + prob.addNames(xp.Namespaces.ROW, names, first_row_ind, last_row_ind) + + XpressDirect._addRows = _addRows + class _xpress_importer_class: # We want to be able to *update* the message that the deferred @@ -355,6 +392,14 @@ def __init__(self, **kwds): self._range_constraints = set() self._python_api_exists = xpress_available + # Tracks whether the current solver objective has quadratic terms, + # so _set_objective can decide whether to use the chgObj fast path. + self._obj_is_quadratic = False + + # Counter-based row-index cache (see _set_instance for details). + self._con_insertion_counter = 0 + self._con_name_to_counter = {} + self._deleted_counters = [] # TODO: this isn't a limit of XPRESS, which implements an SLP # method for NLPs. But it is a limit of *this* interface @@ -888,6 +933,14 @@ def _add_var(self, var): self._pyomo_var_to_solver_var_map[var] = xpress_var self._solver_var_to_pyomo_var_map[xpress_var] = var self._referenced_variables[var] = 0 + # Keep the template var_map in column order so that the + # LinearTemplateRepnVisitor produces Xpress column indices directly. + self._template_var_map[id(var)] = var + # Invalidate any cached template visitor (its var_recorder's var_order + # is derived lazily from _template_var_map, so new variables are picked + # up automatically; but the expanded_templates cache may reference stale + # column indices if the problem is rebuilt, so reset on structural change). + self._template_visitor = None def _set_instance(self, model, kwds={}): self._range_constraints = set() @@ -896,6 +949,22 @@ def _set_instance(self, model, kwds={}): self._solver_con_to_pyomo_con_map = ComponentMap() self._pyomo_var_to_solver_var_map = ComponentMap() self._solver_var_to_pyomo_var_map = ComponentMap() + # Counter-based row-index cache used by the persistent interface to + # delete rows in O(log M) without an O(N) getIndexFromName lookup. + # Each constraint gets a unique insertion counter; the current Xpress + # row index equals (counter - number of previously deleted rows with a + # lower counter), computable via bisect on the sorted deleted list. + self._con_insertion_counter = 0 + self._con_name_to_counter = {} + self._deleted_counters = [] + # var_map used by the template visitor: maps id(var) -> var in Xpress + # column order so that template-compiled column indices equal Xpress + # column indices. Reset here; populated in _add_var and + # _load_linear_problem. + self._template_var_map = {} + # Lazily constructed LinearTemplateRepnVisitor; created on first use + # of a templated constraint in _add_constraint. + self._template_visitor = None try: if model.name is not None: self._solver_model = xpress.problem(name=model.name) @@ -909,11 +978,260 @@ def _set_instance(self, model, kwds={}): "bindings for Xpress?\n\n\t" + "Error message: {0}".format(e) ) raise Exception(msg) - self._add_block(model) + self._load_linear_problem(model) def _add_block(self, block): DirectOrPersistentSolver._add_block(self, block) + def _get_template_visitor(self): + """Return (creating if needed) a LinearTemplateRepnVisitor seeded with + the current variable ordering so its column indices equal Xpress column + positions.""" + if self._template_visitor is None: + sorter = FileDeterminism_to_SortComponents(FileDeterminism.ORDERED) + var_recorder = TemplateVarRecorder(self._template_var_map, sorter) + self._template_visitor = LinearTemplateRepnVisitor( + {}, var_recorder=var_recorder + ) + return self._template_visitor + + def _load_linear_problem(self, model): + """Fast path for models with linear or mixed linear/nonlinear structure. + + Uses :class:`~pyomo.repn.plugins.standard_form.LinearStandardFormCompiler` + to build the full coefficient matrix for all linear constraints and + loads them into Xpress in a single ``loadproblem()`` call, bypassing + the per-constraint Python overhead of ``_add_block``. + + Any nonlinear constraints or objectives identified by the compiler are + added afterwards via the standard ``_add_constraint`` / ``_set_objective`` + per-expression path. + + All solver maps are populated so the persistent interface + (``add_constraint``, ``remove_constraint``, ``add_var``, etc.) continues + to work normally after this fast path. + + Parameters + ---------- + model : Pyomo ConcreteModel + """ + # Compile: mixed_form preserves <=, >=, == direction; + # keep_range_constraints emits range rows as a single 'R' entry. + # set_sense=None keeps objective coefficients as-is so we can set + # the sense ourselves via chgobjsense after loadproblem. + # allow_nonlinear=True: collect nonlinear constraints/objectives into + # repn.nonlinear_constraints / repn.nonlinear_objectives rather than + # raising; we handle them below via _add_constraint / _set_objective. + repn = LinearStandardFormCompiler().write( + model, + mixed_form=True, + keep_range_constraints=True, + set_sense=None, + allow_nonlinear=True, + extra_valid_ctypes=[pyomo.core.base.sos.SOSConstraint], + ) + + if len(repn.objectives) + len(repn.nonlinear_objectives) > 1: + raise ValueError("Solver interface does not support multiple objectives.") + + xprob = self._solver_model + n_cols = len(repn.columns) + + # ------------------------------------------------------------------ + # Objective + # ------------------------------------------------------------------ + if repn.objectives: + obj = repn.objectives[0] + if obj.sense == minimize: + obj_sense = xpress.minimize + elif obj.sense == maximize: + obj_sense = xpress.maximize + else: + raise ValueError( + 'Objective sense is not recognized: {0}'.format(obj.sense) + ) + # c is (n_objectives x n_cols) CSC; convert to CSR for row slicing. + c_csr = repn.c.tocsr() + c_dense = np.asarray(c_csr[0, :].todense()).flatten() + c_offset = float(repn.c_offset[0]) + else: + obj = None + obj_sense = xpress.minimize + c_csr = None + c_dense = np.zeros(n_cols) + c_offset = 0.0 + + # ------------------------------------------------------------------ + # Variable arrays + # ------------------------------------------------------------------ + lb = [] + ub = [] + colnames = [] + entind = [] + coltype_chars = [] + for i, var in enumerate(repn.columns): + bnd_lb, bnd_ub = var.bounds + lb.append(-xpress.infinity if bnd_lb is None else float(bnd_lb)) + ub.append(xpress.infinity if bnd_ub is None else float(bnd_ub)) + colnames.append(self._symbol_map.getSymbol(var, self._labeler)) + if var.is_binary(): + entind.append(i) + coltype_chars.append('B') + elif var.is_integer(): + entind.append(i) + coltype_chars.append('I') + + # ------------------------------------------------------------------ + # Constraint arrays + # ------------------------------------------------------------------ + _rowtype_map = {0: 'E', 1: 'L', -1: 'G', 2: 'R'} + rowtype = [] + rownames = [] + for row_entry in repn.rows: + rowtype.append(_rowtype_map[row_entry.bound_type]) + rownames.append( + self._symbol_map.getSymbol(row_entry.constraint, self._labeler) + ) + + rhs = repn.rhs + rng = repn.rhs_range + + # ------------------------------------------------------------------ + # CSC coefficient matrix + # repn.A is already CSC; pass the numpy arrays directly — loadproblem + # accepts them without requiring a Python-list copy. + # ------------------------------------------------------------------ + A = repn.A.asformat('csc') + + # ------------------------------------------------------------------ + # Load the full LP/MIP in one call + # ------------------------------------------------------------------ + xprob.loadproblem( + model.name or '', + rowtype, + rhs, + rng, + c_dense, + A.indptr, + None, # collen — use indptr[ncol] convention + A.indices, + A.data, + lb, + ub, + colnames=colnames, + rownames=rownames, + coltype=coltype_chars if coltype_chars else None, + entind=entind if entind else None, + ) + + # Objective sense and constant term (OBJRHS). + # chgobj([-1], [-k]) sets OBJRHS = k so the solver adds +k to c^T x. + self._chgObjSense(xprob, obj_sense) + if c_offset != 0.0: + self._chgObj(xprob, [-1], [-c_offset]) + self._obj_is_quadratic = False + + # ------------------------------------------------------------------ + # Retrieve linked xpress.var objects (in column index order) + # ------------------------------------------------------------------ + xpress_vars = xprob.getVariable() + + # ------------------------------------------------------------------ + # Populate variable maps + # ------------------------------------------------------------------ + for var, xv in zip(repn.columns, xpress_vars): + self._pyomo_var_to_solver_var_map[var] = xv + self._solver_var_to_pyomo_var_map[xv] = var + self._referenced_variables[var] = 0 + + # Build _template_var_map in Xpress column order so the + # LinearTemplateRepnVisitor produces column indices that directly + # correspond to Xpress column positions. + self._template_var_map = {id(var): var for var in repn.columns} + self._template_visitor = None + + # LSFC only includes variables that appear in linear constraints or + # objectives. Variables that are unreferenced or appear only in SOS + # constraints must still be added to Xpress so the solver model is + # complete. We detect them by checking against _pyomo_var_to_solver_var_map + # and fall back to _add_var for each missing one. + for var in model.component_data_objects( + pyomo.core.base.var.Var, active=True, descend_into=True + ): + if var not in self._pyomo_var_to_solver_var_map: + self._add_var(var) + + # ------------------------------------------------------------------ + # Populate constraint maps and O(log M) row-index counters + # ------------------------------------------------------------------ + for i, (row_entry, conname) in enumerate(zip(repn.rows, rownames)): + con = row_entry.constraint + self._pyomo_con_to_solver_con_map[con] = conname + self._solver_con_to_pyomo_con_map[conname] = con + self._con_name_to_counter[conname] = i + if row_entry.bound_type == 2: + self._range_constraints.add(conname) + self._con_insertion_counter = len(repn.rows) + + # ------------------------------------------------------------------ + # Populate _vars_referenced_by_con and _referenced_variables + # ------------------------------------------------------------------ + # Convert A to CSR for efficient per-row nonzero-column access. + A_csr = A.tocsr() + for i, row_entry in enumerate(repn.rows): + con = row_entry.constraint + j0, j1 = int(A_csr.indptr[i]), int(A_csr.indptr[i + 1]) + vars_in_con = ComponentSet(repn.columns[j] for j in A_csr.indices[j0:j1]) + self._vars_referenced_by_con[con] = vars_in_con + for var in vars_in_con: + self._referenced_variables[var] += 1 + + # ------------------------------------------------------------------ + # Populate objective reference state + # ------------------------------------------------------------------ + if obj is not None: + j0, j1 = int(c_csr.indptr[0]), int(c_csr.indptr[1]) + obj_vars = ComponentSet(repn.columns[j] for j in c_csr.indices[j0:j1]) + for var in obj_vars: + self._referenced_variables[var] += 1 + self._objective = obj + self._vars_referenced_by_obj = obj_vars + + # ------------------------------------------------------------------ + # Add SOS constraints (variable maps must be populated first) + # ------------------------------------------------------------------ + for sub_block in model.block_data_objects(descend_into=True, active=True): + for sos_con in sub_block.component_data_objects( + ctype=pyomo.core.base.sos.SOSConstraint, + descend_into=False, + active=True, + sort=True, + ): + self._add_sos_constraint(sos_con) + + # ------------------------------------------------------------------ + # Nonlinear constraints and objectives + # ------------------------------------------------------------------ + # Any constraints/objectives the LSFC could not linearize are added + # via the standard per-expression path, which supports quadratic and + # general nonlinear expressions through Xpress' expression-tree API. + for con in repn.nonlinear_constraints: + self._add_constraint(con) + if repn.nonlinear_objectives: + self._set_objective(repn.nonlinear_objectives[0]) + + # LSFC silently skips trivially feasible constant constraints (e.g., + # `0 <= 1 <= 3`). When _skip_trivial_constraints is False we must + # still register them so they appear in the symbol map. + if not self._skip_trivial_constraints: + for con in model.component_data_objects( + pyomo.core.base.constraint.Constraint, active=True, descend_into=True + ): + if con not in self._pyomo_con_to_solver_con_map and ( + con.has_lb() or con.has_ub() + ): + self._add_constraint(con) + def _add_constraint(self, con): if not con.active: return None @@ -923,15 +1241,6 @@ def _add_constraint(self, con): conname = self._symbol_map.getSymbol(con, self._labeler) - if con._linear_canonical_form: - xpress_expr, referenced_vars = self._get_expr_from_pyomo_repn( - con.canonical_form(), self._max_constraint_degree - ) - else: - xpress_expr, referenced_vars = self._get_expr_from_pyomo_expr( - con.body, self._max_constraint_degree - ) - if con.has_lb(): if not is_fixed(con.lower): raise ValueError( @@ -943,51 +1252,200 @@ def _add_constraint(self, con): "Upper bound of constraint {0} is not constant.".format(con) ) - if con.equality: - xpress_con = self._addConstraint( - self._solver_model, - body=xpress_expr, - type=xpress.eq, - rhs=value(con.lower), - name=conname, - ) - elif con.has_lb() and con.has_ub(): - xpress_con = self._addConstraint( - self._solver_model, - body=xpress_expr, - type=xpress.rng, - lb=value(con.lower), - ub=value(con.upper), - name=conname, - ) - self._range_constraints.add(xpress_con) - elif con.has_lb(): - xpress_con = self._addConstraint( - self._solver_model, - body=xpress_expr, - type=xpress.geq, - rhs=value(con.lower), - name=conname, - ) - elif con.has_ub(): - xpress_con = self._addConstraint( - self._solver_model, - body=xpress_expr, - type=xpress.leq, - rhs=value(con.upper), - name=conname, + # ------------------------------------------------------------------ + # Fast path: templated linear constraint + # ------------------------------------------------------------------ + # TemplateConstraintData objects carry a pre-compiled body lambda. + # Evaluating it is much cheaper than walking the expression tree via + # generate_standard_repn, and the column indices it returns are + # exactly the Xpress column positions (because _template_var_map is + # kept in column order). + if hasattr(con, 'template_expr'): + try: + tv = self._get_template_visitor() + offset, col_indices, coeffs, lb, ub = tv.expand_expression( + con, con.template_expr() + ) + except InvalidExpressionError: + pass # nonlinear template — fall through to generate_standard_repn + else: + n = len(col_indices) + referenced_vars = ComponentSet( + tv.var_recorder.var_list[i] for i in col_indices + ) + coeffs = [float(c) for c in coeffs] + const = float(offset) + + if lb == ub and lb is not None: + rowtype = ['E'] + rhs = [lb - const] + rhsrange = None + elif lb is not None and ub is not None: + ub_val = ub - const + lb_val = lb - const + rowtype = ['R'] + rhs = [ub_val] + rhsrange = [ub_val - lb_val] + self._range_constraints.add(conname) + elif lb is not None: + rowtype = ['G'] + rhs = [lb - const] + rhsrange = None + elif ub is not None: + rowtype = ['L'] + rhs = [ub - const] + rhsrange = None + else: + raise ValueError( + "Constraint does not have a lower " + "or an upper bound: {0} \n".format(con) + ) + + self._addRows( + self._solver_model, + rowtype, + rhs, + [0, n], + col_indices, + coeffs, + rhsrange=rhsrange, + names=[conname], + ) + for var in referenced_vars: + self._referenced_variables[var] += 1 + self._vars_referenced_by_con[con] = referenced_vars + self._pyomo_con_to_solver_con_map[con] = conname + self._solver_con_to_pyomo_con_map[conname] = con + self._con_name_to_counter[conname] = self._con_insertion_counter + self._con_insertion_counter += 1 + return + + # ------------------------------------------------------------------ + # General path: generate_standard_repn (linear, quadratic, nonlinear) + # ------------------------------------------------------------------ + if con._linear_canonical_form: + repn = con.canonical_form() + else: + repn = generate_standard_repn(con.body, quadratic=True) + + degree = repn.polynomial_degree() + if (degree is None) or (degree > self._max_constraint_degree): + raise DegreeError( + 'XpressDirect does not support expressions of degree {0}.' + '\nexpr: {1}'.format(degree, con.body) ) + + if repn.quadratic_vars: + # Quadratic constraint: build an Xpress expression tree and use + # _addConstraint, which links a Python constraint object. + try: + xpress_expr, referenced_vars = self._get_expr_from_pyomo_repn( + repn, self._max_constraint_degree + ) + except DegreeError as e: + msg = e.args[0] + msg += '\nexpr: {0}'.format(con.body) + raise DegreeError(msg) + + if con.equality: + self._addConstraint( + self._solver_model, + body=xpress_expr, + type=xpress.eq, + rhs=value(con.lower), + name=conname, + ) + elif con.has_lb() and con.has_ub(): + self._addConstraint( + self._solver_model, + body=xpress_expr, + type=xpress.rng, + lb=value(con.lower), + ub=value(con.upper), + name=conname, + ) + self._range_constraints.add(conname) + elif con.has_lb(): + self._addConstraint( + self._solver_model, + body=xpress_expr, + type=xpress.geq, + rhs=value(con.lower), + name=conname, + ) + elif con.has_ub(): + self._addConstraint( + self._solver_model, + body=xpress_expr, + type=xpress.leq, + rhs=value(con.upper), + name=conname, + ) + else: + raise ValueError( + "Constraint does not have a lower " + "or an upper bound: {0} \n".format(con) + ) else: - raise ValueError( - "Constraint does not have a lower " - "or an upper bound: {0} \n".format(con) + # Linear constraint: bypass expression-tree construction and add + # the row directly using _addRows, which is significantly faster + # for large models. + referenced_vars = ComponentSet(repn.linear_vars) + xpress_vars = [ + self._pyomo_var_to_solver_var_map[v] for v in repn.linear_vars + ] + coeffs = [float(c) for c in repn.linear_coefs] + const = float(repn.constant) + n = len(xpress_vars) + + if con.equality: + rowtype = ['E'] + rhs = [value(con.lower) - const] + rhsrange = None + elif con.has_lb() and con.has_ub(): + ub_val = value(con.upper) - const + lb_val = value(con.lower) - const + rowtype = ['R'] + rhs = [ub_val] + rhsrange = [ub_val - lb_val] + self._range_constraints.add(conname) + elif con.has_lb(): + rowtype = ['G'] + rhs = [value(con.lower) - const] + rhsrange = None + elif con.has_ub(): + rowtype = ['L'] + rhs = [value(con.upper) - const] + rhsrange = None + else: + raise ValueError( + "Constraint does not have a lower " + "or an upper bound: {0} \n".format(con) + ) + + self._addRows( + self._solver_model, + rowtype, + rhs, + [0, n], + xpress_vars, + coeffs, + rhsrange=rhsrange, + names=[conname], ) for var in referenced_vars: self._referenced_variables[var] += 1 self._vars_referenced_by_con[con] = referenced_vars - self._pyomo_con_to_solver_con_map[con] = xpress_con - self._solver_con_to_pyomo_con_map[xpress_con] = con + # Store constraint name string (not the xp.constraint object) so the + # solver maps are stable across row additions/deletions and compatible + # with the name-based getDuals/getSlacks APIs. + self._pyomo_con_to_solver_con_map[con] = conname + self._solver_con_to_pyomo_con_map[conname] = con + # Record the insertion counter so the persistent interface can compute + # the current Xpress row index in O(log M) without getIndexFromName. + self._con_name_to_counter[conname] = self._con_insertion_counter + self._con_insertion_counter += 1 def _add_sos_constraint(self, con): if not con.active: @@ -1045,12 +1503,6 @@ def _xpress_vartype_from_var(self, var): return vartype def _set_objective(self, obj): - if self._objective is not None: - for var in self._vars_referenced_by_obj: - self._referenced_variables[var] -= 1 - self._vars_referenced_by_obj = ComponentSet() - self._objective = None - if obj.active is False: raise ValueError('Cannot add inactive objective to solver.') @@ -1061,14 +1513,79 @@ def _set_objective(self, obj): else: raise ValueError('Objective sense is not recognized: {0}'.format(obj.sense)) - xpress_expr, referenced_vars = self._get_expr_from_pyomo_expr( - obj.expr, self._max_obj_degree - ) + repn = generate_standard_repn(obj.expr, quadratic=True) - for var in referenced_vars: - self._referenced_variables[var] += 1 + degree = repn.polynomial_degree() + if (degree is None) or (degree > self._max_obj_degree): + raise DegreeError( + 'XpressDirect does not support expressions of degree {0}.'.format( + degree + ) + ) + + if repn.quadratic_vars or self._obj_is_quadratic: + # Quadratic objective (or replacing a quadratic one): use + # setObjective which handles both the linear and quadratic parts + # and safely clears any previous quadratic terms. + try: + xpress_expr, referenced_vars = self._get_expr_from_pyomo_repn( + repn, self._max_obj_degree + ) + except DegreeError as e: + msg = e.args[0] + msg += '\nexpr: {0}'.format(obj.expr) + raise DegreeError(msg) + + if self._objective is not None: + for var in self._vars_referenced_by_obj: + self._referenced_variables[var] -= 1 + self._vars_referenced_by_obj = ComponentSet() + self._objective = None + + for var in referenced_vars: + self._referenced_variables[var] += 1 + + self._solver_model.setObjective(xpress_expr, sense=sense) + self._obj_is_quadratic = bool(repn.quadratic_vars) + else: + # Linear objective fast path: build coefficient arrays directly + # and call chgObj, bypassing Xpress expression-tree construction. + referenced_vars = ComponentSet(repn.linear_vars) + new_xpress_vars = [ + self._pyomo_var_to_solver_var_map[v] for v in repn.linear_vars + ] + new_coeffs = [float(c) for c in repn.linear_coefs] + const = float(repn.constant) + + if self._objective is not None: + # Zero out objective coefficients for variables that are no + # longer in the new objective. + vars_to_zero = self._vars_referenced_by_obj - referenced_vars + if vars_to_zero: + zero_xpress_vars = [ + self._pyomo_var_to_solver_var_map[v] for v in vars_to_zero + ] + self._chgObj( + self._solver_model, + zero_xpress_vars, + [0.0] * len(zero_xpress_vars), + ) + for var in self._vars_referenced_by_obj: + self._referenced_variables[var] -= 1 + self._vars_referenced_by_obj = ComponentSet() + self._objective = None + + # chgobj([-1], [k]) sets OBJRHS = -k, so the objective constant + # becomes +k. Pass -const so that OBJRHS = const. + self._chgObj( + self._solver_model, new_xpress_vars + [-1], new_coeffs + [-const] + ) + self._chgObjSense(self._solver_model, sense) + self._obj_is_quadratic = False + + for var in referenced_vars: + self._referenced_variables[var] += 1 - self._solver_model.setObjective(xpress_expr, sense=sense) self._objective = obj self._vars_referenced_by_obj = referenced_vars @@ -1167,9 +1684,15 @@ def _postsolve(self): soln_constraints = soln.constraint if extract_duals or extract_slacks: - xpress_cons = list(self._solver_con_to_pyomo_con_map.keys()) + # Only include regular (non-SOS) constraints, which are + # keyed by name strings in the solver-con map. + xpress_cons = [ + con + for con in self._solver_con_to_pyomo_con_map.keys() + if isinstance(con, str) + ] for con in xpress_cons: - soln_constraints[con.name] = {} + soln_constraints[con] = {} xpress_vars = list(self._solver_var_to_pyomo_var_map.keys()) try: @@ -1209,24 +1732,25 @@ def _postsolve(self): if extract_duals: vals = self._getDuals(xprob, xpress_cons) for val, con in zip(vals, xpress_cons): - soln_constraints[con.name]["Dual"] = val + soln_constraints[con]["Dual"] = val if extract_slacks: for con, val in zip(xpress_cons, slacks): if con in self._range_constraints: ## for xpress, the slack on a range constraint ## is based on the upper bound - lb = con.lb - ub = con.ub + pyomo_con = self._solver_con_to_pyomo_con_map[con] + lb = value(pyomo_con.lower) + ub = value(pyomo_con.upper) ub_s = val expr_val = ub - ub_s lb_s = lb - expr_val if abs(ub_s) > abs(lb_s): - soln_constraints[con.name]["Slack"] = ub_s + soln_constraints[con]["Slack"] = ub_s else: - soln_constraints[con.name]["Slack"] = lb_s + soln_constraints[con]["Slack"] = lb_s else: - soln_constraints[con.name]["Slack"] = val + soln_constraints[con]["Slack"] = val elif self._load_solutions: if have_soln: @@ -1321,8 +1845,8 @@ def _load_slacks(self, cons_to_load=None): if xpress_con in self._range_constraints: ## for xpress, the slack on a range constraint ## is based on the upper bound - lb = xpress_con.lb - ub = xpress_con.ub + lb = value(pyomo_con.lower) + ub = value(pyomo_con.upper) ub_s = val expr_val = ub - ub_s lb_s = lb - expr_val diff --git a/pyomo/solvers/plugins/solvers/xpress_persistent.py b/pyomo/solvers/plugins/solvers/xpress_persistent.py index ea103901ab0..4fa3ce1c188 100644 --- a/pyomo/solvers/plugins/solvers/xpress_persistent.py +++ b/pyomo/solvers/plugins/solvers/xpress_persistent.py @@ -13,7 +13,7 @@ from pyomo.core.expr.numvalue import value, is_fixed import pyomo.core.expr as EXPR from pyomo.opt.base import SolverFactory -import collections +import bisect @SolverFactory.register( @@ -50,12 +50,29 @@ def __init__(self, **kwds): self.set_instance(self._pyomo_model, **kwds) def _remove_constraint(self, solver_con): - self._solver_model.delConstraint(solver_con) + # solver_con is the constraint name string. Use the insertion-counter + # cache to compute the current Xpress row index in O(log M) — avoiding + # an O(N) getIndexFromName scan — then delete the row and record the + # deletion so future index lookups stay accurate. + counter = self._con_name_to_counter.pop(solver_con) + current_idx = counter - bisect.bisect_left(self._deleted_counters, counter) + self._solver_model.delConstraint(current_idx) + bisect.insort(self._deleted_counters, counter) def _remove_sos_constraint(self, solver_sos_con): self._solver_model.delSOS(solver_sos_con) def _remove_var(self, solver_var): + # Keep _template_var_map in sync: removing a column shifts all + # subsequent Xpress column indices down by one. Deleting the entry + # from the ordered dict has exactly the same effect on the remaining + # vars' positions (from enumerate), so their indices stay correct + # after the visitor is rebuilt. Reset _template_visitor so + # TemplateVarRecorder.var_order is recomputed on next use. + pyomo_var = self._solver_var_to_pyomo_var_map.get(solver_var) + if pyomo_var is not None: + self._template_var_map.pop(id(pyomo_var), None) + self._template_visitor = None self._solver_model.delVariable(solver_var) def _warm_start(self): @@ -155,6 +172,10 @@ def _add_column(self, var, obj_coef, constraints, coefficients): self._pyomo_var_to_solver_var_map[var] = xpress_var self._solver_var_to_pyomo_var_map[xpress_var] = var self._referenced_variables[var] = len(coefficients) + # Append to _template_var_map in column order (new column is last) + # and invalidate any cached template visitor. + self._template_var_map[id(var)] = var + self._template_visitor = None def get_xpress_attribute(self, *args): """