Add parallel constraint/variable check to report_numerical_issues (#…

…1385)

* update imports of native_types and pyomo_constant_types (which is deprecated)

* Adding next batch of diagnostics tests

* Next batch of daignsotics tests

* Work around for ASL issue

* Adding more diagnostics checks

* Last unit model diagnostics tests

* Fixing typo

* Improving fix for ASL issue

* Better implementation of fix

* Fixing pylint and Python 3.8 failures

* Removing old implementation of workaround

* Fixing noisy test

* Moving registration of external functions

* Reverting to version that works on Windows

* Trying another way to get binary files

* add parallel variable/constraint checks to report_numerical_issues; change default parallel tolerance to 1e-8

* update tests

* tighten parallel_component_tolerance to 1e-8

* adjust model to make parallel variable test less sensitive

* update test to reflect new default tolerance of 1e-8

* consistent format for displaying parallel tolerance

---------

Co-authored-by: Andrew Lee <andrew.lee@netl.doe.gov>

idaes/core/util/model_diagnostics.py

@@ -286,7 +286,7 @@ def svd_sparse(jacobian, number_singular_values):
 CONFIG.declare(
     "parallel_component_tolerance",
     ConfigValue(
-        default=1e-4,
+        default=1e-8,
         domain=float,
         description="Tolerance for identifying near-parallel Jacobian rows/columns",
     ),
@@ -1179,6 +1179,29 @@ def _collect_numerical_warnings(self, jac=None, nlp=None):
                 self.display_constraints_with_extreme_jacobians.__name__ + "()"
             )
 
+        # Parallel variables and constraints
+        partol = self.config.parallel_component_tolerance
+        par_cons = check_parallel_jacobian(
+            self._model, tolerance=partol, direction="row", jac=jac, nlp=nlp
+        )
+        par_vars = check_parallel_jacobian(
+            self._model, tolerance=partol, direction="column", jac=jac, nlp=nlp
+        )
+        if par_cons:
+            p = "pair" if len(par_cons) == 1 else "pairs"
+            warnings.append(
+                f"WARNING: {len(par_cons)} {p} of constraints are parallel"
+                f" (to tolerance {partol:.1E})"
+            )
+            next_steps.append(self.display_near_parallel_constraints.__name__ + "()")
+        if par_vars:
+            p = "pair" if len(par_vars) == 1 else "pairs"
+            warnings.append(
+                f"WARNING: {len(par_vars)} {p} of variables are parallel"
+                f" (to tolerance {partol:.1E})"
+            )
+            next_steps.append(self.display_near_parallel_variables.__name__ + "()")
+
         return warnings, next_steps
 
     def _collect_numerical_cautions(self, jac=None, nlp=None):
@@ -1441,7 +1464,6 @@ def report_numerical_issues(self, stream=None):
             lines_list=next_steps,
             title="Suggested next steps:",
             line_if_empty=f"If you still have issues converging your model consider:\n"
-            f"{TAB*2}display_near_parallel_constraints()\n{TAB*2}display_near_parallel_variables()"
             f"\n{TAB*2}prepare_degeneracy_hunter()\n{TAB*2}prepare_svd_toolbox()",
             footer="=",
         )
@@ -3619,14 +3641,25 @@ def ipopt_solve_halt_on_error(model, options=None):
     )
 
 
