Clip and scale calculate_variable_from_constraint

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)

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
@@ -36,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)
@@ -240,7 +245,80 @@ 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, 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.
         # This tests that the algorithm cleanly continues
@@ -413,6 +491,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