Use scaling factors in solve_strongly_connected_components

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]

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")

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)