diff --git a/pyomo/util/calc_var_value.py b/pyomo/util/calc_var_value.py index c9a58b68f6e..738a71d0779 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" @@ -231,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) @@ -258,7 +335,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( @@ -286,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: + # Why isn't this the Armijo condition? + if ( + fkp1 is not None + and (fkp1 * sf_constraint) ** 2 < c1 * (fk * sf_constraint) ** 2 + ): # found an alpha value with sufficient reduction # continue to the next step fk = fkp1 @@ -310,7 +391,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..cab2e8ce5df 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 @@ -24,6 +25,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 @@ -36,6 +38,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) @@ -240,7 +246,84 @@ 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) + + # 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 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) + + 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"\.", + str(err), + ) + assert out is not None + + # 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 +496,35 @@ 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): + + 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 `" + + SCIENTIFIC_NOTATION_REGEX + + "` 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