-def check_parallel_jacobian(model, tolerance: float = 1e-4, direction: str = "row"):
+def check_parallel_jacobian(
+    model,
+    tolerance: float = 1e-4,
+    direction: str = "row",
+    jac=None,
+    nlp=None,
+):
     """
     Check for near-parallel rows or columns in the Jacobian.
 
     Near-parallel rows or columns indicate a potential degeneracy in the model,
     as this means that the associated constraints or variables are (near)
     duplicates of each other.
 
+    For efficiency, the ``jac`` and ``nlp`` arguments may be provided if they are
+    already available. If these are provided, the provided model is not used. If
+    either ``jac`` or ``nlp`` is not provided, a Jacobian and ``PyomoNLP`` are
+    computed using the model.
+
     This method is based on work published in:
 
     Klotz, E., Identification, Assessment, and Correction of Ill-Conditioning and
@@ -3637,6 +3670,8 @@ def check_parallel_jacobian(model, tolerance: float = 1e-4, direction: str = "ro
         model: model to be analysed
         tolerance: tolerance to use to determine if constraints/variables are parallel
         direction: 'row' (default, constraints) or 'column' (variables)
+        jac: model Jacobian as a ``scipy.sparse.coo_matrix``, optional
+        nlp: ``PyomoNLP`` of model, optional
 
     Returns:
         list of 2-tuples containing parallel Pyomo components
@@ -3651,7 +3686,8 @@ def check_parallel_jacobian(model, tolerance: float = 1e-4, direction: str = "ro
             "Must be 'row' or 'column'."
         )
 
-    jac, nlp = get_jacobian(model, scaled=False)
+    if jac is None or nlp is None:
+        jac, nlp = get_jacobian(model, scaled=False)
 
     # Get vectors that we will check, and the Pyomo components
     # they correspond to.

idaes/core/util/tests/test_model_diagnostics.py

@@ -975,10 +975,10 @@ def test_display_near_parallel_variables(self):
         model.v3 = Var()
         model.v4 = Var()
 
-        model.c1 = Constraint(expr=model.v1 == model.v2 - 0.99999 * model.v4)
-        model.c2 = Constraint(expr=model.v1 + 1.00001 * model.v4 == 1e-8 * model.v3)
+        model.c1 = Constraint(expr=1e-8 * model.v1 == 1e-8 * model.v2 - 1e-8 * model.v4)
+        model.c2 = Constraint(expr=1e-8 * model.v1 + 1e-8 * model.v4 == model.v3)
         model.c3 = Constraint(
-            expr=1e8 * (model.v1 + model.v4) + 1e10 * model.v2 == 1e-6 * model.v3
+            expr=1e3 * (model.v1 + model.v4) + 1e3 * model.v2 == model.v3
         )
 
         dt = DiagnosticsToolbox(model=model)
@@ -1121,7 +1121,7 @@ def test_collect_numerical_warnings_jacobian(self):
 
         warnings, next_steps = dt._collect_numerical_warnings()
 
-        assert len(warnings) == 3
+        assert len(warnings) == 4
         assert (
             "WARNING: 2 Variables with extreme Jacobian values (<1.0E-08 or >1.0E+08)"
             in warnings
@@ -1132,7 +1132,7 @@ def test_collect_numerical_warnings_jacobian(self):
         )
         assert "WARNING: 1 Constraint with large residuals (>1.0E-05)" in warnings
 
-        assert len(next_steps) == 3
+        assert len(next_steps) == 4
         assert "display_variables_with_extreme_jacobians()" in next_steps
         assert "display_constraints_with_extreme_jacobians()" in next_steps
         assert "display_constraints_with_large_residuals()" in next_steps
@@ -1338,8 +1338,7 @@ def test_report_numerical_issues_ok(self):
 Suggested next steps:
 
     If you still have issues converging your model consider:
-        display_near_parallel_constraints()
-        display_near_parallel_variables()
+
         prepare_degeneracy_hunter()
         prepare_svd_toolbox()
 
@@ -1369,9 +1368,11 @@ def test_report_numerical_issues_exactly_singular(self):
     Jacobian Condition Number: Undefined (Exactly Singular)
 
 ------------------------------------------------------------------------------------
-1 WARNINGS
+3 WARNINGS
 
     WARNING: 2 Constraints with large residuals (>1.0E-05)
+    WARNING: 1 pair of constraints are parallel (to tolerance 1.0E-08)
+    WARNING: 1 pair of variables are parallel (to tolerance 1.0E-08)
 
 ------------------------------------------------------------------------------------
 0 Cautions
@@ -1382,6 +1383,8 @@ def test_report_numerical_issues_exactly_singular(self):
 Suggested next steps:
 
     display_constraints_with_large_residuals()
+    display_near_parallel_constraints()
+    display_near_parallel_variables()
 
 ====================================================================================
 """
@@ -1433,8 +1436,8 @@ def test_report_numerical_issues_jacobian(self):
         model.v2 = Var(initialize=0)
         model.v3 = Var(initialize=0)
 
-        model.c1 = Constraint(expr=model.v1 == model.v2)
-        model.c2 = Constraint(expr=model.v1 == 1e-8 * model.v3)
+        model.c1 = Constraint(expr=1e-2 * model.v1 == model.v2)
+        model.c2 = Constraint(expr=1e-2 * model.v1 == 1e-8 * model.v3)
         model.c3 = Constraint(expr=1e8 * model.v1 + 1e10 * model.v2 == 1e-6 * model.v3)
 
         dt = DiagnosticsToolbox(model=model)
@@ -1445,14 +1448,15 @@ def test_report_numerical_issues_jacobian(self):
         expected = """====================================================================================
 Model Statistics
 
-    Jacobian Condition Number: 1.407E+18
+    Jacobian Condition Number: 1.118E+18
 
 ------------------------------------------------------------------------------------
-3 WARNINGS
+4 WARNINGS
 
     WARNING: 1 Constraint with large residuals (>1.0E-05)
     WARNING: 2 Variables with extreme Jacobian values (<1.0E-08 or >1.0E+08)
     WARNING: 1 Constraint with extreme Jacobian values (<1.0E-08 or >1.0E+08)
+    WARNING: 3 pairs of variables are parallel (to tolerance 1.0E-08)
 
 ------------------------------------------------------------------------------------
 4 Cautions
@@ -1468,6 +1472,7 @@ def test_report_numerical_issues_jacobian(self):
     display_constraints_with_large_residuals()
     display_variables_with_extreme_jacobians()
     display_constraints_with_extreme_jacobians()
+    display_near_parallel_variables()
 
 ====================================================================================
 """