From 0dd64c75155fa5b1761c5dfa4cb1e20ea1e70943 Mon Sep 17 00:00:00 2001 From: Doug A Date: Fri, 10 Jul 2026 14:59:37 -0400 Subject: [PATCH 1/6] clipping and scaling --- pyomo/util/calc_var_value.py | 84 +++++++++++++++++++++++-- pyomo/util/tests/test_calc_var_value.py | 54 +++++++++++++++- 2 files changed, 131 insertions(+), 7 deletions(-) diff --git a/pyomo/util/calc_var_value.py b/pyomo/util/calc_var_value.py index c9a58b68f6e..c3cd71f96c4 100644 --- a/pyomo/util/calc_var_value.py +++ b/pyomo/util/calc_var_value.py @@ -11,6 +11,7 @@ from pyomo.common.numeric_types import native_numeric_types, native_complex_types, value from pyomo.core.expr.calculus.derivatives import differentiate from pyomo.core.base.constraint import Constraint +from pyomo.core.base.suffix import SuffixFinder import logging @@ -24,6 +25,53 @@ } +def _clip_variable_to_bounds(variable, expr, eps): + """ + Function to try to clip a variable back within its bounds. + If the expression residual is less than eps, then the new + value of the variable is retained. If the value is greater + than eps, then the original, bound-violating value of + variable is retained. + + Parameters: + ----------- + variable: :py:class:`VarData` + The variable to clip to within its bounds + expr: :py:class:`NumericExpression` + The expression to evaluate the constraint residual + eps: `float` + The tolerance to use to determine constraint feasibility. + + Returns: + -------- + None + + """ + # Should we check against variable domain here as well? + variable_value = value(variable) + if ( + variable.lb is not None + and variable.ub is not None + and variable.lb > variable.ub + ): + raise ValueError( + f"Variable {variable} has inconsistent bounds. " + f"It has a lower bound of {variable.lb} and an " + f"upper bound of {variable.ub}." + ) + elif variable.lb is not None and variable_value < variable.lb: + variable.set_value(variable.lb) + if abs(value(expr)) < eps: + return + # Else proceed to set the variable to its old value + elif variable.ub is not None and variable_value > variable.ub: + variable.set_value(variable.ub) + if abs(value(expr)) < eps: + return + # Else proceed to set the variable to its old value + variable.set_value(variable_value) + + def calculate_variable_from_constraint( variable, constraint, @@ -32,6 +80,8 @@ def calculate_variable_from_constraint( linesearch=True, alpha_min=1e-8, diff_mode=None, + scale_problem=False, + clip_to_bounds=False, ): """Calculate the variable value given a specified equality constraint @@ -72,6 +122,16 @@ def calculate_variable_from_constraint( diff_mode: :py:enum:`pyomo.core.expr.calculus.derivatives.Modes` The mode to use to differentiate the expression. If unspecified, defaults to `Modes.sympy` + scale_problem: `bool` + Whether to take variable and constraint scaling into account + when calculating if the constraint residual tolerance is + satisfied or if the initial derivative is zero. + clip_to_bounds: `bool` + If terminating at a point outside a variable's bounds, set the + variable to its closest bound and see whether the constraint + residual termination condition is satisfied. This method is + applied only to constraints that are not solved by the shortcut + methods used to quickly solve linear constraints. Returns: -------- @@ -97,6 +157,16 @@ def calculate_variable_from_constraint( if lower != upper: raise ValueError(f"Constraint '{constraint}' must be an equality constraint") + if scale_problem: + finder = SuffixFinder('scaling_factor', 1.0, constraint.model()) + + sf_variable = finder.find(variable) + sf_constraint = finder.find(constraint) + eps /= sf_constraint + else: + sf_variable = 1 + sf_constraint = 1 + _invalid_types = set(native_complex_types) _invalid_types.add(type(None)) @@ -219,7 +289,7 @@ def calculate_variable_from_constraint( else: fp0 = value(expr_deriv) - if abs(value(fp0)) < 1e-12: + if abs(value(fp0) * sf_constraint / sf_variable) < 1e-12: raise ValueError( f"Initial value for variable '{variable}' results in a derivative " f"value for constraint '{constraint}' that is very close to zero.\n" @@ -258,7 +328,7 @@ def calculate_variable_from_constraint( else: fpk = value(expr_deriv) - if abs(fpk) < 1e-12: + if abs(fpk * sf_constraint / sf_variable) < 1e-12: # TODO: should this raise a ValueError or a new # DerivativeError (subclassing ArithmeticError)? raise RuntimeError( @@ -310,7 +380,9 @@ def calculate_variable_from_constraint( f"variable '{variable}' using constraint '{constraint}'; " f"remaining residual = {residual}." ) - # - # Re-set the variable value to trigger any warnings WRT the final - # variable state - variable.set_value(variable.value) + if clip_to_bounds: + _clip_variable_to_bounds(variable, expr, eps) + else: + # Re-set the variable value to trigger any warnings WRT the final + # variable state + variable.set_value(variable.value) diff --git a/pyomo/util/tests/test_calc_var_value.py b/pyomo/util/tests/test_calc_var_value.py index e1f7b5bb2f2..756422afb5d 100644 --- a/pyomo/util/tests/test_calc_var_value.py +++ b/pyomo/util/tests/test_calc_var_value.py @@ -24,6 +24,7 @@ exp, NonNegativeReals, Binary, + Suffix, ) from pyomo.util.calc_var_value import calculate_variable_from_constraint from pyomo.core.expr.calculus.diff_with_sympy import differentiate_available @@ -240,7 +241,26 @@ def test_nonlinear(self): m.x, m.f, linesearch=False, diff_mode=mode ) - # Calculate the bubble point of Benzene. THe first step + # Test whether scaling is being correctly applied + m.scaling_factor = Suffix(direction=Suffix.EXPORT) + m.scaling_factor[m.e] = 1e4 + for mode in all_diff_modes: + m.x.set_value(3.1) + calculate_variable_from_constraint( + m.x, m.e, eps=1e-4, diff_mode=mode, scale_problem=False + ) + # Without scaling, we should stop at x = 3 + 1.028e-5 + assert value(m.x) - 3 > 1e-5 + + for mode in all_diff_modes: + m.x.set_value(3.1) + calculate_variable_from_constraint( + m.x, m.e, eps=1e-6, diff_mode=mode, scale_problem=True + ) + # With scaling, we should stop at x = 3 + 5.288-11 + self.assertAlmostEqual(value(m.x), 3, places=9) + + # Calculate the bubble point of Benzene. The first step # computed by calculate_variable_from_constraint will make the # second term become complex, and the evaluation will fail. # This tests that the algorithm cleanly continues @@ -413,6 +433,34 @@ def test_warn_final_value_nonlinear(self): ) self.assertAlmostEqual(value(m.x), 3.5, 3) + @unittest.skipUnless(differentiate_available, "this test requires sympy") + def test_clip_bounds_nonlinear(self): + # This test was written with the assistance of Google Gemini 3.5 Flash + # (in order to generate the regex to match the value in the log). + m = ConcreteModel() + # The first step of Newton's method will overshoot the root + # and subsequent steps will be negative. + m.x = Var(initialize=1, bounds=(0, None)) + m.c = Constraint(expr=-1.5 * m.x**2 + 3.5 * m.x == 0) + + with LoggingIntercept() as LOG: + calculate_variable_from_constraint(m.x, m.c, eps=1e-6, linesearch=False) + self.assertRegex( + LOG.getvalue().strip(), + r"Setting Var 'x' to a numeric value `-?(?:\d+(?:\.\d*)?|\.\d+)[eE][+-]?\d+` outside the " + r"bounds \(0, None\).", + ) + self.assertAlmostEqual(value(m.x), 0) + assert value(m.x) < 0 + + m.x.set_value(1) + with LoggingIntercept() as LOG: + calculate_variable_from_constraint( + m.x, m.c, eps=1e-6, linesearch=False, clip_to_bounds=True + ) + assert len(LOG.getvalue().strip()) == 0 + assert value(m.x) == 0 + @unittest.skipUnless(differentiate_available, "this test requires sympy") def test_nonlinear_overflow(self): # Regression check to make sure calculate_variable_from_constraint @@ -466,3 +514,7 @@ def test_linea_search_overflow(self): calculate_variable_from_constraint(m.x, m.con) self.assertAlmostEqual(value(m.x), 0) + + +if __name__ == "__main__": + unittest.main() From 69af65d9be36eb7dad868c756e6fdc3b3eec0b0c Mon Sep 17 00:00:00 2001 From: Doug A Date: Mon, 13 Jul 2026 10:52:34 -0400 Subject: [PATCH 2/6] more tests --- pyomo/util/calc_var_value.py | 23 ++++++--- pyomo/util/tests/test_calc_var_value.py | 65 +++++++++++++++++++++++-- 2 files changed, 79 insertions(+), 9 deletions(-) diff --git a/pyomo/util/calc_var_value.py b/pyomo/util/calc_var_value.py index c3cd71f96c4..b6855307c58 100644 --- a/pyomo/util/calc_var_value.py +++ b/pyomo/util/calc_var_value.py @@ -301,11 +301,18 @@ def calculate_variable_from_constraint( while abs(fk) > eps and iter_left: iter_left -= 1 if not iter_left: - raise IterationLimitError( - f"Iteration limit (%s) reached solving for variable '{variable}' " - f"using constraint '{constraint}'; remaining residual = %s" - % (iterlim, value(expr)) - ) + if scale_problem: + raise IterationLimitError( + f"Iteration limit (%s) reached solving for variable '{variable}' " + f"using constraint '{constraint}'; remaining (scaled) residual = %s" + % (iterlim, sf_constraint * value(expr)) + ) + else: + raise IterationLimitError( + f"Iteration limit (%s) reached solving for variable '{variable}' " + f"using constraint '{constraint}'; remaining residual = %s" + % (iterlim, value(expr)) + ) # compute step xk = value(variable) @@ -356,7 +363,11 @@ def calculate_variable_from_constraint( if fkp1.__class__ in _invalid_types: # We cannot perform computations on complex numbers fkp1 = None - if fkp1 is not None and fkp1**2 < c1 * fk**2: + if ( + fkp1 is not None + and (fkp1 * sf_constraint / sf_variable) ** 2 + < c1 * (fk * sf_constraint / sf_variable) ** 2 + ): # found an alpha value with sufficient reduction # continue to the next step fk = fkp1 diff --git a/pyomo/util/tests/test_calc_var_value.py b/pyomo/util/tests/test_calc_var_value.py index 756422afb5d..5d3dbea8406 100644 --- a/pyomo/util/tests/test_calc_var_value.py +++ b/pyomo/util/tests/test_calc_var_value.py @@ -37,6 +37,10 @@ differentiate.Modes.reverse_numeric, ] +# This regex was written with Google Gemini 3.5 to match an arbitrary +# floating point number in scientific notation +SCIENTIFIC_NOTATION_REGEX = r"[+-]?(?:\d+(?:\.\d*)?)[eE][+-]?\d+" + def sum_sq(args, fixed, fgh): f = sum(arg**2 for arg in args) @@ -260,6 +264,60 @@ def test_nonlinear(self): # With scaling, we should stop at x = 3 + 5.288-11 self.assertAlmostEqual(value(m.x), 3, places=9) + # Check whether scaling factors are being used when calculating + # problem derivatives + # Here, the variable scaling causes a failure because the scaled + # derivative is near-zero + m.rhs = Param(initialize=3, mutable=True) + m.g = Constraint(expr=m.x**2 + m.x == m.rhs) + m.scaling_factor[m.x] = 1e16 + for mode in all_diff_modes: + m.x.set_value(1) + with self.assertRaisesRegex( + ValueError, + "Initial value for variable 'x' results in a " + "derivative value for constraint 'g' that is very close to zero.", + ): + calculate_variable_from_constraint( + m.x, m.g, diff_mode=mode, scale_problem=True + ) + + # Here, the constraint scaling factor causes a failure because the + # scaled derivative is near-zero + del m.scaling_factor[m.x] + # Set rhs to large value to make sure we don't terminate due to a small + # constraint residual in any of the shortcut steps + m.rhs.set_value(1e16) + m.scaling_factor[m.g] = 1e-16 + for mode in all_diff_modes: + m.x.set_value(1) + with self.assertRaisesRegex( + ValueError, + "Initial value for variable 'x' results in a " + "derivative value for constraint 'g' that is very close to zero.", + ): + calculate_variable_from_constraint( + m.x, m.g, diff_mode=mode, scale_problem=True + ) + # Here, because both the variable and constraint are scaled, we don't terminate due to + # a zero derivative. It still raises an error, however, because we're asking for an + # unreasonable amount of precision + m.scaling_factor[m.x] = 1e16 + m.scaling_factor[m.g] = 1e16 + m.rhs.set_value(3) + for mode in all_diff_modes: + m.x.set_value(1) + with self.assertRaisesRegex( + IterationLimitError, + "Linesearch iteration limit reached solving for variable 'x' using " + "constraint 'g'; remaining residual = " + + SCIENTIFIC_NOTATION_REGEX + + r"\.", + ): + calculate_variable_from_constraint( + m.x, m.g, diff_mode=mode, scale_problem=True + ) + # Calculate the bubble point of Benzene. The first step # computed by calculate_variable_from_constraint will make the # second term become complex, and the evaluation will fail. @@ -435,8 +493,7 @@ def test_warn_final_value_nonlinear(self): @unittest.skipUnless(differentiate_available, "this test requires sympy") def test_clip_bounds_nonlinear(self): - # This test was written with the assistance of Google Gemini 3.5 Flash - # (in order to generate the regex to match the value in the log). + m = ConcreteModel() # The first step of Newton's method will overshoot the root # and subsequent steps will be negative. @@ -447,7 +504,9 @@ def test_clip_bounds_nonlinear(self): calculate_variable_from_constraint(m.x, m.c, eps=1e-6, linesearch=False) self.assertRegex( LOG.getvalue().strip(), - r"Setting Var 'x' to a numeric value `-?(?:\d+(?:\.\d*)?|\.\d+)[eE][+-]?\d+` outside the " + r"Setting Var 'x' to a numeric value `" + + SCIENTIFIC_NOTATION_REGEX + + "` outside the " r"bounds \(0, None\).", ) self.assertAlmostEqual(value(m.x), 0) From efb1fb3951a67345f371665de9cc03589a12b6b4 Mon Sep 17 00:00:00 2001 From: Doug A Date: Mon, 13 Jul 2026 11:25:59 -0400 Subject: [PATCH 3/6] fix line search decrease condition --- pyomo/util/calc_var_value.py | 4 ++-- pyomo/util/tests/test_calc_var_value.py | 4 ---- 2 files changed, 2 insertions(+), 6 deletions(-) diff --git a/pyomo/util/calc_var_value.py b/pyomo/util/calc_var_value.py index b6855307c58..738a71d0779 100644 --- a/pyomo/util/calc_var_value.py +++ b/pyomo/util/calc_var_value.py @@ -363,10 +363,10 @@ def calculate_variable_from_constraint( if fkp1.__class__ in _invalid_types: # We cannot perform computations on complex numbers fkp1 = None + # Why isn't this the Armijo condition? if ( fkp1 is not None - and (fkp1 * sf_constraint / sf_variable) ** 2 - < c1 * (fk * sf_constraint / sf_variable) ** 2 + and (fkp1 * sf_constraint) ** 2 < c1 * (fk * sf_constraint) ** 2 ): # found an alpha value with sufficient reduction # continue to the next step diff --git a/pyomo/util/tests/test_calc_var_value.py b/pyomo/util/tests/test_calc_var_value.py index 5d3dbea8406..1bf0ce63a7b 100644 --- a/pyomo/util/tests/test_calc_var_value.py +++ b/pyomo/util/tests/test_calc_var_value.py @@ -573,7 +573,3 @@ def test_linea_search_overflow(self): calculate_variable_from_constraint(m.x, m.con) self.assertAlmostEqual(value(m.x), 0) - - -if __name__ == "__main__": - unittest.main() From 5724e8a69a764c45b61341fd5138cec33c916a5d Mon Sep 17 00:00:00 2001 From: Doug A Date: Tue, 14 Jul 2026 10:44:03 -0400 Subject: [PATCH 4/6] copy_scaling_factors --- pyomo/util/subsystems.py | 56 ++++++++++- pyomo/util/tests/test_subsystems.py | 145 ++++++++++++++++++++++++++++ 2 files changed, 200 insertions(+), 1 deletion(-) diff --git a/pyomo/util/subsystems.py b/pyomo/util/subsystems.py index 49111858be9..8fb3271ad35 100644 --- a/pyomo/util/subsystems.py +++ b/pyomo/util/subsystems.py @@ -7,16 +7,18 @@ # software. This software is distributed under the 3-clause BSD License. # ____________________________________________________________________________________ -from pyomo.core.base.block import Block +from pyomo.core.base.block import Block, BlockData from pyomo.core.base.reference import Reference from pyomo.core.expr.visitor import identify_variables from pyomo.common.collections import ComponentSet, ComponentMap from pyomo.common.modeling import unique_component_name from pyomo.util.vars_from_expressions import get_vars_from_components +from pyomo.core.base.var import Var from pyomo.core.base.constraint import Constraint from pyomo.core.base.expression import Expression from pyomo.core.base.objective import Objective from pyomo.core.base.external import ExternalFunction +from pyomo.core.base.suffix import Suffix, SuffixFinder from pyomo.core.expr.visitor import StreamBasedExpressionVisitor from pyomo.core.expr.numeric_expr import ExternalFunctionExpression from pyomo.core.expr.numvalue import native_types, NumericValue @@ -392,3 +394,55 @@ def __next__(self): ) return inputs, outputs + + +# Based on code that John Eslick originally wrote for IDAES. +################################################################################# +# The Institute for the Design of Advanced Energy Systems Integrated Platform +# Framework (IDAES IP) was produced under the DOE Institute for the +# Design of Advanced Energy Systems (IDAES). +# +# Copyright (c) 2018-2026 by the software owners: The Regents of the +# University of California, through Lawrence Berkeley National Laboratory, +# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon +# University, West Virginia University Research Corporation, et al. +# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md +# for full copyright and license information. +################################################################################# + + +def copy_scaling_factors_to_block( + source_block: BlockData, destination_block: BlockData, overwrite=True +): + """ + Copies variable and constraint scaling factors from one block heirarchy to + another. Copying scaling factors is useful when creating a temporary block + to solve a subproblem because the .nl writer will not use any scaling factors + on the original model (see Pyomo #3800). + + Arguments + --------- + source_block: BlockData + Scaling factors will be copied from this BlockData. Any scaling factors + from any child blocks will also be copied, but *not* from any parent block(s) + destination_block: BlockData + Scaling factors will be copied to the scaling_factor Suffix on this BlockData. + If no such suffix exists, it will be created. + overwrite: Bool + Whether to overwrite any scaling factors that may already exist on destination_block. + + Returns + ------- + None + + """ + if not hasattr(destination_block, "scaling_factor"): + destination_block.scaling_factor = Suffix(direction=Suffix.EXPORT) + finder = SuffixFinder('scaling_factor', None, source_block) + for comp in destination_block.component_data_objects((Var, Constraint)): + sf = finder.find(comp) + if sf is not None: + if not overwrite and comp in destination_block.scaling_factor: + # Component already has a scaling factor, so don't overwrite it. + continue + destination_block.scaling_factor[comp] = sf diff --git a/pyomo/util/tests/test_subsystems.py b/pyomo/util/tests/test_subsystems.py index 327b05d3e6c..59c12f3ed1d 100644 --- a/pyomo/util/tests/test_subsystems.py +++ b/pyomo/util/tests/test_subsystems.py @@ -17,6 +17,7 @@ ParamSweeper, identify_external_functions, add_local_external_functions, + copy_scaling_factors_to_block, ) from pyomo.common.gsl import find_GSL @@ -754,5 +755,149 @@ def test_with_exception(self): self.assertFalse(m.v4.fixed) +# Code originally based on example in IDAES. +################################################################################# +# The Institute for the Design of Advanced Energy Systems Integrated Platform +# Framework (IDAES IP) was produced under the DOE Institute for the +# Design of Advanced Energy Systems (IDAES). +# +# Copyright (c) 2018-2026 by the software owners: The Regents of the +# University of California, through Lawrence Berkeley National Laboratory, +# National Technology & Engineering Solutions of Sandia, LLC, Carnegie Mellon +# University, West Virginia University Research Corporation, et al. +# All rights reserved. Please see the files COPYRIGHT.md and LICENSE.md +# for full copyright and license information. +################################################################################# + + +def _create_model(): + m = pyo.ConcreteModel() + m.s = pyo.Set(initialize=[1, 2, 3, 4]) + + m.v = pyo.Var(m.s) + + @m.Constraint(m.s) + def c(b, i): + return b.v[i] == i + + m.b = pyo.Block(m.s) + + for bd in m.b.values(): + bd.v2 = pyo.Var() + + return m + + +def _scale_model(m): + for bd in m.b.values(): + bd.scaling_factor = pyo.Suffix(direction=pyo.Suffix.EXPORT) + bd.scaling_factor[bd.v2] = 10 + + m.scaling_factor = pyo.Suffix(direction=pyo.Suffix.EXPORT) + for i in m.s: + m.scaling_factor[m.v[i]] = 5 * i + m.scaling_factor[m.c[i]] = 5**i + + return m + + +class TestCopyScalingFactorsToBlock(unittest.TestCase): + def test_copy_scaling_factors_to_block(self): + m = _create_model() + _scale_model(m) + + blk = create_subsystem_block( + constraints=[m.c[1], m.c[3]], + variables=[m.v[1], m.b[2].v2], # m.v[3] should be an input variable + ) + + assert not hasattr(blk, "scaling_factor") + copy_scaling_factors_to_block(m, blk) + assert hasattr(blk, "scaling_factor") + assert len(blk.scaling_factor) == 5 + assert blk.scaling_factor[m.v[1]] == 5 + assert blk.scaling_factor[m.v[3]] == 15 + assert blk.scaling_factor[m.b[2].v2] == 10 + assert blk.scaling_factor[m.c[1]] == 5 + assert blk.scaling_factor[m.c[3]] == 125 + + def test_copy_scaling_factors_preexisting_scaling_factor(self): + # Now test what happens if a scaling factor suffix already exits + m = _create_model() + _scale_model(m) + + blk = create_subsystem_block( + constraints=[m.c[1], m.c[3]], + variables=[m.v[1], m.b[2].v2], # m.v[3] should be an input variable + ) + blk.scaling_factor = pyo.Suffix(direction=pyo.Suffix.EXPORT) + blk.scaling_factor[m.b[2].v2] = 137 + # If a scaling factor doesn't exist on the original model + # for a variable, it should leave the existing scaling factor + # on the temporary block + del m.b[2].scaling_factor[m.b[2].v2] + copy_scaling_factors_to_block(m, blk) + assert hasattr(blk, "scaling_factor") + assert len(blk.scaling_factor) == 5 + assert blk.scaling_factor[m.v[1]] == 5 + assert blk.scaling_factor[m.v[3]] == 15 + assert blk.scaling_factor[m.b[2].v2] == 137 + assert blk.scaling_factor[m.c[1]] == 5 + assert blk.scaling_factor[m.c[3]] == 125 + + def test_copy_scaling_factors_overwrite_False(self): + m = _create_model() + _scale_model(m) + + blk = create_subsystem_block( + constraints=[m.c[1], m.c[3]], + variables=[m.v[1], m.b[2].v2], # m.v[3] should be an input variable + ) + blk.scaling_factor = pyo.Suffix(direction=pyo.Suffix.EXPORT) + blk.scaling_factor[m.b[2].v2] = 137 + + copy_scaling_factors_to_block(m, blk, overwrite=False) + assert hasattr(blk, "scaling_factor") + assert len(blk.scaling_factor) == 5 + assert blk.scaling_factor[m.v[1]] == 5 + assert blk.scaling_factor[m.v[3]] == 15 + assert blk.scaling_factor[m.b[2].v2] == 137 + assert blk.scaling_factor[m.c[1]] == 5 + assert blk.scaling_factor[m.c[3]] == 125 + + def test_copy_scaling_factors_overwrite_true(self): + m = _create_model() + _scale_model(m) + + blk = create_subsystem_block( + constraints=[m.c[1], m.c[3]], + variables=[m.v[1], m.b[2].v2], # m.v[3] should be an input variable + ) + blk.scaling_factor = pyo.Suffix(direction=pyo.Suffix.EXPORT) + blk.scaling_factor[m.b[2].v2] = 137 + + copy_scaling_factors_to_block(m, blk, overwrite=True) + assert hasattr(blk, "scaling_factor") + assert len(blk.scaling_factor) == 5 + assert blk.scaling_factor[m.v[1]] == 5 + assert blk.scaling_factor[m.v[3]] == 15 + assert blk.scaling_factor[m.b[2].v2] == 10 + assert blk.scaling_factor[m.c[1]] == 5 + assert blk.scaling_factor[m.c[3]] == 125 + + def test_copy_scaling_factors_to_block_from_subblock(self): + m = _create_model() + _scale_model(m) + + blk = create_subsystem_block( + constraints=[m.c[1], m.c[3]], + variables=[m.v[1], m.b[2].v2], # m.v[3] should be an input variable + ) + # Only the scaling suffix on m.b[2] should be copied. + copy_scaling_factors_to_block(m.b[2], blk) + assert len(blk.scaling_factor) == 1 + assert blk.scaling_factor[m.b[2].v2] == 10 + + if __name__ == '__main__': unittest.main() From 56798a6af4a8644e483973d98c4e1d0347d965cc Mon Sep 17 00:00:00 2001 From: Doug A Date: Tue, 14 Jul 2026 15:26:26 -0400 Subject: [PATCH 5/6] copy scaling factors to temp blocks --- .../contrib/incidence_analysis/scc_solver.py | 21 +++- .../tests/test_scc_solver.py | 106 +++++++++++++++++- 2 files changed, 121 insertions(+), 6 deletions(-) diff --git a/pyomo/contrib/incidence_analysis/scc_solver.py b/pyomo/contrib/incidence_analysis/scc_solver.py index 9971f8c8bc7..e0aef22a9b6 100644 --- a/pyomo/contrib/incidence_analysis/scc_solver.py +++ b/pyomo/contrib/incidence_analysis/scc_solver.py @@ -11,7 +11,11 @@ from pyomo.core.base.constraint import Constraint from pyomo.util.calc_var_value import calculate_variable_from_constraint -from pyomo.util.subsystems import TemporarySubsystemManager, generate_subsystem_blocks +from pyomo.util.subsystems import ( + generate_subsystem_blocks, + TemporarySubsystemManager, + copy_scaling_factors_to_block, +) from pyomo.contrib.incidence_analysis.interface import ( IncidenceGraphInterface, _generate_variables_in_constraints, @@ -89,7 +93,13 @@ def generate_strongly_connected_components( def solve_strongly_connected_components( - block, *, solver=None, solve_kwds=None, use_calc_var=True, calc_var_kwds=None + block, + *, + solver=None, + solve_kwds=None, + use_calc_var=True, + calc_var_kwds=None, + copy_scaling_factors=False, ): """Solve a square system of variables and equality constraints by solving strongly connected components individually. @@ -119,7 +129,8 @@ def solve_strongly_connected_components( square system solves calc_var_kwds: Dictionary Keyword arguments for calculate_variable_from_constraint - + copy_scaling_factors: Bool + Whether to copy scaling factors to the temporary blocks before solving them. Returns ------- List of results objects returned by each call to solve @@ -129,6 +140,8 @@ def solve_strongly_connected_components( solve_kwds = {} if calc_var_kwds is None: calc_var_kwds = {} + if copy_scaling_factors and "scale_problem" not in calc_var_kwds: + calc_var_kwds["scale_problem"] = True igraph = IncidenceGraphInterface( block, @@ -154,6 +167,8 @@ def solve_strongly_connected_components( scc.vars[0], scc.cons[0], **calc_var_kwds ) else: + if copy_scaling_factors: + copy_scaling_factors_to_block(block.model(), scc) if solver is None: var_names = [var.name for var in scc.vars.values()][:10] con_names = [con.name for con in scc.cons.values()][:10] diff --git a/pyomo/contrib/incidence_analysis/tests/test_scc_solver.py b/pyomo/contrib/incidence_analysis/tests/test_scc_solver.py index 52c3da8eaed..bb9afdfd35c 100644 --- a/pyomo/contrib/incidence_analysis/tests/test_scc_solver.py +++ b/pyomo/contrib/incidence_analysis/tests/test_scc_solver.py @@ -12,6 +12,8 @@ from pyomo.common.dependencies import networkx_available from pyomo.common.dependencies import scipy_available from pyomo.common.collections import ComponentSet, ComponentMap +from pyomo.util.calc_var_value import calculate_variable_from_constraint +from pyomo.util.subsystems import copy_scaling_factors_to_block from pyomo.contrib.incidence_analysis.scc_solver import ( TemporarySubsystemManager, generate_strongly_connected_components, @@ -21,9 +23,26 @@ make_gas_expansion_model, make_dynamic_model, ) +from pyomo.contrib.solver.solvers.ipopt import Ipopt import pyomo.common.unittest as unittest +def _make_example_model(): + m = pyo.ConcreteModel() + m.w = pyo.Var(initialize=1) + m.x = pyo.Var(initialize=2) + m.y = pyo.Var(initialize=3) + m.z = pyo.Var(initialize=4) + + m.eqn1 = pyo.Constraint(expr=m.w**2 == 4) + m.eqn2 = pyo.Constraint(expr=m.w + m.x + m.y == 3) + m.eqn3 = pyo.Constraint(expr=m.w + 3 * m.x - m.y == 5) + m.eqn4 = pyo.Constraint(expr=m.z == m.y**2 + m.y - 2) + + # This model has a solution w=2, x=1, y=0, z=-2 + return m + + @unittest.skipUnless(scipy_available, "SciPy is not available") @unittest.skipUnless(networkx_available, "NetworkX is not available") class TestGenerateSCC(unittest.TestCase): @@ -194,7 +213,7 @@ def test_dynamic_backward_disc_without_initial_conditions(self): with TemporarySubsystemManager(to_fix=inputs): # We have a much easier time testing the SCCs generated # in this test. - t = m.time[i + 2] + t = m.time.at(i + 2) t_prev = m.time.prev(t) con_set = ComponentSet( @@ -261,7 +280,7 @@ def test_dynamic_backward_with_inputs(self): generate_strongly_connected_components(constraints, variables) ): with TemporarySubsystemManager(to_fix=inputs): - t = m.time[i + 2] + t = m.time.at(i + 2) t_prev = m.time.prev(t) con_set = ComponentSet( @@ -324,7 +343,7 @@ def test_dynamic_forward_disc(self): # algebraic -> derivative -> differential -> algebraic -> ... idx = i // 3 mod = i % 3 - t = m.time[idx + 1] + t = m.time.at(idx + 1) if t != time.last(): t_next = m.time.next(t) @@ -498,6 +517,87 @@ def test_with_inequalities(self): self.assertAlmostEqual(m.x[2].value, 3.21642835) self.assertEqual(m.x[3].value, 1.0) + def _meta_copy_scaling_factors(self, copy_scaling_factors, scale_problem): + """ + Parameterized test function to ensure that submodel settings are + getting set as expected. + """ + if copy_scaling_factors is None: + expected_copy_scaling_factors = False + else: + expected_copy_scaling_factors = copy_scaling_factors + + if scale_problem is None: + expected_scale_problem = expected_copy_scaling_factors + else: + expected_scale_problem = scale_problem + + m = _make_example_model() + m.scaling_factor = pyo.Suffix(direction=pyo.Suffix.EXPORT) + m.scaling_factor[m.x] = 5 + m.scaling_factor[m.eqn3] = 47 + + mock_calculate_variable_from_constraint = unittest.mock.MagicMock( + spec=calculate_variable_from_constraint + ) + + def wrap_cvfc(variable, constraint, **kwargs): + assert kwargs["scale_problem"] == expected_scale_problem + return calculate_variable_from_constraint(variable, constraint, kwargs) + + mock_calculate_variable_from_constraint + + mock_IPOPT = unittest.mock.MagicMock(Ipopt, autospec=True) + + def solve_two_by_two(blk, **kwargs): + if expected_copy_scaling_factors is True: + assert len(blk.scaling_factor) == 2 + assert blk.scaling_factor[m.x] == 5 + assert blk.scaling_factor[m.eqn3] == 47 + else: + assert not hasattr(blk, "scaling_factor") + m.x.set_value(1) + m.y.set_value(0) + # Creating a Results object or a mocked version + # thereof doesn't provide much benefit here. + return "Solve Successful" + + mock_IPOPT.solve.side_effect = solve_two_by_two + + if scale_problem is None: + calc_var_kwds = {} + else: + calc_var_kwds = {"scale_problem": scale_problem} + + with ( + unittest.mock.patch( + "pyomo.contrib.incidence_analysis.scc_solver.calculate_variable_from_constraint", + new=mock_calculate_variable_from_constraint, + ) as mock_calc_var_from_con, + unittest.mock.patch( + "pyomo.contrib.incidence_analysis.scc_solver.copy_scaling_factors_to_block", + spec=copy_scaling_factors_to_block, + wraps=copy_scaling_factors_to_block, + ) as mock_copy_sf, + ): + solve_strongly_connected_components( + m, + solver=mock_IPOPT, + calc_var_kwds=calc_var_kwds, + copy_scaling_factors=copy_scaling_factors, + ) + assert mock_calc_var_from_con.call_count == 2 + if expected_copy_scaling_factors is True: + assert mock_copy_sf.call_count == 1 + + def test_copy_scaling_factors(self): + for copy_scaling_factors in [True, False, None]: + for scale_problem in [True, False, None]: + self._meta_copy_scaling_factors( + copy_scaling_factors=copy_scaling_factors, + scale_problem=scale_problem, + ) + @unittest.skipUnless(scipy_available, "SciPy is not available") @unittest.skipUnless(networkx_available, "NetworkX is not available") From f10283dc816c564decb1719d538274ba412e46e6 Mon Sep 17 00:00:00 2001 From: Doug A Date: Mon, 20 Jul 2026 12:11:01 -0400 Subject: [PATCH 6/6] typo --- pyomo/util/subsystems.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyomo/util/subsystems.py b/pyomo/util/subsystems.py index 8fb3271ad35..80db72d54c4 100644 --- a/pyomo/util/subsystems.py +++ b/pyomo/util/subsystems.py @@ -415,7 +415,7 @@ def copy_scaling_factors_to_block( source_block: BlockData, destination_block: BlockData, overwrite=True ): """ - Copies variable and constraint scaling factors from one block heirarchy to + Copies variable and constraint scaling factors from one block hierarchy to another. Copying scaling factors is useful when creating a temporary block to solve a subproblem because the .nl writer will not use any scaling factors on the original model (see Pyomo #3800).