diff --git a/pyomo/repn/plugins/standard_form.py b/pyomo/repn/plugins/standard_form.py index 70b1b74b8cc..162fe86271d 100644 --- a/pyomo/repn/plugins/standard_form.py +++ b/pyomo/repn/plugins/standard_form.py @@ -58,6 +58,9 @@ RowEntry = collections.namedtuple('RowEntry', ['constraint', 'bound_type']) +_inf = float('inf') +_ninf = -_inf + # TODO: make a proper base class class LinearStandardFormInfo: @@ -82,15 +85,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] @@ -113,11 +132,14 @@ class LinearStandardFormInfo: """ - 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 + ): 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 @@ -180,6 +202,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( @@ -412,10 +447,16 @@ def write(self, model): # 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 = [] @@ -453,6 +494,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 == _ninf: + lb = None + if ub == _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 +522,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 +544,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 +553,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 +589,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 +693,15 @@ 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, ) 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..80114b222f5 100644 --- a/pyomo/repn/tests/test_standard_form.py +++ b/pyomo/repn/tests/test_standard_form.py @@ -351,6 +351,116 @@ 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_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. + """ + # --- 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_kernel_compatibility(self): + 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) + 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 7289f922532..019230dfe29 100644 --- a/pyomo/repn/util.py +++ b/pyomo/repn/util.py @@ -891,7 +891,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) @@ -908,9 +913,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