Add diagnostics tests for all unit models (#1375)

* 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

---------

Co-authored-by: Robert Parker <robbybparker@gmail.com>
Co-authored-by: Keith Beattie <ksbeattie@lbl.gov>

idaes/__init__.py

@@ -22,6 +22,9 @@
 import os
 import copy
 import logging
+from typing import Optional, List
+
+from pyomo.common.fileutils import find_library
 
 from . import config
 from .ver import __version__  # noqa
@@ -83,6 +86,29 @@ def _handle_optional_compat_activation(
 _log.debug("'idaes' logger debug test")
 
 
+# TODO: Remove once AMPL bug is fixed
+# TODO: https://github.com/ampl/asl/issues/13
+# There appears to be a bug in the ASL which causes terminal failures
+# if you try to create multiple ASL structs with different external
+# functions in the same process. This causes pytest to crash during testing.
+# To avoid this, register all known external functions at initialization.
+def _ensure_external_functions_libs_in_env(
+    ext_funcs: List[str], var_name: str = "AMPLFUNC", sep: str = "\n"
+):
+    libraries_str = os.environ.get(var_name, "")
+    libraries = [lib for lib in libraries_str.split(sep) if lib.strip()]
+    for func_name in ext_funcs:
+        lib: Optional[str] = find_library(os.path.join(bin_directory, func_name))
+        if lib is not None and lib not in libraries:
+            libraries.append(lib)
+    os.environ[var_name] = sep.join(libraries)
+
+
+_ensure_external_functions_libs_in_env(
+    ["cubic_roots", "general_helmholtz_external", "functions"]
+)
+
+
 def _create_data_dir():
     """Create the IDAES directory to store data files in."""
     config.create_dir(data_directory)

idaes/core/util/model_diagnostics.py

@@ -997,16 +997,25 @@ def display_near_parallel_variables(self, stream=None):
     # TODO: Block triangularization analysis
     # Number and size of blocks, polynomial degree of 1x1 blocks, simple pivot check of moderate sized sub-blocks?
 
-    def _collect_structural_warnings(self):
+    def _collect_structural_warnings(
+        self, ignore_evaluation_errors=False, ignore_unit_consistency=False
+    ):
         """
         Runs checks for structural warnings and returns two lists.
 
+        Args:
+            ignore_evaluation_errors - ignore potential evaluation error warnings
+            ignore_unit_consistency - ignore unit consistency warnings
+
         Returns:
             warnings - list of warning messages from structural analysis
             next_steps - list of suggested next steps to further investigate warnings
 
         """
-        uc = identify_inconsistent_units(self._model)
+        if not ignore_unit_consistency:
+            uc = identify_inconsistent_units(self._model)
+        else:
+            uc = []
         uc_var, uc_con, oc_var, oc_con = self.get_dulmage_mendelsohn_partition()
 
         # Collect warnings
@@ -1040,12 +1049,15 @@ def _collect_structural_warnings(self):
         if any(len(x) > 0 for x in [oc_var, oc_con]):
             next_steps.append(self.display_overconstrained_set.__name__ + "()")
 
-        eval_warnings = self._collect_potential_eval_errors()
-        if len(eval_warnings) > 0:
-            warnings.append(
-                f"WARNING: Found {len(eval_warnings)} potential evaluation errors."
-            )
-            next_steps.append(self.display_potential_evaluation_errors.__name__ + "()")
+        if not ignore_evaluation_errors:
+            eval_warnings = self._collect_potential_eval_errors()
+            if len(eval_warnings) > 0:
+                warnings.append(
+                    f"WARNING: Found {len(eval_warnings)} potential evaluation errors."
+                )
+                next_steps.append(
+                    self.display_potential_evaluation_errors.__name__ + "()"
+                )
 
         return warnings, next_steps
 
@@ -1289,16 +1301,27 @@ def _collect_numerical_cautions(self, jac=None, nlp=None):
 
         return cautions
 
-    def assert_no_structural_warnings(self):
+    def assert_no_structural_warnings(
+        self,
+        ignore_evaluation_errors: bool = False,
+        ignore_unit_consistency: bool = False,
+    ):
         """
         Checks for structural warnings in the model and raises an AssertionError
         if any are found.
 
+        Args:
+            ignore_evaluation_errors - ignore potential evaluation error warnings
+            ignore_unit_consistency - ignore unit consistency warnings
+
         Raises:
             AssertionError if any warnings are identified by structural analysis.
 
         """
-        warnings, _ = self._collect_structural_warnings()
+        warnings, _ = self._collect_structural_warnings(
+            ignore_evaluation_errors=ignore_evaluation_errors,
+            ignore_unit_consistency=ignore_unit_consistency,
+        )
         if len(warnings) > 0:
             raise AssertionError(f"Structural issues found ({len(warnings)}).")
 
@@ -1387,9 +1410,17 @@ def report_numerical_issues(self, stream=None):
         cautions = self._collect_numerical_cautions(jac=jac, nlp=nlp)
 
         stats = []
-        stats.append(
-            f"Jacobian Condition Number: {jacobian_cond(jac=jac, scaled=False):.3E}"
-        )
+        try:
+            stats.append(
+                f"Jacobian Condition Number: {jacobian_cond(jac=jac, scaled=False):.3E}"
+            )
+        except RuntimeError as err:
+            if "Factor is exactly singular" in str(err):
+                _log.info(err)
+                stats.append("Jacobian Condition Number: Undefined (Exactly Singular)")
+            else:
+                raise
+
         _write_report_section(
             stream=stream, lines_list=stats, title="Model Statistics", header="="
         )

idaes/core/util/tests/test_model_diagnostics.py

@@ -33,6 +33,7 @@
     acos,
     sqrt,
     Objective,
+    PositiveIntegers,
     Set,
     SolverFactory,
     Suffix,
@@ -1342,6 +1343,46 @@ def test_report_numerical_issues_ok(self):
         prepare_degeneracy_hunter()
         prepare_svd_toolbox()
 
+====================================================================================
+"""
+
+        assert stream.getvalue() == expected
+
+    @pytest.mark.component
+    def test_report_numerical_issues_exactly_singular(self):
+        m = ConcreteModel()
+        m.x = Var([1, 2], initialize=1.0)
+        m.eq = Constraint(PositiveIntegers)
+        m.eq[1] = m.x[1] * m.x[2] == 1.5
+        m.eq[2] = m.x[2] * m.x[1] == 1.5
+        m.obj = Objective(expr=m.x[1] ** 2 + 2 * m.x[2] ** 2)
+
+        dt = DiagnosticsToolbox(m)
+        dt.report_numerical_issues()
+
+        stream = StringIO()
+        dt.report_numerical_issues(stream)
+
+        expected = """====================================================================================
+Model Statistics
+
+    Jacobian Condition Number: Undefined (Exactly Singular)
+
+------------------------------------------------------------------------------------
+1 WARNINGS
+
+    WARNING: 2 Constraints with large residuals (>1.0E-05)
+
+------------------------------------------------------------------------------------
+0 Cautions
+
+    No cautions found!
+
+------------------------------------------------------------------------------------
+Suggested next steps:
+
+    display_constraints_with_large_residuals()
+
 ====================================================================================
 """
 
@@ -2241,7 +2282,7 @@ def test_run_ipopt_with_stats(self):
         assert iters == 1
         assert iters_in_restoration == 0
         assert iters_w_regularization == 0
-        assert time < 0.01
+        assert isinstance(time, float)
 
     @pytest.mark.component
     @pytest.mark.solver

idaes/models/properties/examples/saponification_thermo.py

@@ -83,15 +83,15 @@ def build(self):
 
         # Heat capacity of water
         self.cp_mol = Param(
-            mutable=False,
+            mutable=True,
             initialize=75.327,
             doc="Molar heat capacity of water [J/mol.K]",
             units=units.J / units.mol / units.K,
         )
 
         # Density of water
         self.dens_mol = Param(
-            mutable=False,
+            mutable=True,
             initialize=55388.0,
             doc="Molar density of water [mol/m^3]",
             units=units.mol / units.m**3,

idaes/models/unit_models/heat_exchanger.py

@@ -814,11 +814,11 @@ def _get_performance_contents(self, time_point=0):
         var_dict["Heat Duty"] = self.heat_duty[time_point]
         if self.config.flow_pattern == HeatExchangerFlowPattern.crossflow:
             var_dict["Crossflow Factor"] = self.crossflow_factor[time_point]
+        var_dict["Delta T In"] = self.delta_temperature_in[time_point]
+        var_dict["Delta T Out"] = self.delta_temperature_out[time_point]
 
         expr_dict = {}
         expr_dict["Delta T Driving"] = self.delta_temperature[time_point]
-        expr_dict["Delta T In"] = self.delta_temperature_in[time_point]
-        expr_dict["Delta T Out"] = self.delta_temperature_out[time_point]
 
         return {"vars": var_dict, "exprs": expr_dict}
 

idaes/models/unit_models/tests/test_cstr.py

@@ -18,7 +18,6 @@
 import pytest
 
 from pyomo.environ import check_optimal_termination, ConcreteModel, units, value
-from pyomo.util.check_units import assert_units_consistent, assert_units_equivalent
 
 from idaes.core import (
     FlowsheetBlock,
@@ -34,7 +33,6 @@
     SaponificationReactionParameterBlock,
 )
 from idaes.core.util.model_statistics import (
-    degrees_of_freedom,
     number_variables,
     number_total_constraints,
     number_unused_variables,
@@ -50,6 +48,8 @@
     SingleControlVolumeUnitInitializer,
     InitializationStatus,
 )
+from idaes.core.util import DiagnosticsToolbox
+
 
 # -----------------------------------------------------------------------------
 # Get default solver for testing
@@ -109,8 +109,8 @@ def sapon(self):
         m.fs.unit.inlet.conc_mol_comp[0, "H2O"].fix(55388.0)
         m.fs.unit.inlet.conc_mol_comp[0, "NaOH"].fix(100.0)
         m.fs.unit.inlet.conc_mol_comp[0, "EthylAcetate"].fix(100.0)
-        m.fs.unit.inlet.conc_mol_comp[0, "SodiumAcetate"].fix(0.0)
-        m.fs.unit.inlet.conc_mol_comp[0, "Ethanol"].fix(0.0)
+        m.fs.unit.inlet.conc_mol_comp[0, "SodiumAcetate"].fix(1e-8)
+        m.fs.unit.inlet.conc_mol_comp[0, "Ethanol"].fix(1e-8)
 
         m.fs.unit.inlet.temperature.fix(303.15)
         m.fs.unit.inlet.pressure.fix(101325.0)
@@ -149,15 +149,9 @@ def test_build(self, sapon):
         assert number_unused_variables(sapon) == 0
 
     @pytest.mark.component
-    def test_units(self, sapon):
-        assert_units_consistent(sapon)
-        assert_units_equivalent(sapon.fs.unit.volume[0], units.m**3)
-        assert_units_equivalent(sapon.fs.unit.heat_duty[0], units.W)
-        assert_units_equivalent(sapon.fs.unit.deltaP[0], units.Pa)
-
-    @pytest.mark.unit
-    def test_dof(self, sapon):
-        assert degrees_of_freedom(sapon) == 0
+    def test_structural_issues(self, sapon):
+        dt = DiagnosticsToolbox(sapon)
+        dt.assert_no_structural_warnings()
 
     @pytest.mark.solver
     @pytest.mark.skipif(solver is None, reason="Solver not available")
@@ -248,6 +242,13 @@ def test_conservation(self, sapon):
             <= 1e-3
         )
 
+    @pytest.mark.solver
+    @pytest.mark.skipif(solver is None, reason="Solver not available")
+    @pytest.mark.component
+    def test_numerical_issues(self, sapon):
+        dt = DiagnosticsToolbox(sapon)
+        dt.assert_no_numerical_warnings()
+
     @pytest.mark.ui
     @pytest.mark.unit
     def test_get_performance_contents(self, sapon):

idaes/models/unit_models/tests/test_equilibrium_reactor.py

@@ -33,7 +33,6 @@
     SaponificationReactionParameterBlock,
 )
 from idaes.core.util.model_statistics import (
-    degrees_of_freedom,
     number_variables,
     number_total_constraints,
     number_unused_variables,
@@ -44,14 +43,12 @@
     initialization_tester,
 )
 from idaes.core.solvers import get_solver
-from pyomo.util.check_units import assert_units_consistent, assert_units_equivalent
-from idaes.core.util.exceptions import InitializationError
-
 from idaes.core.initialization import (
     BlockTriangularizationInitializer,
     SingleControlVolumeUnitInitializer,
     InitializationStatus,
 )
+from idaes.core.util import DiagnosticsToolbox
 
 # -----------------------------------------------------------------------------
 # Get default solver for testing
@@ -116,8 +113,8 @@ def sapon(self):
         m.fs.unit.inlet.conc_mol_comp[0, "H2O"].fix(55388.0)
         m.fs.unit.inlet.conc_mol_comp[0, "NaOH"].fix(100.0)
         m.fs.unit.inlet.conc_mol_comp[0, "EthylAcetate"].fix(100.0)
-        m.fs.unit.inlet.conc_mol_comp[0, "SodiumAcetate"].fix(0.0)
-        m.fs.unit.inlet.conc_mol_comp[0, "Ethanol"].fix(0.0)
+        m.fs.unit.inlet.conc_mol_comp[0, "SodiumAcetate"].fix(1e-8)
+        m.fs.unit.inlet.conc_mol_comp[0, "Ethanol"].fix(1e-8)
 
         m.fs.unit.inlet.temperature.fix(303.15)
         m.fs.unit.inlet.pressure.fix(101325.0)
@@ -153,14 +150,9 @@ def test_build(self, sapon):
         assert number_unused_variables(sapon) == 0
 
     @pytest.mark.component
-    def test_units(self, sapon):
-        assert_units_consistent(sapon)
-        assert_units_equivalent(sapon.fs.unit.heat_duty[0], units.W)
-        assert_units_equivalent(sapon.fs.unit.deltaP[0], units.Pa)
-
-    @pytest.mark.unit
-    def test_dof(self, sapon):
-        assert degrees_of_freedom(sapon) == 0
+    def test_structural_issues(self, sapon):
+        dt = DiagnosticsToolbox(sapon)
+        dt.assert_no_structural_warnings()
 
     @pytest.mark.solver
     @pytest.mark.skipif(solver is None, reason="Solver not available")
@@ -241,6 +233,13 @@ def test_conservation(self, sapon):
             <= 1e-1
         )
 
+    @pytest.mark.solver
+    @pytest.mark.skipif(solver is None, reason="Solver not available")
+    @pytest.mark.component
+    def test_numerical_issues(self, sapon):
+        dt = DiagnosticsToolbox(sapon)
+        dt.assert_no_numerical_warnings()
+
     @pytest.mark.ui
     @pytest.mark.unit
     def test_get_performance_contents(self, sapon):