From 0dd64c75155fa5b1761c5dfa4cb1e20ea1e70943 Mon Sep 17 00:00:00 2001 From: Doug A Date: Fri, 10 Jul 2026 14:59:37 -0400 Subject: [PATCH 1/5] 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/5] 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/5] 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 dbf30f8c0b07e8433faedf3c7d9a742443f5cd80 Mon Sep 17 00:00:00 2001 From: Doug A Date: Tue, 21 Jul 2026 10:22:20 -0400 Subject: [PATCH 4/5] exception not raised on some platforms --- pyomo/util/tests/test_calc_var_value.py | 23 ++++++++++++++--------- 1 file changed, 14 insertions(+), 9 deletions(-) diff --git a/pyomo/util/tests/test_calc_var_value.py b/pyomo/util/tests/test_calc_var_value.py index 1bf0ce63a7b..0d8908746d1 100644 --- a/pyomo/util/tests/test_calc_var_value.py +++ b/pyomo/util/tests/test_calc_var_value.py @@ -8,6 +8,7 @@ # ____________________________________________________________________________________ import logging +import re from io import StringIO import pyomo.common.unittest as unittest @@ -300,23 +301,27 @@ def test_nonlinear(self): 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 + # a zero derivative. It still raises an error on some platforms, 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"\.", - ): + + try: calculate_variable_from_constraint( m.x, m.g, diff_mode=mode, scale_problem=True ) + except IterationLimitError as err: + out = re.match( + "Linesearch iteration limit reached solving for variable 'x' using " + "constraint 'g'; remaining residual = " + + SCIENTIFIC_NOTATION_REGEX + + r"\.", + err.message, + ) + assert out is not None # Calculate the bubble point of Benzene. The first step # computed by calculate_variable_from_constraint will make the From 922c74d305f003b1ecd4a766d6a23e41c4952ad9 Mon Sep 17 00:00:00 2001 From: Doug A Date: Wed, 22 Jul 2026 13:50:12 -0400 Subject: [PATCH 5/5] message was deprecated --- pyomo/util/tests/test_calc_var_value.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyomo/util/tests/test_calc_var_value.py b/pyomo/util/tests/test_calc_var_value.py index 0d8908746d1..cab2e8ce5df 100644 --- a/pyomo/util/tests/test_calc_var_value.py +++ b/pyomo/util/tests/test_calc_var_value.py @@ -319,7 +319,7 @@ def test_nonlinear(self): "constraint 'g'; remaining residual = " + SCIENTIFIC_NOTATION_REGEX + r"\.", - err.message, + str(err), ) assert out is not None