From 3e403b78e3aaf5d6c69b1995c05f859548d5a097 Mon Sep 17 00:00:00 2001
From: Andrew Lee <andrew.lee@netl.doe.gov>
Date: Tue, 5 Nov 2024 12:23:07 -0500
Subject: [PATCH 1/2] Diagnostics tool for ill-posed constraints (#1454)

* Working on cancellation detection

* Fixing typo

* Improving efficiency

* More efficiency and testing

* Reducing duplicated code in tests

* Finishing testing

* Adding additioanl tests and clean up

* Running pylint

* Adding constraint walker to diagnostics toolbox

* Addressing simple review comments

* Accumulate node rather than str(node)

* Clean up some tests

* Supporting ranged expressions

* Apply suggestions from code review

* Regorganising and fixing pylint issues

* Fixing false positives for linking constraints

* Fixing unused variable warning

* Fixing bug related to external functions with string arguments

* Addressing comments

* removing debugging prints

* Addressing more comments

* Debugging issue with nested external functions

* Adding todo

* Relax test tolerance due to Windows failures

* Adding catches for string arguments/values

* Test for nested external functions

* Fixing bug in walker logic for string arguments to external funcitons

* Address double counting of constraints

* Fixing typo

* Bumping Pyomo version

* Adding zero tolerance

* Adding test for zero tolerance

* Fixing typo

* Logic to skip scaled equality constraints

* Handling mole fraction type constraints

* Catching negations

* Fixing typos

* Limit of cancelling combinations

* Fixing typos

* Displaying cancelling terms to user

* More detail of mismatches

* Adding hooks and warnings about new config arguments

* Expression aware to string method

* Moving expression writer to util/misc

* Typos and comments

* Using StrEnum
---
 idaes/core/scaling/custom_scaler_base.py      |    4 +-
 idaes/core/util/misc.py                       |   66 +
 idaes/core/util/model_diagnostics.py          |  791 +++++++++-
 idaes/core/util/tests/test_misc.py            |   56 +-
 .../core/util/tests/test_model_diagnostics.py | 1377 ++++++++++++++++-
 setup.py                                      |    2 +-
 6 files changed, 2256 insertions(+), 40 deletions(-)

diff --git a/idaes/core/scaling/custom_scaler_base.py b/idaes/core/scaling/custom_scaler_base.py
index 8e94cc7c5a..a6b311086d 100644
--- a/idaes/core/scaling/custom_scaler_base.py
+++ b/idaes/core/scaling/custom_scaler_base.py
@@ -16,7 +16,6 @@
 Author: Andrew Lee
 """
 from copy import copy
-from enum import Enum
 
 from pyomo.environ import ComponentMap, units, value
 from pyomo.core.base.units_container import UnitsError
@@ -26,6 +25,7 @@
 from idaes.core.scaling.scaling_base import CONFIG, ScalerBase
 from idaes.core.scaling.util import get_scaling_factor, NominalValueExtractionVisitor
 import idaes.logger as idaeslog
+from idaes.core.util.misc import StrEnum
 
 # Set up logger
 _log = idaeslog.getLogger(__name__)
@@ -41,7 +41,7 @@
 }
 
 
-class ConstraintScalingScheme(str, Enum):
+class ConstraintScalingScheme(StrEnum):
     """
     Schemes available for calculating constraint scaling factors.
 
diff --git a/idaes/core/util/misc.py b/idaes/core/util/misc.py
index 5daa73676d..78185256cd 100644
--- a/idaes/core/util/misc.py
+++ b/idaes/core/util/misc.py
@@ -15,9 +15,12 @@
 This module contains miscellaneous utility functions for use in IDAES models.
 """
 from enum import Enum
+import sys
 
 import pyomo.environ as pyo
 from pyomo.common.config import ConfigBlock
+from pyomo.core.expr.visitor import _ToStringVisitor
+from pyomo.core.base.expression import ExpressionData
 
 import idaes.logger as idaeslog
 
@@ -177,3 +180,66 @@ class StrEnum(str, Enum):
 
     def __str__(self):
         return str(self.value)
+
+
+class _ToExprStringVisitor(_ToStringVisitor):
+    """
+    Derived version of the Pyomo _ToStringVisitor class
+    which checks for named Expressions in the expression tree
+    and prints their name instead of expanding the expression tree.
+    """
+
+    def visiting_potential_leaf(self, node):
+        # Check if node is a named Expression
+        if isinstance(node, ExpressionData):
+            # If it is, return the name of the Expression
+            return True, node.name
+
+        # Otherwise, continue descending as normal
+        return super().visiting_potential_leaf(node)
+
+
+def compact_expression_to_string(expr):
+    """Return a compact string representation of an expression.
+
+    Unlike the normal Pyomo string representations, this function will
+    identify any Expressions that appear within the expression tree and
+    print the name of these rather than expanding the expression tree.
+
+    Args:
+        expr: ExpressionBase. The root node of an expression tree.
+
+    Returns:
+        A string representation for the expression.
+    """
+    # Create and execute the visitor pattern
+    #
+    visitor = _ToExprStringVisitor(verbose=False, smap=None)
+    return visitor.dfs_postorder_stack(expr)
+
+
+def print_compact_form(expr, stream=None):
+    """
+    Writes a compact string representation of the given component or
+    expression to the specified stream.
+
+    Unlike the normal Pyomo string representations, this function will
+    identify any Expressions that appear within the expression tree and
+    print the name of these rather than expanding the expression tree.
+
+    Args:
+        expr: component or expression to print
+        stream: StringIO object to write to. Default is stdout.
+
+    Returns:
+        None
+    """
+    if stream is None:
+        stream = sys.stdout
+
+    if hasattr(expr, "expr"):
+        # We have a Constraint or Expression
+        # We want to print the expression, not the object name
+        expr = expr.expr
+
+    stream.write(compact_expression_to_string(expr))
diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py
index c7f5335c26..f701677d47 100644
--- a/idaes/core/util/model_diagnostics.py
+++ b/idaes/core/util/model_diagnostics.py
@@ -20,10 +20,11 @@
 from operator import itemgetter
 import sys
 from inspect import signature
-from math import log, isclose, inf, isfinite
+from math import log, log10, isclose, inf, isfinite
 import json
 from typing import List
 import logging
+from itertools import combinations, chain
 
 import numpy as np
 from scipy.linalg import svd
@@ -35,6 +36,7 @@
     Integers,
     Block,
     check_optimal_termination,
+    ComponentMap,
     ConcreteModel,
     Constraint,
     Expression,
@@ -58,6 +60,7 @@
 from pyomo.core.base.block import BlockData
 from pyomo.core.base.var import VarData
 from pyomo.core.base.constraint import ConstraintData
+from pyomo.core.base.expression import ExpressionData
 from pyomo.repn.standard_repn import (  # pylint: disable=no-name-in-module
     generate_standard_repn,
 )
@@ -67,6 +70,8 @@
     ConfigValue,
     document_kwargs_from_configdict,
     PositiveInt,
+    NonNegativeFloat,
+    NonNegativeInt,
 )
 from pyomo.util.check_units import identify_inconsistent_units
 from pyomo.contrib.incidence_analysis import IncidenceGraphInterface
@@ -78,8 +83,12 @@
 from pyomo.common.deprecation import deprecation_warning
 from pyomo.common.errors import PyomoException
 from pyomo.common.tempfiles import TempfileManager
+from pyomo.core import expr as EXPR
+from pyomo.common.numeric_types import native_types
+from pyomo.core.base.units_container import _PyomoUnit
 
 from idaes.core.solvers.get_solver import get_solver
+from idaes.core.util.misc import compact_expression_to_string
 from idaes.core.util.model_statistics import (
     activated_blocks_set,
     deactivated_blocks_set,
@@ -180,7 +189,6 @@ def svd_sparse(jacobian, number_singular_values):
     """
     u, s, vT = svds(jacobian, k=number_singular_values, which="SM")
 
-    print(u, s, vT, number_singular_values)
     return u, s, vT.transpose()
 
 
@@ -189,7 +197,7 @@ def svd_sparse(jacobian, number_singular_values):
     "variable_bounds_absolute_tolerance",
     ConfigValue(
         default=1e-4,
-        domain=float,
+        domain=NonNegativeFloat,
         description="Absolute tolerance for considering a variable to be close "
         "to its bounds.",
     ),
@@ -198,7 +206,7 @@ def svd_sparse(jacobian, number_singular_values):
     "variable_bounds_relative_tolerance",
     ConfigValue(
         default=1e-4,
-        domain=float,
+        domain=NonNegativeFloat,
         description="Relative tolerance for considering a variable to be close "
         "to its bounds.",
     ),
@@ -207,7 +215,7 @@ def svd_sparse(jacobian, number_singular_values):
     "variable_bounds_violation_tolerance",
     ConfigValue(
         default=0,
-        domain=float,
+        domain=NonNegativeFloat,
         description="Absolute tolerance for considering a variable to violate its bounds.",
         doc="Absolute tolerance for considering a variable to violate its bounds. "
         "Some solvers relax bounds on variables thus allowing a small violation to be "
@@ -218,15 +226,47 @@ def svd_sparse(jacobian, number_singular_values):
     "constraint_residual_tolerance",
     ConfigValue(
         default=1e-5,
-        domain=float,
+        domain=NonNegativeFloat,
         description="Absolute tolerance to use when checking constraint residuals.",
     ),
 )
+CONFIG.declare(
+    "constraint_term_mismatch_tolerance",
+    ConfigValue(
+        default=1e6,
+        domain=NonNegativeFloat,
+        description="Magnitude difference to use when checking for mismatched additive terms in constraints.",
+    ),
+)
+CONFIG.declare(
+    "constraint_term_cancellation_tolerance",
+    ConfigValue(
+        default=1e-4,
+        domain=NonNegativeFloat,
+        description="Absolute tolerance to use when checking for canceling additive terms in constraints.",
+    ),
+)
+CONFIG.declare(
+    "max_canceling_terms",
+    ConfigValue(
+        default=5,
+        domain=NonNegativeInt,
+        description="Maximum number of terms to consider when looking for canceling combinations in expressions.",
+    ),
+)
+CONFIG.declare(
+    "constraint_term_zero_tolerance",
+    ConfigValue(
+        default=1e-10,
+        domain=NonNegativeFloat,
+        description="Absolute tolerance to use when determining if a constraint term is equal to zero.",
+    ),
+)
 CONFIG.declare(
     "variable_large_value_tolerance",
     ConfigValue(
         default=1e4,
-        domain=float,
+        domain=NonNegativeFloat,
         description="Absolute tolerance for considering a value to be large.",
     ),
 )
@@ -234,7 +274,7 @@ def svd_sparse(jacobian, number_singular_values):
     "variable_small_value_tolerance",
     ConfigValue(
         default=1e-4,
-        domain=float,
+        domain=NonNegativeFloat,
         description="Absolute tolerance for considering a value to be small.",
     ),
 )
@@ -242,7 +282,7 @@ def svd_sparse(jacobian, number_singular_values):
     "variable_zero_value_tolerance",
     ConfigValue(
         default=1e-8,
-        domain=float,
+        domain=NonNegativeFloat,
         description="Absolute tolerance for considering a value to be near to zero.",
     ),
 )
@@ -250,7 +290,7 @@ def svd_sparse(jacobian, number_singular_values):
     "jacobian_large_value_caution",
     ConfigValue(
         default=1e4,
-        domain=float,
+        domain=NonNegativeFloat,
         description="Tolerance for raising a caution for large Jacobian values.",
     ),
 )
@@ -258,7 +298,7 @@ def svd_sparse(jacobian, number_singular_values):
     "jacobian_large_value_warning",
     ConfigValue(
         default=1e8,
-        domain=float,
+        domain=NonNegativeFloat,
         description="Tolerance for raising a warning for large Jacobian values.",
     ),
 )
@@ -266,7 +306,7 @@ def svd_sparse(jacobian, number_singular_values):
     "jacobian_small_value_caution",
     ConfigValue(
         default=1e-4,
-        domain=float,
+        domain=NonNegativeFloat,
         description="Tolerance for raising a caution for small Jacobian values.",
     ),
 )
@@ -274,7 +314,7 @@ def svd_sparse(jacobian, number_singular_values):
     "jacobian_small_value_warning",
     ConfigValue(
         default=1e-8,
-        domain=float,
+        domain=NonNegativeFloat,
         description="Tolerance for raising a warning for small Jacobian values.",
     ),
 )
@@ -290,7 +330,7 @@ def svd_sparse(jacobian, number_singular_values):
     "parallel_component_tolerance",
     ConfigValue(
         default=1e-8,
-        domain=float,
+        domain=NonNegativeFloat,
         description="Tolerance for identifying near-parallel Jacobian rows/columns",
     ),
 )
@@ -298,7 +338,7 @@ def svd_sparse(jacobian, number_singular_values):
     "absolute_feasibility_tolerance",
     ConfigValue(
         default=1e-6,
-        domain=float,
+        domain=NonNegativeFloat,
         description="Feasibility tolerance for identifying infeasible constraints and bounds",
     ),
 )
@@ -336,7 +376,7 @@ def svd_sparse(jacobian, number_singular_values):
     "singular_value_tolerance",
     ConfigValue(
         default=1e-6,
-        domain=float,
+        domain=NonNegativeFloat,
         description="Tolerance for defining a small singular value",
     ),
 )
@@ -344,7 +384,7 @@ def svd_sparse(jacobian, number_singular_values):
     "size_cutoff_in_singular_vector",
     ConfigValue(
         default=0.1,
-        domain=float,
+        domain=NonNegativeFloat,
         description="Size below which to ignore constraints and variables in "
         "the singular vector",
     ),
@@ -371,7 +411,7 @@ def svd_sparse(jacobian, number_singular_values):
     "M",  # TODO: Need better name
     ConfigValue(
         default=1e5,
-        domain=float,
+        domain=NonNegativeFloat,
         description="Maximum value for nu in MILP models.",
     ),
 )
@@ -379,7 +419,7 @@ def svd_sparse(jacobian, number_singular_values):
     "m_small",  # TODO: Need better name
     ConfigValue(
         default=1e-5,
-        domain=float,
+        domain=NonNegativeFloat,
         description="Smallest value for nu to be considered non-zero in MILP models.",
     ),
 )
@@ -387,7 +427,7 @@ def svd_sparse(jacobian, number_singular_values):
     "trivial_constraint_tolerance",
     ConfigValue(
         default=1e-6,
-        domain=float,
+        domain=NonNegativeFloat,
         description="Tolerance for identifying non-zero rows in Jacobian.",
     ),
 )
@@ -1052,6 +1092,239 @@ 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_constraint_mismatches(self, descend_into=True):
+        """
+        Call ConstraintTermAnalysisVisitor on all Constraints in model to walk expression
+        tree and collect any instances of sum expressions with mismatched terms or potential
+        cancellations.
+
+        Args:
+            descend_into: whether to descend_into child_blocks
+
+        Returns:
+            List of strings summarising constraints with mismatched terms
+            List of strings summarising constraints with cancellations
+            List of strings with constraint names where constraint contains no free variables
+        """
+        walker = ConstraintTermAnalysisVisitor(
+            term_mismatch_tolerance=self.config.constraint_term_mismatch_tolerance,
+            term_cancellation_tolerance=self.config.constraint_term_cancellation_tolerance,
+            term_zero_tolerance=self.config.constraint_term_zero_tolerance,
+            # for the high level summary, we only need to know if there are any cancellations,
+            # but don't need to find all of them
+            max_cancellations_per_node=1,
+            max_canceling_terms=self.config.max_canceling_terms,
+        )
+
+        mismatch = []
+        cancellation = []
+        constant = []
+
+        for c in self._model.component_data_objects(
+            Constraint, descend_into=descend_into
+        ):
+            _, expr_mismatch, expr_cancellation, expr_constant, _ = (
+                walker.walk_expression(c.expr)
+            )
+
+            if len(expr_mismatch) > 0:
+                mismatch.append(f"{c.name}: {len(expr_mismatch)} mismatched term(s)")
+
+            if len(expr_cancellation) > 0:
+                cancellation.append(
+                    f"{c.name}: {len(expr_cancellation)} potential canceling term(s)"
+                )
+
+            if expr_constant:
+                constant.append(c.name)
+
+        return mismatch, cancellation, constant
+
+    def display_constraints_with_mismatched_terms(self, stream=None):
+        """
+        Display constraints in model which contain additive terms of significantly different magnitude.
+
+        Args:
+            stream: I/O object to write report to (default = stdout)
+
+        Returns:
+            None
+
+        """
+        if stream is None:
+            stream = sys.stdout
+
+        mismatch, _, _ = self._collect_constraint_mismatches()
+
+        # Write the output
+        _write_report_section(
+            stream=stream,
+            lines_list=mismatch,
+            title="The following constraints have mismatched terms:",
+            end_line="Call display_problematic_constraint_terms(constraint) for more information.",
+            header="=",
+            footer="=",
+        )
+
+    def display_constraints_with_canceling_terms(self, stream=None):
+        """
+        Display constraints in model which contain additive terms which potentially cancel each other.
+
+        Note that this method looks a the current state of the constraint, and will flag terms as
+        cancelling if you have a form A == B + C where C is significantly smaller than A and B. In some
+        cases this behavior is intended, as C is a correction term which happens to be very
+        small at the current state. However, you should review these constraints to determine whether
+        the correction term is important for the situation you are modeling and consider removing the
+        term if it will never be significant.
+
+        Args:
+            stream: I/O object to write report to (default = stdout)
+
+        Returns:
+            None
+
+        """
+        if stream is None:
+            stream = sys.stdout
+
+        _, canceling, _ = self._collect_constraint_mismatches()
+
+        # Write the output
+        _write_report_section(
+            stream=stream,
+            lines_list=canceling,
+            title="The following constraints have canceling terms:",
+            end_line="Call display_problematic_constraint_terms(constraint) for more information.",
+            header="=",
+            footer="=",
+        )
+
+    def display_problematic_constraint_terms(
+        self, constraint, max_cancellations: int = 5, stream=None
+    ):
+        """
+        Display a summary of potentially problematic terms in a given constraint.
+
+        Note that this method looks a the current state of the constraint, and will flag terms as
+        cancelling if you have a form A == B + C where C is significantly smaller than A and B. In some
+        cases this behavior is intended, as C is a correction term which happens to be very
+        small at the current state. However, you should review these constraints to determine whether
+        the correction term is important for the situation you are modeling and consider removing the
+        term if it will never be significant.
+
+        Args:
+            constraint: ConstraintData object to be examined
+            max_cancellations: maximum number of cancellations per node before termination.
+                None = find all cancellations.
+            stream: I/O object to write report to (default = stdout)
+
+        Returns:
+            None
+
+        """
+        if stream is None:
+            stream = sys.stdout
+
+        # Check that constraint is of correct type to give useful error message
+        if not isinstance(constraint, ConstraintData):
+            # Wrong type, check if it is an indexed constraint
+            if isinstance(constraint, Constraint):
+                raise TypeError(
+                    f"{constraint.name} is an IndexedConstraint. Please provide "
+                    f"an individual element of {constraint.name} (ConstraintData) "
+                    "to be examined for problematic terms."
+                )
+            else:
+                # Not a constraint
+                raise TypeError(
+                    f"{constraint.name} is not an instance of a Pyomo Constraint."
+                )
+
+        walker = ConstraintTermAnalysisVisitor(
+            term_mismatch_tolerance=self.config.constraint_term_mismatch_tolerance,
+            term_cancellation_tolerance=self.config.constraint_term_cancellation_tolerance,
+            term_zero_tolerance=self.config.constraint_term_zero_tolerance,
+            max_cancellations_per_node=max_cancellations,
+            max_canceling_terms=self.config.max_canceling_terms,
+        )
+
+        _, expr_mismatch, expr_cancellation, _, tripped = walker.walk_expression(
+            constraint.expr
+        )
+
+        # Combine mismatches and cancellations into a summary list
+        issues = []
+        for k, v in expr_mismatch.items():
+            tag = " "
+            if isinstance(k, ExpressionData):
+                # For clarity, if the problem node is a named Expression, note this in output
+                tag = " Expression "
+            # Want to show full expression node plus largest and smallest magnitudes
+            issues.append(
+                f"Mismatch in{tag}{compact_expression_to_string(k)} (Max {v[0]}, Min {v[1]})"
+            )
+        # Collect summary of cancelling terms for user
+        # Walker gives us back a list of nodes with cancelling terms
+        for k, v in expr_cancellation.items():
+            # Each node may have multiple cancellations, these are given as separate tuples
+            tag = " "
+            if isinstance(k, ExpressionData):
+                # For clarity, if the problem node is a named Expression, note this in output
+                tag = " Expression "
+            for i in v:
+                # For each cancellation, iterate over contributing terms and write a summary
+                terms = ""
+                for j in i:
+                    if len(terms) > 0:
+                        terms += ", "
+                    # +1 to switch from 0-index to 1-index
+                    terms += f"{j[0]+1} ({j[1]})"
+                issues.append(
+                    f"Cancellation in{tag}{compact_expression_to_string(k)}. Terms {terms}"
+                )
+
+        if tripped:
+            end_line = (
+                f"Number of canceling terms per node limited to {max_cancellations}."
+            )
+        else:
+            end_line = None
+
+        # Write the output
+        _write_report_section(
+            stream=stream,
+            lines_list=issues,
+            title=f"The following terms in {constraint.name} are potentially problematic:",
+            end_line=end_line,
+            header="=",
+            footer="=",
+        )
+
+    def display_constraints_with_no_free_variables(self, stream=None):
+        """
+        Display constraints in model which contain no free variables.
+
+        Args:
+            stream: I/O object to write report to (default = stdout)
+
+        Returns:
+            None
+
+        """
+        if stream is None:
+            stream = sys.stdout
+
+        _, _, constant = self._collect_constraint_mismatches()
+
+        # Write the output
+        _write_report_section(
+            stream=stream,
+            lines_list=constant,
+            title="The following constraints have no free variables:",
+            header="=",
+            footer="=",
+        )
+
     def _collect_structural_warnings(
         self, ignore_evaluation_errors=False, ignore_unit_consistency=False
     ):
@@ -1335,6 +1608,28 @@ def _collect_numerical_cautions(self, jac=None, nlp=None):
                 cstring = "Variable"
             cautions.append(f"Caution: {len(none_value)} {cstring} with None value")
 
+        # Constraints with possible ill-posed terms
+        mismatch, cancellation, constant = self._collect_constraint_mismatches()
+        if len(mismatch) > 0:
+            cstring = "Constraints"
+            if len(mismatch) == 1:
+                cstring = "Constraint"
+            cautions.append(f"Caution: {len(mismatch)} {cstring} with mismatched terms")
+        if len(cancellation) > 0:
+            cstring = "Constraints"
+            if len(cancellation) == 1:
+                cstring = "Constraint"
+            cautions.append(
+                f"Caution: {len(cancellation)} {cstring} with potential cancellation of terms"
+            )
+        if len(constant) > 0:
+            cstring = "Constraints"
+            if len(constant) == 1:
+                cstring = "Constraint"
+            cautions.append(
+                f"Caution: {len(constant)} {cstring} with no free variables"
+            )
+
         # Extreme Jacobian rows and columns
         jac_col = extreme_jacobian_columns(
             jac=jac,
@@ -2777,8 +3072,6 @@ def eq_degenerate(m_dh, v):
                     # This variable does not appear in any constraint
                     return Constraint.Skip
 
-        m_dh.pprint()
-
         m_dh.degenerate = Constraint(m_dh.V, rule=eq_degenerate)
 
         # When y_pos = 1, nu >= m_small
@@ -4122,3 +4415,457 @@ def _collect_model_statistics(model):
     )
 
     return stats
+
+
+class ConstraintTermAnalysisVisitor(EXPR.StreamBasedExpressionVisitor):
+    """
+    Expression walker for checking Constraints for problematic terms.
+
+    This walker will walk the expression and look for summation terms
+    with mismatched magnitudes or potential cancellations.
+
+    Args:
+        term_mismatch_tolerance: tolerance to use when determining mismatched
+            terms
+        term_cancellation_tolerance: tolerance to use when identifying
+            possible cancellation of terms
+        term_zero_tolerance: tolerance for considering terms equal to zero
+        max_canceling_terms: maximum number of terms to consider when looking
+            for canceling combinations (None = consider all possible combinations)
+        max_cancellations_per_node: maximum number of cancellations to collect
+            for a single node. Collection will terminate when this many cancellations
+            have been identified (None = collect all cancellations)
+
+    Returns:
+        list of values for top-level summation terms
+        list of terms with mismatched magnitudes
+        list of terms with potential cancellations
+        bool indicating whether expression is a constant
+    """
+
+    def __init__(
+        self,
+        term_mismatch_tolerance: float = 1e6,
+        term_cancellation_tolerance: float = 1e-4,
+        term_zero_tolerance: float = 1e-10,
+        max_canceling_terms: int = 4,
+        max_cancellations_per_node: int = 5,
+    ):
+        super().__init__()
+
+        # Tolerance attributes
+        self._log_mm_tol = log10(term_mismatch_tolerance)
+        self._sum_tol = term_cancellation_tolerance
+        self._zero_tolerance = term_zero_tolerance
+        self._max_canceling_terms = max_canceling_terms
+        self._max_cancellations_per_node = max_cancellations_per_node
+
+        # Placeholders for collecting results
+        self.canceling_terms = ComponentMap()
+        self.mismatched_terms = ComponentMap()
+
+        # Flag for if cancellation collection hit limit
+        self._cancellation_tripped = False
+
+    def _get_value_for_sum_subexpression(self, child_data):
+        # child_data is a tuple, with the 0-th element being the node values
+        if isinstance(child_data[0][0], str):
+            # Values may be a list containing a string in some cases (e.g. external functions)
+            # Return the string in this case
+            return child_data[0][0]
+        return sum(i for i in child_data[0])
+
+    def _generate_combinations(self, inputs, equality=False):
+        # We want to test all combinations of terms for cancellation
+
+        # The number of combinations we check depends on whether this is an (in)equality
+        # expression or a sum node deeper in the expression tree.
+        # We expect (in)equalities to generally sum to 0 (0 == expr) thus we want to
+        # check if any subset of the sum terms sum to zero (i.e. are any terms unnecessary).
+        # For other sum nodes, we need to check for any combination of terms.
+
+        # Maximum number of terms to include in combinations
+        max_comb = len(inputs)
+        if equality:
+            # Subtract 1 if (in)equality node
+            max_comb += -1
+        # We also have a limit on the maximum number of terms to consider
+        if self._max_canceling_terms is not None:
+            max_comb = min(max_comb, self._max_canceling_terms)
+
+        # Single terms cannot cancel, thus we want all combinations of length 2 to max terms
+        # Note the need for +1 due to way range works
+        combo_range = range(2, max_comb + 1)
+
+        # Collect combinations of terms in an iterator
+        # In order to identify the terms in each cancellation, we will pair each value
+        # with its index in the input set as a tuple using enumerate
+        for i in chain.from_iterable(
+            combinations(enumerate(inputs), r) for r in combo_range
+        ):
+            # Yield each combination in the set
+            yield i
+
+    def _check_sum_cancellations(self, values_list, equality=False):
+        # First, strip any terms with value 0 as they do not contribute to cancellation
+        # We do this to keep the number of possible combinations as small as possible
+        stripped = [i for i in values_list if abs(i) >= self._zero_tolerance]
+
+        cancellations = []
+
+        if len(stripped) == 0:
+            # If the stripped list is empty, there are no non-zero terms
+            # We can stop here and return False as there are no possible cancellations
+            return cancellations
+
+        # For scaling of tolerance, we want to compare to the largest absolute value of
+        # the input values
+        max_value = abs(max(stripped, key=abs))
+
+        for i in self._generate_combinations(stripped, equality):
+            # Generate combinations will return a list of combinations of input terms
+            # each element of the list will be a 2-tuple representing a term in the
+            # input list with the first value being the position in the input set and
+            # the second being the value
+
+            # Check if the sum of values in the combination are below tolerance
+            if abs(sum(j[1] for j in i)) <= self._sum_tol * max_value:
+                # If so, record combination as canceling
+                cancellations.append(i)
+
+            # Terminate loop if we have reached the max cancellations to collect
+            if len(cancellations) >= self._max_cancellations_per_node:
+                self._cancellation_tripped = True
+                break
+
+        return cancellations
+
+    def _perform_checks(self, node, child_data):
+        # Perform checks for problematic expressions
+        # First, need to check to see if any child data is a list
+        # This indicates a sum expression
+        const = True
+
+        for d in child_data:
+            # We will check for canceling terms here, rather than the sum itself, to handle special cases
+            # We want to look for cases where a sum term results in a value much smaller
+            # than the terms of the sum
+            # Each element of child_data is a tuple where the 0-th element is the node values
+            if isinstance(d[0][0], str):
+                # Values may be a list containing a string in some cases (e.g. external functions)
+                # Skip if this is the case
+                pass
+            else:
+                for c in self._check_sum_cancellations(d[0]):
+                    if node in self.canceling_terms.keys():
+                        self.canceling_terms[node].append(c)
+                    else:
+                        self.canceling_terms[node] = [c]
+
+            # Expression is not constant if any child is not constant
+            # Element 1 is a bool indicating if the child node is constant
+            if not d[1]:
+                const = False
+
+        # Return any problematic terms found
+        return const
+
+    def _check_base_type(self, node):
+        if isinstance(node, VarData):
+            const = node.fixed
+        else:
+            const = True
+        return [value(node)], const, False
+
+    def _check_equality_expression(self, node, child_data):
+        # (In)equality expressions are a special case of sum expressions
+        # child_data has two elements; 0 is the LHS of the (in)equality and 1 is the RHS
+        # Each of these then contains three elements; 0 is a list of values for the sum components,
+        # 1 is a bool indicating if the child node term is constant, and 2 is a bool indicating if
+        # the child ode is a sum expression.
+
+        # First, to check for cancellation we need to negate one side of the expression
+        # mdata will contain the new set of child_data with negated values
+        mdata = []
+        # child_data[0][0] has the values of the LHS of the (in)equality, and we will negate these
+        vals = []
+        for j in child_data[0][0]:
+            vals.append(-j)
+        # Append the negated values along with whether the node is constant to mdata
+        mdata.append((vals, child_data[0][1]))
+        # child_data[1] is the RHS, so we take this as it appears and append to mdata
+        mdata.append(child_data[1])
+
+        # Next, call the method to check the sum expression
+        vals, const, _ = self._check_sum_expression(node, mdata)
+
+        # Next, we need to check for canceling terms.
+        # child_data[x][2] indicates if a node is a sum expression or not
+        if not child_data[0][2] and not child_data[1][2]:
+            # If both sides are not sum expressions, we have the form a == b
+            # Simple lLinking constraints do not need to be checked
+            pass
+        # Next, we can ignore any term that has already been flagged as mismatched
+        elif node in self.mismatched_terms.keys():
+            pass
+        # We can also ignore any case where one side of the (in)equality is constant
+        # I.e. if either child_node[x][1] is True
+        elif any(d[1] for d in mdata):
+            pass
+        else:
+            # Check for cancellation
+            # First, collect terms from both sides
+            # Note: outer loop comes first in list comprehension
+            t = [i for d in mdata for i in d[0]]
+
+            # Then check for cancellations
+            for c in self._check_sum_cancellations(t, equality=True):
+                if node in self.canceling_terms.keys():
+                    self.canceling_terms[node].append(c)
+                else:
+                    self.canceling_terms[node] = [c]
+
+        return vals, const, False
+
+    def _check_product_expr(self, node, child_data):
+        # We expect child_data to be a tuple of len 2 (a*b)
+        # If both or neither a and b are sums (xor), handle like other expressions
+        if not child_data[0][2] ^ child_data[1][2]:
+            return self._check_general_expr(node, child_data)
+        else:
+            # Here we have the case of a sum and a multiplier
+            # For this case, we will make the result look like a sum with the
+            # multiplier applied
+            # This is important for scaling of the sum terms, as cancellation should be
+            # Checked on the scaled value
+
+            # First, check if both terms are constants - if so we can just get the value of
+            # this node and move on.
+            if child_data[0][1] and child_data[1][1]:
+                return self._check_general_expr(node, child_data)
+
+            # The length of the values (child_data[i][0]) of the multiplier will be 1
+            # We can just iterate over all terms in both value lists and multiply
+            # each pair
+            vals = [i * j for i in child_data[0][0] for j in child_data[1][0]]
+
+            # Return the list of values, not constant, is a sum expression (apparent)
+            return vals, False, True
+
+    def _check_div_expr(self, node, child_data):
+        # We expect child_data to be a tuple of len 2 (a/b)
+        # If the numerator is not a sum, handle like general expression
+        # child_data[0] is numerator, child_data[1] is denominator
+        if not child_data[0][2]:
+            return self._check_general_expr(node, child_data)
+        else:
+            # If the numerator is a sum, we will treat this as if the division
+            # were applied to each term in the numerator separately
+            # This is important for scaling of the sum terms, as cancellation should be
+            # Checked on the scaled value
+
+            # First, check if the numerator is constant - if so we can just get the value of
+            # this node and move on.
+            if child_data[0][1]:
+                return self._check_general_expr(node, child_data)
+
+            # Next, we need to get the value of the denominator
+            denom = self._get_value_for_sum_subexpression(child_data[1])
+
+            try:
+                vals = [i / denom for i in child_data[0][0]]
+            except ZeroDivisionError:
+                raise ZeroDivisionError(
+                    f"Error in ConstraintTermAnalysisVisitor: found division with denominator of 0 "
+                    f"({str(node)})."
+                )
+
+            # Return the list of values, not constant, is a sum expression (apparent)
+            return vals, False, True
+
+    def _check_negation_expr(self, node, child_data):
+        # Here we need to defer checking of cancellations too due to different ways
+        # these can appear in an expression.
+        # We will simply negate all values on the child node values (child_data[0][0])
+        # and pass on the rest.
+        vals = [-i for i in child_data[0][0]]
+
+        # Return the list of values, not constant, is a sum expression (apparent)
+        return vals, child_data[0][1], child_data[0][2]
+
+    def _check_general_expr(self, node, child_data):
+        const = self._perform_checks(node, child_data)
+
+        try:
+            # pylint: disable-next=protected-access
+            val = node._apply_operation(
+                list(map(self._get_value_for_sum_subexpression, child_data))
+            )
+        except ValueError:
+            raise ValueError(
+                f"Error in ConstraintTermAnalysisVisitor: error evaluating {str(node)}."
+            )
+        except ZeroDivisionError:
+            raise ZeroDivisionError(
+                f"Error in ConstraintTermAnalysisVisitor: found division with denominator of 0 "
+                f"({str(node)})."
+            )
+        except Exception as err:
+            # Catch and re-raise any other error that occurs
+            _log.exception(f"Unexpected {err=}, {type(err)=}")
+            raise
+
+        return [val], const, False
+
+    def _check_other_expression(self, node, child_data):
+        const = self._perform_checks(node, child_data)
+
+        # First, need to get value of input terms, which may be sub-expressions
+        input_mag = []
+        for i in child_data:
+            input_mag.append(self._get_value_for_sum_subexpression(i))
+
+        # Next, create a copy of the function with expected magnitudes as inputs
+        newfunc = node.create_node_with_local_data(input_mag)
+
+        # Evaluate new function and return the value along with check results
+        return [value(newfunc)], const, False
+
+    def _check_ranged_expression(self, node, child_data):
+        # child_data should have 3 elements, LHS, middle term, and RHS
+        lhs_vals, lhs_const, _ = self._check_equality_expression(node, child_data[:2])
+        rhs_vals, rhs_const, _ = self._check_equality_expression(node, child_data[1:])
+
+        # Constant is a bit vague in terms ranged expressions.
+        # We will call the term constant only if all parts are constant
+        const = lhs_const and rhs_const
+
+        # For values, we need to avoid double counting the middle term
+        # Also for sign convention, we will negate the outer terms
+        vals = lhs_vals + [-rhs_vals[1]]
+
+        return vals, const, False
+
+    def _check_sum_expression(self, node, child_data):
+        # Sum expressions need special handling
+        # For sums, collect all child values into a list
+        vals = []
+        # We will check for cancellation in this node at the next level
+        # Pyomo is generally good at simplifying compound sums
+        const = True
+        # Collect data from child nodes
+        for d in child_data:
+            vals.append(self._get_value_for_sum_subexpression(d))
+
+            # Expression is not constant if any child is not constant
+            # Element 1 is a bool indicating if node is constant
+            if not d[1]:
+                const = False
+
+        # Check for mismatched terms
+        if len(vals) > 1:
+            absvals = [abs(v) for v in vals if abs(v) >= self._zero_tolerance]
+            if len(absvals) > 0:
+                vl = max(absvals)
+                vs = min(absvals)
+                if vl != vs:
+                    if vs == 0:
+                        diff = log10(vl)
+                    else:
+                        diff = log10(vl / vs)
+
+                    if diff >= self._log_mm_tol:
+                        self.mismatched_terms[node] = (vl, vs)
+
+        return vals, const, True
+
+    node_type_method_map = {
+        EXPR.EqualityExpression: _check_equality_expression,
+        EXPR.InequalityExpression: _check_equality_expression,
+        EXPR.RangedExpression: _check_ranged_expression,
+        EXPR.SumExpression: _check_sum_expression,
+        EXPR.NPV_SumExpression: _check_sum_expression,
+        EXPR.ProductExpression: _check_product_expr,
+        EXPR.MonomialTermExpression: _check_product_expr,
+        EXPR.NPV_ProductExpression: _check_product_expr,
+        EXPR.DivisionExpression: _check_div_expr,
+        EXPR.NPV_DivisionExpression: _check_div_expr,
+        EXPR.PowExpression: _check_general_expr,
+        EXPR.NPV_PowExpression: _check_general_expr,
+        EXPR.NegationExpression: _check_negation_expr,
+        EXPR.NPV_NegationExpression: _check_negation_expr,
+        EXPR.AbsExpression: _check_general_expr,
+        EXPR.NPV_AbsExpression: _check_general_expr,
+        EXPR.UnaryFunctionExpression: _check_general_expr,
+        EXPR.NPV_UnaryFunctionExpression: _check_general_expr,
+        EXPR.Expr_ifExpression: _check_other_expression,
+        EXPR.ExternalFunctionExpression: _check_other_expression,
+        EXPR.NPV_ExternalFunctionExpression: _check_other_expression,
+        EXPR.LinearExpression: _check_sum_expression,
+    }
+
+    def exitNode(self, node, data):
+        """
+        Method to call when exiting node to check for potential issues.
+        """
+        # Return [node values], constant (bool), sum expression (bool)
+        # first check if the node is a leaf
+        nodetype = type(node)
+
+        if nodetype in native_types:
+            return [node], True, False
+
+        node_func = self.node_type_method_map.get(nodetype, None)
+        if node_func is not None:
+            return node_func(self, node, data)
+
+        if not node.is_expression_type():
+            # this is a leaf, but not a native type
+            if nodetype is _PyomoUnit:
+                # Unit have no value, so return 1 as a placeholder
+                return [1], True, False
+
+            # Var or Param
+            return self._check_base_type(node)
+            # might want to add other common types here
+
+        # not a leaf - check if it is a named expression
+        if (
+            hasattr(node, "is_named_expression_type")
+            and node.is_named_expression_type()
+        ):
+            return self._check_other_expression(node, data)
+
+        raise TypeError(
+            f"An unhandled expression node type: {str(nodetype)} was encountered while "
+            f"analyzing constraint terms {str(node)}"
+        )
+
+    def walk_expression(self, expr):
+        """
+        Main method to call to walk an expression and return analysis.
+
+        Args:
+            expr - expression to be analyzed
+
+        Returns:
+            list of values of top-level additive terms
+            ComponentSet containing any mismatched terms
+            ComponentSet containing any canceling terms
+            Bool indicating whether expression is a constant
+        """
+        # Create new holders for collected terms
+        self.canceling_terms = ComponentMap()
+        self.mismatched_terms = ComponentMap()
+
+        # Call parent walk_expression method
+        vals, const, _ = super().walk_expression(expr)
+
+        # Return results
+        return (
+            vals,
+            self.mismatched_terms,
+            self.canceling_terms,
+            const,
+            self._cancellation_tripped,
+        )
diff --git a/idaes/core/util/tests/test_misc.py b/idaes/core/util/tests/test_misc.py
index 594b790e12..363fcd8df0 100644
--- a/idaes/core/util/tests/test_misc.py
+++ b/idaes/core/util/tests/test_misc.py
@@ -14,16 +14,19 @@
 This module contains miscellaneous utility functions for use in IDAES models.
 """
 
+from io import StringIO
 import pytest
 import re
 
-from pyomo.environ import ConcreteModel, Set, Block, Var, units
+from pyomo.environ import ConcreteModel, Constraint, Expression, Set, Block, Var, units
 from pyomo.common.config import ConfigBlock
 from pyomo.core.base.units_container import UnitsError
 
 from idaes.core.util.misc import (
     add_object_reference,
     set_param_from_config,
+    compact_expression_to_string,
+    print_compact_form,
 )
 import idaes.logger as idaeslog
 
@@ -361,3 +364,54 @@ def test_unitted_default(self, caplog):
             "b no units provided for parameter test_param - assuming "
             "default units" in caplog.text
         )
+
+
+class TestToExprStringVisitor:
+    @pytest.mark.unit
+    def test_compact_expression_to_string(self):
+        m = ConcreteModel()
+        m.v1 = Var()
+        m.v2 = Var()
+        m.v3 = Var()
+
+        m.e1 = Expression(expr=m.v1 + m.v2)
+        m.c1 = Constraint(expr=0 == m.v3 * m.e1)
+
+        str = compact_expression_to_string(m.c1.expr)
+        expected = "v3*e1  ==  0"
+
+        assert str == expected
+
+    @pytest.mark.unit
+    def test_print_compact_form_expr(self):
+        m = ConcreteModel()
+        m.v1 = Var()
+        m.v2 = Var()
+        m.v3 = Var()
+
+        m.e1 = Expression(expr=m.v1 + m.v2)
+        m.c1 = Constraint(expr=0 == m.v3 * m.e1)
+
+        expected = "v3*e1  ==  0"
+
+        stream = StringIO()
+        print_compact_form(m.c1.expr, stream=stream)
+
+        assert stream.getvalue() == expected
+
+    @pytest.mark.unit
+    def test_print_compact_form_component(self):
+        m = ConcreteModel()
+        m.v1 = Var()
+        m.v2 = Var()
+        m.v3 = Var()
+
+        m.e1 = Expression(expr=m.v1 + m.v2)
+        m.c1 = Constraint(expr=0 == m.v3 * m.e1)
+
+        expected = "v3*e1  ==  0"
+
+        stream = StringIO()
+        print_compact_form(m.c1, stream=stream)
+
+        assert stream.getvalue() == expected
diff --git a/idaes/core/util/tests/test_model_diagnostics.py b/idaes/core/util/tests/test_model_diagnostics.py
index ec5b671e34..e10f8a7c57 100644
--- a/idaes/core/util/tests/test_model_diagnostics.py
+++ b/idaes/core/util/tests/test_model_diagnostics.py
@@ -13,52 +13,67 @@
 """
 This module contains model diagnostic utility functions for use in IDAES (Pyomo) models.
 """
+from copy import deepcopy
 from io import StringIO
 import math
-import numpy as np
-import pytest
-import re
 import os
-from copy import deepcopy
+import re
+from unittest import TestCase
 
+import numpy as np
 from pandas import DataFrame
+import pytest
 
 from pyomo.environ import (
     Block,
     ConcreteModel,
     Constraint,
     Expression,
-    log,
-    tan,
-    asin,
-    acos,
-    sqrt,
+    Expr_if,
     Objective,
     PositiveIntegers,
     Set,
     SolverFactory,
-    Suffix,
-    TransformationFactory,
     units,
     value,
     Var,
     assert_optimal_termination,
     Param,
     Integers,
+    exp,
+    log,
+    log10,
+    sin,
+    cos,
+    tan,
+    asin,
+    acos,
+    atan,
+    sqrt,
+    sinh,
+    cosh,
+    tanh,
+    asinh,
+    acosh,
+    atanh,
 )
 from pyomo.common.collections import ComponentSet
 from pyomo.contrib.pynumero.asl import AmplInterface
 from pyomo.contrib.pynumero.interfaces.pyomo_nlp import PyomoNLP
 from pyomo.common.fileutils import this_file_dir
 from pyomo.common.tempfiles import TempfileManager
+from pyomo.core import expr as EXPR
 
 import idaes.core.util.scaling as iscale
 import idaes.logger as idaeslog
 from idaes.core.solvers import get_solver
 from idaes.core import FlowsheetBlock
 from idaes.core.util.testing import PhysicalParameterTestBlock
-from unittest import TestCase
-
+from idaes.models.properties.modular_properties.eos.ceos_common import (
+    cubic_roots_available,
+    CubicThermoExpressions,
+    CubicType as CubicEoS,
+)
 from idaes.core.util.model_diagnostics import (
     DiagnosticsToolbox,
     SVDToolbox,
@@ -81,6 +96,7 @@
     _collect_model_statistics,
     check_parallel_jacobian,
     compute_ill_conditioning_certificate,
+    ConstraintTermAnalysisVisitor,
 )
 from idaes.core.util.parameter_sweep import (
     SequentialSweepRunner,
@@ -90,6 +106,7 @@
     UniformSampling,
 )
 from idaes.core.util.testing import _enable_scip_solver_for_testing
+from idaes.models.properties import iapws95
 
 
 __author__ = "Alex Dowling, Douglas Allan, Andrew Lee"
@@ -999,6 +1016,232 @@ def test_display_near_parallel_variables(self):
     v1, v4
     v2, v4
 
+====================================================================================
+"""
+
+        assert stream.getvalue() == expected
+
+    @pytest.mark.component
+    def test_collect_constraint_mismatches(self):
+        m = ConcreteModel()
+        m.v1 = Var(initialize=2)
+        m.v2 = Var(initialize=3)
+
+        # Constraint with no free variables
+        m.c1 = Constraint(expr=m.v1 == m.v2)
+        m.v1.fix()
+        m.v2.fix()
+
+        # Constraint with mismatched terms
+        m.v3 = Var(initialize=10)
+        m.v4 = Var(initialize=10)
+        m.v5 = Var(initialize=1e-6)
+        m.c2 = Constraint(expr=m.v3 == m.v4 + m.v5)
+
+        # Constraint with cancellation
+        m.v6 = Var(initialize=10)
+        m.c3 = Constraint(expr=m.v6 == 10 + m.v3 - m.v4)
+
+        dt = DiagnosticsToolbox(model=m)
+
+        mismatch, cancellation, constant = dt._collect_constraint_mismatches()
+
+        assert mismatch == ["c2: 1 mismatched term(s)"]
+        assert cancellation == [
+            "c2: 1 potential canceling term(s)",
+            "c3: 1 potential canceling term(s)",
+        ]
+        assert constant == ["c1"]
+
+    @pytest.mark.component
+    def test_display_problematic_constraint_terms(self):
+        m = ConcreteModel()
+
+        # Constraint with mismatched terms
+        m.v3 = Var(initialize=10)
+        m.v4 = Var(initialize=10)
+        m.v5 = Var(initialize=1e-6)
+        m.c2 = Constraint(expr=m.v3 == m.v4 + m.v5)
+
+        dt = DiagnosticsToolbox(model=m)
+
+        stream = StringIO()
+        dt.display_problematic_constraint_terms(m.c2, stream=stream)
+
+        expected = """====================================================================================
+The following terms in c2 are potentially problematic:
+
+    Mismatch in v4 + v5 (Max 10, Min 1e-06)
+    Cancellation in v3  ==  v4 + v5. Terms 1 (-10), 2 (10)
+
+====================================================================================
+"""
+
+        assert stream.getvalue() == expected
+
+    @pytest.mark.component
+    def test_display_problematic_constraint_terms_limited(self):
+        m = ConcreteModel()
+        m.v1 = Var(initialize=10)
+        m.v2 = Var(initialize=-10)
+        m.v4 = Var(initialize=10)
+
+        # Constraint with cancellation
+        m.v6 = Var(initialize=10)
+        m.c1 = Constraint(expr=m.v6 == 10 - m.v4 + m.v1 + m.v2)
+
+        dt = DiagnosticsToolbox(model=m)
+
+        stream = StringIO()
+        dt.display_problematic_constraint_terms(m.c1, stream=stream)
+
+        expected = """====================================================================================
+The following terms in c1 are potentially problematic:
+
+    Cancellation in v6  ==  10 - v4 + v1 + v2. Terms 1 (-10), 2 (10)
+    Cancellation in v6  ==  10 - v4 + v1 + v2. Terms 1 (-10), 4 (10)
+    Cancellation in v6  ==  10 - v4 + v1 + v2. Terms 2 (10), 3 (-10)
+    Cancellation in v6  ==  10 - v4 + v1 + v2. Terms 2 (10), 5 (-10)
+    Cancellation in v6  ==  10 - v4 + v1 + v2. Terms 3 (-10), 4 (10)
+
+Number of canceling terms per node limited to 5.
+====================================================================================
+"""
+
+        assert stream.getvalue() == expected
+
+    @pytest.mark.component
+    def test_display_problematic_constraint_terms_indexed_error(self):
+        m = ConcreteModel()
+        m.s = Set(initialize=[1, 2])
+        m.v1 = Var(m.s, initialize=2)
+        m.v2 = Var(m.s, initialize=3)
+
+        def _c_rule(b, i):
+            return b.v1[i] == b.v2[i]
+
+        m.c1 = Constraint(m.s, rule=_c_rule)
+
+        dt = DiagnosticsToolbox(model=m)
+
+        # Check for indexed constraints
+        with pytest.raises(
+            TypeError,
+            match=re.escape(
+                "c1 is an IndexedConstraint. Please provide "
+                "an individual element of c1 (ConstraintData) "
+                "to be examined for problematic terms."
+            ),
+        ):
+            dt.display_problematic_constraint_terms(m.c1)
+
+        # Check for not a constraints
+        with pytest.raises(
+            TypeError, match="v1 is not an instance of a Pyomo Constraint."
+        ):
+            dt.display_problematic_constraint_terms(m.v1)
+
+    @pytest.mark.component
+    def test_display_problematic_constraint_terms_named_expression(self):
+        m = ConcreteModel()
+
+        # Constraint with mismatched terms
+        m.v1 = Var(initialize=10)
+        m.v2 = Var(initialize=10)
+        m.v3 = Var(initialize=1e-6)
+
+        m.e1 = Expression(expr=(m.v1 - m.v2))
+
+        m.c2 = Constraint(expr=m.v3 == m.e1**2)
+
+        dt = DiagnosticsToolbox(model=m)
+
+        stream = StringIO()
+        dt.display_problematic_constraint_terms(m.c2, stream=stream)
+
+        expected = """====================================================================================
+The following terms in c2 are potentially problematic:
+
+    Cancellation in Expression e1. Terms 1 (10), 2 (-10)
+
+====================================================================================
+"""
+
+        assert stream.getvalue() == expected
+
+    @pytest.mark.component
+    def test_display_constraints_with_mismatched_terms(self):
+        m = ConcreteModel()
+        # Constraint with mismatched terms
+        m.v3 = Var(initialize=10)
+        m.v4 = Var(initialize=10)
+        m.v5 = Var(initialize=1e-6)
+        m.c2 = Constraint(expr=m.v3 == m.v4 + m.v5)
+
+        dt = DiagnosticsToolbox(model=m)
+
+        stream = StringIO()
+        dt.display_constraints_with_mismatched_terms(stream=stream)
+
+        expected = """====================================================================================
+The following constraints have mismatched terms:
+
+    c2: 1 mismatched term(s)
+
+Call display_problematic_constraint_terms(constraint) for more information.
+====================================================================================
+"""
+
+        assert stream.getvalue() == expected
+
+    @pytest.mark.component
+    def test_display_constraints_with_canceling_terms(self):
+        m = ConcreteModel()
+        # Constraint with mismatched terms
+        m.v3 = Var(initialize=10)
+        m.v4 = Var(initialize=10)
+
+        # Constraint with cancellation
+        m.v6 = Var(initialize=10)
+        m.c3 = Constraint(expr=m.v6 == 10 + m.v3 - m.v4)
+
+        dt = DiagnosticsToolbox(model=m)
+
+        stream = StringIO()
+        dt.display_constraints_with_canceling_terms(stream=stream)
+
+        expected = """====================================================================================
+The following constraints have canceling terms:
+
+    c3: 1 potential canceling term(s)
+
+Call display_problematic_constraint_terms(constraint) for more information.
+====================================================================================
+"""
+
+        assert stream.getvalue() == expected
+
+    @pytest.mark.component
+    def test_display_constraints_with_no_free_variables(self):
+        m = ConcreteModel()
+        m.v1 = Var(initialize=2)
+        m.v2 = Var(initialize=3)
+
+        # Constraint with no free variables
+        m.c1 = Constraint(expr=m.v1 == m.v2)
+        m.v1.fix()
+        m.v2.fix()
+
+        dt = DiagnosticsToolbox(model=m)
+
+        stream = StringIO()
+        dt.display_constraints_with_no_free_variables(stream=stream)
+
+        expected = """====================================================================================
+The following constraints have no free variables:
+
+    c1
+
 ====================================================================================
 """
 
@@ -1151,6 +1394,7 @@ def test_collect_numerical_cautions(self, model):
         dt = DiagnosticsToolbox(model=model.b)
 
         cautions = dt._collect_numerical_cautions()
+
         assert len(cautions) == 5
         assert (
             "Caution: 2 Variables with value close to their bounds (abs=1.0E-04, rel=1.0E-04)"
@@ -1440,6 +1684,58 @@ def test_report_numerical_issues(self, model):
     compute_infeasibility_explanation()
     display_variables_at_or_outside_bounds()
 
+====================================================================================
+"""
+
+        assert stream.getvalue() == expected
+
+    @pytest.mark.component
+    def test_report_numerical_issues_cancellation(self):
+        model = ConcreteModel()
+
+        model.v1 = Var(initialize=1)
+        model.v2 = Var(initialize=2)
+        model.v3 = Var(initialize=1e-8)
+
+        # Non-problematic constraints
+        model.c1 = Constraint(expr=2 * model.v1 == model.v2)
+
+        # Problematic constraints
+        model.c10 = Constraint(expr=0 == exp(model.v1 - 0.5 * model.v2))
+        model.c11 = Constraint(expr=0 == model.v1 - 0.5 * model.v2 + model.v3)
+
+        dt = DiagnosticsToolbox(model=model)
+
+        stream = StringIO()
+        dt.report_numerical_issues(stream)
+
+        expected = """====================================================================================
+Model Statistics
+
+    Jacobian Condition Number: Undefined (Exactly Singular)
+
+------------------------------------------------------------------------------------
+3 WARNINGS
+
+    WARNING: 1 Constraint 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)
+
+------------------------------------------------------------------------------------
+3 Cautions
+
+    Caution: 1 Variable with value close to zero (tol=1.0E-08)
+    Caution: 1 Constraint with mismatched terms
+    Caution: 1 Constraint with potential cancellation of terms
+
+------------------------------------------------------------------------------------
+Suggested next steps:
+
+    display_constraints_with_large_residuals()
+    compute_infeasibility_explanation()
+    display_near_parallel_constraints()
+    display_near_parallel_variables()
+
 ====================================================================================
 """
 
@@ -2267,7 +2563,6 @@ def test_build_model(self, model):
         ca = IpoptConvergenceAnalysis(model)
 
         clone = ca._build_model()
-        clone.pprint()
 
         assert clone is not model
         assert isinstance(clone.v1, Var)
@@ -4053,3 +4348,1057 @@ def test_output(self, model):
 Constraints / bounds in guards for stability:
 """
         assert expected in stream.getvalue()
+
+
+class TestConstraintTermAnalysisVisitor:
+    @pytest.mark.unit
+    def test_sum_combinations(self):
+        # Check method to generate sums of all combinations of terms
+        # excludes single term sums
+        terms = [1, 2, 3, 4, 5]
+        visitor = ConstraintTermAnalysisVisitor(max_canceling_terms=None)
+        sums = [i for i in visitor._generate_combinations(terms)]
+
+        print(sums)
+
+        expected = [
+            ((0, 1), (1, 2)),
+            ((0, 1), (2, 3)),
+            ((0, 1), (3, 4)),
+            ((0, 1), (4, 5)),
+            ((1, 2), (2, 3)),
+            ((1, 2), (3, 4)),
+            ((1, 2), (4, 5)),
+            ((2, 3), (3, 4)),
+            ((2, 3), (4, 5)),
+            ((3, 4), (4, 5)),
+            ((0, 1), (1, 2), (2, 3)),
+            ((0, 1), (1, 2), (3, 4)),
+            ((0, 1), (1, 2), (4, 5)),
+            ((0, 1), (2, 3), (3, 4)),
+            ((0, 1), (2, 3), (4, 5)),
+            ((0, 1), (3, 4), (4, 5)),
+            ((1, 2), (2, 3), (3, 4)),
+            ((1, 2), (2, 3), (4, 5)),
+            ((1, 2), (3, 4), (4, 5)),
+            ((2, 3), (3, 4), (4, 5)),
+            ((0, 1), (1, 2), (2, 3), (3, 4)),
+            ((0, 1), (1, 2), (2, 3), (4, 5)),
+            ((0, 1), (1, 2), (3, 4), (4, 5)),
+            ((0, 1), (2, 3), (3, 4), (4, 5)),
+            ((1, 2), (2, 3), (3, 4), (4, 5)),
+            ((0, 1), (1, 2), (2, 3), (3, 4), (4, 5)),
+        ]
+
+        assert sums == expected
+
+    @pytest.mark.unit
+    def test_sum_combinations_limited_default(self):
+        # Check method to generate sums of all combinations of terms
+        # excludes single term sums
+        terms = [1, 2, 3, 4, 5]
+        visitor = ConstraintTermAnalysisVisitor()
+        sums = [i for i in visitor._generate_combinations(terms)]
+
+        expected = expected = [
+            ((0, 1), (1, 2)),
+            ((0, 1), (2, 3)),
+            ((0, 1), (3, 4)),
+            ((0, 1), (4, 5)),
+            ((1, 2), (2, 3)),
+            ((1, 2), (3, 4)),
+            ((1, 2), (4, 5)),
+            ((2, 3), (3, 4)),
+            ((2, 3), (4, 5)),
+            ((3, 4), (4, 5)),
+            ((0, 1), (1, 2), (2, 3)),
+            ((0, 1), (1, 2), (3, 4)),
+            ((0, 1), (1, 2), (4, 5)),
+            ((0, 1), (2, 3), (3, 4)),
+            ((0, 1), (2, 3), (4, 5)),
+            ((0, 1), (3, 4), (4, 5)),
+            ((1, 2), (2, 3), (3, 4)),
+            ((1, 2), (2, 3), (4, 5)),
+            ((1, 2), (3, 4), (4, 5)),
+            ((2, 3), (3, 4), (4, 5)),
+            ((0, 1), (1, 2), (2, 3), (3, 4)),
+            ((0, 1), (1, 2), (2, 3), (4, 5)),
+            ((0, 1), (1, 2), (3, 4), (4, 5)),
+            ((0, 1), (2, 3), (3, 4), (4, 5)),
+            ((1, 2), (2, 3), (3, 4), (4, 5)),
+            # 5 term combination should be excluded
+        ]
+
+        assert sums == expected
+
+    @pytest.mark.unit
+    def test_check_sum_cancellations(self):
+        terms = [1, -2, 3, -4, 5]
+        visitor = ConstraintTermAnalysisVisitor()
+        cancellations = visitor._check_sum_cancellations(terms)
+
+        # We expect to canceling combinations
+        # Results should be a list with each entry being a tuple containing a
+        # canceling combination
+        # In turn, each element of the tuple should be a 2-tuple with the
+        # position of the term in the input list and its value
+        expected = [((0, 1), (2, 3), (3, -4)), ((0, 1), (1, -2), (3, -4), (4, 5))]
+
+        assert cancellations == expected
+
+    @pytest.mark.unit
+    def test_int(self):
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=7)
+
+        assert vv == [7]
+        assert len(mm) == 0
+        assert len(cc) == 0
+        assert k
+        assert not tr
+
+    @pytest.mark.unit
+    def test_float(self):
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=7.7)
+
+        assert vv == [7.7]
+        assert len(mm) == 0
+        assert len(cc) == 0
+        assert k
+        assert not tr
+
+    @pytest.mark.unit
+    def test_scalar_param(self):
+        m = ConcreteModel()
+        m.scalar_param = Param(initialize=1, mutable=True)
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(
+            expr=m.scalar_param
+        )
+
+        assert vv == [1]
+        assert len(mm) == 0
+        assert len(cc) == 0
+        assert k
+        assert not tr
+
+    @pytest.mark.unit
+    def test_indexed_param(self):
+        m = ConcreteModel()
+        m.indexed_param = Param(["a", "b"], initialize=1, mutable=True)
+
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(
+            expr=m.indexed_param["a"]
+        )
+        assert vv == [1]
+        assert len(mm) == 0
+        assert len(cc) == 0
+        assert k
+        assert not tr
+
+    @pytest.mark.unit
+    def test_scalar_var(self):
+        m = ConcreteModel()
+        m.scalar_var = Var(initialize=1)
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(
+            expr=m.scalar_var
+        )
+
+        assert vv == [1]
+        assert len(mm) == 0
+        assert len(cc) == 0
+        assert not k
+        assert not tr
+
+    @pytest.mark.unit
+    def test_indexed_var(self):
+        m = ConcreteModel()
+        m.indexed_var = Var(["a", "b"], initialize=1)
+
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(
+            expr=m.indexed_var["a"]
+        )
+        assert vv == [1]
+        assert len(mm) == 0
+        assert len(cc) == 0
+        assert not k
+        assert not tr
+
+    @pytest.mark.unit
+    def test_scalar_var_fixed(self):
+        m = ConcreteModel()
+        m.scalar_var = Var(initialize=1)
+        m.scalar_var.fix()
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(
+            expr=m.scalar_var
+        )
+
+        assert vv == [1]
+        assert len(mm) == 0
+        assert len(cc) == 0
+        assert k
+        assert not tr
+
+    @pytest.mark.unit
+    def test_indexed_var_fixed(self):
+        m = ConcreteModel()
+        m.indexed_var = Var(["a", "b"], initialize=1)
+        m.indexed_var.fix()
+
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(
+            expr=m.indexed_var["a"]
+        )
+        assert vv == [1]
+        assert len(mm) == 0
+        assert len(cc) == 0
+        assert k
+        assert not tr
+
+    @pytest.mark.unit
+    def test_equality_expr(self):
+        m = ConcreteModel()
+        m.v1 = Var(initialize=1e-7)
+        m.v2 = Var(initialize=1e7)
+
+        expr = m.v1 == m.v2
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr)
+
+        assert vv == [-1e-7, 1e7]
+        assert expr in mm
+        assert len(mm) == 1
+        assert len(cc) == 0
+        assert not k
+        assert not tr
+
+    @pytest.mark.unit
+    def test_equality_expr_constant(self):
+        m = ConcreteModel()
+        m.v1 = Var(initialize=1e-7)
+        m.v2 = Var(initialize=1e7)
+
+        # Fix v1, not constant yet as v2 still free
+        m.v1.fix()
+
+        expr = m.v1 == m.v2
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr)
+
+        assert vv == [-1e-7, 1e7]
+        assert expr in mm
+        assert len(mm) == 1
+        assert len(cc) == 0
+        assert not k
+        assert not tr
+
+        # Fix v2, now constant
+        m.v2.fix()
+
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr)
+
+        assert vv == [-1e-7, 1e7]
+        assert expr in mm
+        assert len(mm) == 1
+        assert len(cc) == 0
+        assert k
+        assert not tr
+
+    @pytest.mark.unit
+    def test_sum_expr(self):
+        m = ConcreteModel()
+        m.v1 = Var(["a", "b", "c"], initialize=1e7)
+        m.v1["a"].set_value(1e-7)
+
+        expr = sum(m.v1[i] for i in m.v1)
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr)
+
+        assert vv == [1e-7, 1e7, 1e7]
+        assert expr in mm
+        assert len(mm) == 1
+        assert len(cc) == 0
+        assert not k
+        assert not tr
+
+    @pytest.mark.unit
+    def test_sum_expr_constant(self):
+        m = ConcreteModel()
+        m.v1 = Var(["a", "b", "c"], initialize=1e7)
+        m.v1["a"].set_value(1e-7)
+        m.v1.fix()
+
+        expr = sum(m.v1[i] for i in m.v1)
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr)
+
+        assert vv == [1e-7, 1e7, 1e7]
+        assert expr in mm
+        assert len(mm) == 1
+        assert len(cc) == 0
+        assert k
+        assert not tr
+
+    @pytest.fixture(scope="class")
+    def model(self):
+        m = ConcreteModel()
+        m.v1 = Var(initialize=2)
+        m.v2 = Var(initialize=3)
+        m.v3 = Var(initialize=5.0000001)
+        m.v4 = Var(initialize=5)
+
+        return m
+
+    @pytest.mark.unit
+    def test_product_expr(self, model):
+        m = ConcreteModel()
+        expr = model.v1 * model.v2 * model.v3
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr)
+
+        assert vv == [pytest.approx(30.0000006, rel=1e-8)]
+        assert len(mm) == 0
+        assert len(cc) == 0
+        assert not k
+        assert not tr
+
+    @pytest.mark.unit
+    def test_product_sum_expr(self, model):
+        expr = (model.v1 + model.v2) * (model.v3 + model.v4)
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr)
+
+        assert vv == [pytest.approx((2 + 3) * (5.0000001 + 5), rel=1e-8)]
+        assert len(mm) == 0
+        assert len(cc) == 0
+        assert not k
+        assert not tr
+
+    @pytest.mark.unit
+    def test_product_sum_expr_w_negation(self, model):
+        expr = (model.v1 + model.v2) * (model.v3 - model.v4)
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr)
+
+        assert vv == [pytest.approx(0.0000005, rel=1e-8)]
+        assert len(mm) == 0
+        assert expr in cc
+        assert len(cc) == 1
+        assert not k
+        assert not tr
+
+    @pytest.mark.unit
+    def test_division_expr(self, model):
+        expr = model.v1 / model.v2 / model.v3
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr)
+
+        assert vv == [pytest.approx(2 / 3 / 5.0000001, rel=1e-8)]
+        assert len(mm) == 0
+        assert len(cc) == 0
+        assert not k
+        assert not tr
+
+    @pytest.mark.unit
+    def test_division_sum_expr(self, model):
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(
+            expr=(model.v1 + model.v2) / (model.v3 + model.v4)
+        )
+
+        assert vv == [
+            pytest.approx(2 / (5.0000001 + 5), rel=1e-8),
+            pytest.approx(3 / (5.0000001 + 5), rel=1e-8),
+        ]
+        assert len(mm) == 0
+        assert len(cc) == 0
+        assert not k
+        assert not tr
+
+    @pytest.mark.unit
+    def test_division_sum_expr_w_negation(self, model):
+        expr = (model.v1 - model.v2) / (model.v3 - model.v4)
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr)
+
+        assert vv == [
+            pytest.approx(2 / 0.0000001, rel=1e-8),
+            pytest.approx(-3 / 0.0000001, rel=1e-8),
+        ]
+        assert len(mm) == 0
+        # Check for cancellation should be deferred
+        assert len(cc) == 0
+        assert not k
+        assert not tr
+
+    @pytest.mark.unit
+    def test_division_sum_expr_w_zero_denominator(self):
+        m = ConcreteModel()
+        m.v1 = Var(initialize=2)
+        m.v2 = Var(initialize=3)
+        m.v3 = Var(initialize=5)
+        m.v4 = Var(initialize=5)
+
+        expr = (m.v1 - m.v2) / (m.v3 - m.v4)
+        with pytest.raises(
+            ZeroDivisionError,
+            match=re.escape(
+                "Error in ConstraintTermAnalysisVisitor: found division with denominator of 0 "
+                "((v1 - v2)/(v3 - v4))."
+            ),
+        ):
+            ConstraintTermAnalysisVisitor().walk_expression(expr=expr)
+
+    @pytest.mark.unit
+    def test_division_sum_expr_w_negation_deferred(self):
+        m = ConcreteModel()
+        m.v1 = Var(initialize=2)
+        m.v2 = Var(initialize=2)
+        m.v3 = Var(initialize=5.0000001)
+        m.v4 = Var(initialize=5)
+
+        expr = ((m.v1 - m.v2) / (m.v3 - m.v4)) ** 2
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr)
+
+        assert vv == [pytest.approx(0, abs=1e-8)]
+        assert len(mm) == 0
+        assert expr in cc
+        assert len(cc) == 1
+        assert not k
+        assert not tr
+
+    @pytest.mark.unit
+    def test_division_expr_error(self):
+        m = ConcreteModel()
+        m.v1 = Var(initialize=2)
+        m.v2 = Var(initialize=0)
+
+        with pytest.raises(
+            ZeroDivisionError,
+            match=re.escape(
+                "Error in ConstraintTermAnalysisVisitor: found division with "
+                "denominator of 0 (v1/v2)."
+            ),
+        ):
+            ConstraintTermAnalysisVisitor().walk_expression(expr=m.v1 / m.v2)
+
+    @pytest.mark.unit
+    def test_pow_expr(self, model):
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(
+            expr=model.v1**model.v2
+        )
+
+        assert vv == [8]
+        assert len(mm) == 0
+        assert len(cc) == 0
+        assert not k
+        assert not tr
+
+    @pytest.mark.unit
+    def test_pow_sum_expr(self, model):
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(
+            expr=(model.v1 + model.v2) ** (model.v3 + model.v4)
+        )
+
+        assert vv == [pytest.approx(5**10.0000001, rel=1e-8)]
+        assert len(mm) == 0
+        assert len(cc) == 0
+        assert not k
+        assert not tr
+
+    @pytest.mark.unit
+    def test_pow_sum_expr_w_negation(self, model):
+        expr = (model.v1 + model.v2) ** (model.v3 - model.v4)
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr)
+
+        assert vv == [pytest.approx((2 + 3) ** (0.0000001), rel=1e-8)]
+        assert len(mm) == 0
+        assert expr in cc
+        assert len(cc) == 1
+        assert not k
+        assert not tr
+
+    @pytest.fixture(scope="class")
+    def func_model(self):
+        m = ConcreteModel()
+        m.v1 = Var(initialize=1)
+        m.v2 = Var(initialize=0.99999)
+        m.v3 = Var(initialize=-100)
+
+        m.v2.fix()
+
+        return m
+
+    @pytest.mark.unit
+    def test_negation_expr(self, func_model):
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(
+            expr=-func_model.v1
+        )
+
+        assert vv == [-1]
+        assert len(mm) == 0
+        assert len(cc) == 0
+        assert not k
+        assert not tr
+
+    @pytest.mark.unit
+    def test_negation_sum_expr(self, func_model):
+        expr = -(func_model.v1 - func_model.v2)
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr)
+
+        # Checking the cancellation should be deferred
+        # Expect to get two values back
+        assert vv == [pytest.approx(-1, rel=1e-8), pytest.approx(0.99999, rel=1e-8)]
+        assert len(mm) == 0
+        # Check is deferred, so no cancellations identified
+        assert len(cc) == 0
+        assert not k
+        assert not tr
+
+    # acosh has bounds that don't fit with other functions - we will assume we caught enough
+    func_list = [
+        exp,
+        log,
+        log10,
+        sqrt,
+        sin,
+        cos,
+        tan,
+        asin,
+        acos,
+        atan,
+        sinh,
+        cosh,
+        tanh,
+        asinh,
+        atanh,
+    ]
+    func_error_list = [log, log10, sqrt, asin, acos, acosh, atanh]
+
+    @pytest.mark.unit
+    @pytest.mark.parametrize("func", func_list)
+    def test_func_expr(self, func_model, func):
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(
+            expr=func(func_model.v2)
+        )
+
+        assert vv == [pytest.approx(value(func(0.99999)), rel=1e-8)]
+        assert len(mm) == 0
+        assert len(cc) == 0
+        assert k
+        assert not tr
+
+    @pytest.mark.unit
+    @pytest.mark.parametrize("func", func_list)
+    def test_func_sum_expr(self, func_model, func):
+        expr = func(func_model.v1 - func_model.v2)
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr)
+
+        assert vv == [pytest.approx(value(func(0.00001)), rel=1e-8)]
+        assert len(mm) == 0
+        assert expr in cc
+        assert len(cc) == 1
+        assert not k
+        assert not tr
+
+    @pytest.mark.unit
+    @pytest.mark.parametrize("func", func_list)
+    def test_func_sum_expr_w_negation(self, func_model, func):
+        expr = func(func_model.v1 - func_model.v2)
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr)
+
+        assert vv == [pytest.approx(value(func(0.00001)), rel=1e-8)]
+        assert len(mm) == 0
+        assert expr in cc
+        assert len(cc) == 1
+        assert not k
+        assert not tr
+
+    @pytest.mark.unit
+    @pytest.mark.parametrize("func", func_error_list)
+    def test_func_expr_error(self, func_model, func):
+        with pytest.raises(
+            ValueError,
+            match=re.escape(
+                "Error in ConstraintTermAnalysisVisitor: error evaluating "
+            ),
+        ):
+            ConstraintTermAnalysisVisitor().walk_expression(expr=func(func_model.v3))
+
+    @pytest.mark.unit
+    def test_expr_if(self, model):
+        model.exprif = Expr_if(
+            IF=model.v1,
+            THEN=model.v1 + model.v2,
+            ELSE=model.v3 - model.v4,
+        )
+
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(
+            expr=model.exprif
+        )
+
+        assert vv == [pytest.approx(5, rel=1e-8)]
+        assert len(mm) == 0
+        assert model.exprif in cc
+        assert len(cc) == 1
+        assert not k
+        assert not tr
+
+    @pytest.mark.unit
+    @pytest.mark.skipif(
+        not AmplInterface.available(), reason="pynumero_ASL is not available"
+    )
+    @pytest.mark.skipif(not cubic_roots_available, reason="Cubic roots not available")
+    def test_ext_func(self):
+        # Use the cubic root external function to test
+        m = ConcreteModel()
+        m.a = Var(initialize=1)
+        m.b = Var(initialize=1)
+
+        m.expr_write = CubicThermoExpressions(m)
+        Z = m.expr_write.z_liq(eos=CubicEoS.PR, A=m.a, B=m.b)
+
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=Z)
+
+        assert vv == [pytest.approx(-2.1149075414767577, rel=1e-8)]
+        assert len(mm) == 0
+        assert Z in cc
+        assert len(cc) == 1
+        assert not k
+        assert not tr
+
+        # Check that model state did not change
+        assert value(m.a) == 1
+        assert value(m.b) == 1
+        assert value(Z) == pytest.approx(-2.1149075414767577, rel=1e-8)
+
+    @pytest.mark.unit
+    def test_equality_sum_expr(self):
+        m = ConcreteModel()
+        m.v1 = Var(["a", "b", "c"], initialize=1e7)
+        m.v2 = Var(initialize=1e-7)
+
+        expr = m.v2 == sum(m.v1[i] for i in m.v1)
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr)
+
+        assert vv == [-1e-7, 3e7]
+        assert expr in mm
+        assert len(mm) == 1
+        assert len(cc) == 0
+        assert not k
+        assert not tr
+
+    @pytest.mark.unit
+    def test_inequality_sum_expr(self):
+        m = ConcreteModel()
+        m.v1 = Var(["a", "b", "c"], initialize=1e7)
+        m.v2 = Var(initialize=1e-7)
+
+        expr = m.v2 <= sum(m.v1[i] for i in m.v1)
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr)
+
+        assert vv == [-1e-7, 3e7]
+        assert expr in mm
+        assert len(mm) == 1
+        assert len(cc) == 0
+        assert not k
+        assert not tr
+
+    @pytest.mark.unit
+    def test_compound_equality_expr_1(self):
+        m = ConcreteModel()
+        m.v1 = Var(["a", "b", "c"], initialize=1e7)
+        m.v2 = Var(initialize=1e-7)
+
+        expr = 6 * m.v2 == 8 * sum(m.v1[i] for i in m.v1)
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr)
+
+        assert vv == [-6e-7, 2.4e8]
+        assert expr in mm
+        assert len(mm) == 1
+        assert len(cc) == 0
+        assert not k
+        assert not tr
+
+    @pytest.mark.unit
+    def test_ranged_expr(self):
+        m = ConcreteModel()
+        m.v1 = Var(initialize=1e7)
+        m.v2 = Var(initialize=1e-7)
+        m.v3 = Var(initialize=1e7)
+
+        m.expr = EXPR.RangedExpression(args=(m.v1, m.v2, m.v3), strict=True)
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=m.expr)
+
+        assert vv == [-1e7, 1e-7, -1e7]
+        assert m.expr in mm
+        assert len(mm) == 1
+        assert len(cc) == 0
+        assert not k
+        assert not tr
+
+        # Fix v1 and v2 to make first two terms constant
+        m.v1.fix()
+        m.v2.fix()
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=m.expr)
+
+        # Should not be flagged as constant due to v3
+        assert vv == [-1e7, 1e-7, -1e7]
+        assert m.expr in mm
+        assert len(mm) == 1
+        assert len(cc) == 0
+        assert not k
+        assert not tr
+
+        # Fix v3 to make all terms constant
+        m.v3.fix()
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=m.expr)
+
+        # Should now be constant
+        assert vv == [-1e7, 1e-7, -1e7]
+        assert m.expr in mm
+        assert len(mm) == 1
+        assert len(cc) == 0
+        assert k
+        assert not tr
+
+    @pytest.mark.unit
+    def test_compound_equality_expr_2(self):
+        m = ConcreteModel()
+        m.v1 = Var(["a", "b", "c"], initialize=1e7)
+        m.v2 = Var(initialize=1e-7)
+        m.v3 = Var(initialize=1e3)
+
+        # Set this small so we get two mismatched warnings
+        m.v1["a"].set_value(1e-7)
+
+        expr1 = sum(m.v1[i] for i in m.v1)
+        expr = 6 * m.v2 == 8 * expr1 + m.v3
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr)
+
+        assert vv == [-6e-7, pytest.approx(8 * (2e7 + 1e-7) + 1000, rel=1e-8)]
+        assert expr in mm
+        assert expr1 in mm
+        assert len(mm) == 2
+        assert len(cc) == 0
+        assert not k
+        assert not tr
+
+    @pytest.mark.unit
+    def test_canceling_sum_expr(self):
+        m = ConcreteModel()
+        m.v1 = Var(initialize=2)
+        m.v2 = Var(initialize=2)
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(
+            expr=m.v1 - m.v2
+        )
+
+        assert vv == [2, -2]
+        assert len(mm) == 0
+        # We do not check cancellation at the sum, so this should be empty
+        assert len(cc) == 0
+        assert not k
+        assert not tr
+
+    @pytest.mark.unit
+    def test_expr_w_canceling_sum(self):
+        m = ConcreteModel()
+        m.v1 = Var(initialize=2)
+        m.v2 = Var(initialize=2)
+        m.v3 = Var(initialize=3)
+
+        expr = m.v3 * (m.v1 - m.v2)
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr)
+
+        assert vv == [6, -6]
+        assert len(mm) == 0
+        # Check for cancellation should be deferred
+        assert len(cc) == 0
+        assert not k
+        assert not tr
+
+    @pytest.mark.unit
+    def test_expr_w_deferred_canceling_sum(self):
+        m = ConcreteModel()
+        m.v1 = Var(initialize=2)
+        m.v2 = Var(initialize=2)
+        m.v3 = Var(initialize=3)
+
+        expr = (m.v3 * (m.v1 - m.v2)) ** 2
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr)
+
+        assert vv == [0]
+        assert len(mm) == 0
+        # We should get a warning about canceling sums here
+        assert expr in cc
+        assert len(cc) == 1
+        assert not k
+        assert not tr
+
+        # Check for tolerance of sum cancellation
+        m.v2.set_value(2.00022)
+
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor(
+            term_cancellation_tolerance=1e-4
+        ).walk_expression(expr=expr)
+
+        assert vv == [pytest.approx((3 * -0.00022) ** 2, rel=1e-8)]
+        assert len(mm) == 0
+        # This should pass as the difference is greater than tol
+        assert len(cc) == 0
+        assert not k
+        assert not tr
+
+        # Change tolerance so it should identify cancellation
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor(
+            term_cancellation_tolerance=1e-3
+        ).walk_expression(expr=expr)
+
+        assert vv == [pytest.approx((3 * -0.00022) ** 2, rel=1e-8)]
+        assert len(mm) == 0
+        assert expr in cc
+        assert len(cc) == 1
+        assert not k
+        assert not tr
+
+    @pytest.mark.unit
+    def test_canceling_equality_expr_safe(self):
+        m = ConcreteModel()
+        m.v1 = Var(initialize=1e7)
+        m.v2 = Var(initialize=1e7)
+
+        # This is a standard constraint form, so should have no issues despite cancellation
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(
+            expr=0 == m.v1 - m.v2
+        )
+
+        assert vv == [0, pytest.approx(0, abs=1e-12)]
+        assert len(mm) == 0
+        assert len(cc) == 0
+        assert not k
+        assert not tr
+
+    @pytest.mark.unit
+    def test_canceling_equality_expr_zero_term(self):
+        m = ConcreteModel()
+        m.v1 = Var(initialize=1e7)
+        m.v2 = Var(initialize=1e7)
+        m.v3 = Var(initialize=0)
+
+        expr = m.v3 == m.v1 - m.v2
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr)
+
+        assert vv == [0, pytest.approx(0, abs=1e-12)]
+        assert len(mm) == 0
+        # No canceling terms as v3=0 is ignored thus we have a=b form
+        assert len(cc) == 0
+        assert not k
+        assert not tr
+
+        # Set v3 above zero tolerance
+        m.v3.set_value(1e-4)
+
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr)
+        assert vv == [-1e-4, pytest.approx(0, abs=1e-12)]
+        assert len(mm) == 0
+        assert expr in cc
+        assert len(cc) == 1
+        assert not k
+        assert not tr
+
+    @pytest.mark.unit
+    def test_canceling_equality_expr_compound(self):
+        m = ConcreteModel()
+        m.v1 = Var(["a", "b"], initialize=5e6)
+        m.v2 = Var(initialize=1e7)
+        m.v3 = Var(initialize=0)
+
+        expr = m.v3 == sum(v for v in m.v1.values()) - m.v2
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr)
+
+        assert vv == [0, pytest.approx(0, abs=1e-12)]
+        assert len(mm) == 0
+        # No canceling terms as v3=0 is ignored thus we have a=b form
+        assert len(cc) == 0
+        assert not k
+        assert not tr
+
+        # Set v3 above zero tolerance
+        m.v3.set_value(1e-4)
+
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr)
+        assert vv == [-1e-4, pytest.approx(0, abs=1e-12)]
+        assert len(mm) == 0
+        assert expr in cc
+        assert len(cc) == 1
+        assert not k
+        assert not tr
+
+    # Double check for a+eps=c form gets flagged in some way
+    @pytest.mark.unit
+    def test_canceling_equality_expr_canceling_sides(self):
+        m = ConcreteModel()
+        m.v1 = Var(initialize=1)
+        m.v2 = Var(initialize=1)
+        m.v3 = Var(initialize=1e-8)
+
+        expr1 = m.v1 + m.v3
+        expr = m.v2 == expr1
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr)
+
+        assert vv == [-1, pytest.approx(1 + 1e-8, abs=1e-8)]
+        assert expr1 in mm
+        assert len(mm) == 1
+        assert expr in cc
+        assert len(cc) == 1
+        assert not k
+        assert not tr
+
+    # Check to make sure simple linking constraints are not flagged as canceling
+    @pytest.mark.unit
+    def test_linking_equality_expr(self):
+        m = ConcreteModel()
+        m.v1 = Var(initialize=1)
+        m.v2 = Var(initialize=1)
+
+        expr = m.v1 == m.v2
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr)
+
+        assert vv == [-1, 1]
+        assert len(mm) == 0
+        assert len(cc) == 0
+        assert not k
+        assert not tr
+
+    @pytest.mark.unit
+    def test_linking_equality_expr_compound(self):
+        m = ConcreteModel()
+        m.v1 = Var(initialize=1)
+        m.v2 = Var(initialize=1)
+        m.v3 = Var(initialize=1)
+
+        expr = m.v1 == m.v2 * m.v3
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr)
+
+        assert vv == [-1, 1]
+        assert len(mm) == 0
+        assert len(cc) == 0
+        assert not k
+        assert not tr
+
+    @pytest.mark.component
+    def test_external_function_w_string_argument(self):
+        m = ConcreteModel()
+        m.properties = iapws95.Iapws95ParameterBlock()
+        m.state = m.properties.build_state_block([0])
+
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(
+            expr=m.state[0].enth_mol
+        )
+
+        assert vv == [pytest.approx(1.1021387e-2, rel=1e-6)]
+        assert len(mm) == 0
+        assert len(cc) == 0
+        assert not k
+        assert not tr
+
+        # Test nested external functions
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(
+            expr=m.state[0].temperature
+        )
+
+        assert vv == [pytest.approx(270.4877, rel=1e-6)]
+        assert len(mm) == 0
+        assert len(cc) == 0
+        assert not k
+        assert not tr
+
+    @pytest.mark.unit
+    def test_zero_tolerance(self):
+        m = ConcreteModel()
+        m.v1 = Var(initialize=1e-12)
+        m.v2 = Var(initialize=1)
+        m.v3 = Var(initialize=1e-12)
+
+        expr1 = m.v1 - m.v3
+        expr2 = expr1 + 1
+        expr = m.v2 == expr2
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr)
+
+        assert vv == [-1, pytest.approx(1, abs=1e-8)]
+        # We expect no mismatches, as smallest terms are below zero tolerance
+        assert len(mm) == 0
+        # We expect no canceling terms, as v1 and v3 are ignored (zero tolerance), leaving
+        # 1 == 1
+        assert len(cc) == 0
+        assert not k
+        assert not tr
+
+        # Set v1 and v3 above zero tolerance
+        m.v1.set_value(1e-8)
+        m.v3.set_value(1e-8)
+
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr)
+
+        assert vv == [-1, pytest.approx(1, abs=1e-8)]
+        assert expr2 in mm
+        assert len(mm) == 1
+        assert expr in cc
+        assert len(cc) == 1
+        assert not k
+        assert not tr
+
+    # Test to make sure scaled constraints are not flagged as issues
+    @pytest.mark.unit
+    def test_scaled_equality_expr_product(self):
+        m = ConcreteModel()
+        m.v1 = Var(initialize=1)
+        m.v2 = Var(initialize=1)
+        m.v3 = Var(initialize=2)
+
+        expr = 0 == m.v3 * (m.v1 - m.v2)
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr)
+
+        assert vv == [0, 0]
+        assert len(mm) == 0
+        assert len(cc) == 0
+        assert not k
+        assert not tr
+
+    @pytest.mark.unit
+    def test_scaled_equality_expr_division(self):
+        m = ConcreteModel()
+        m.v1 = Var(initialize=1)
+        m.v2 = Var(initialize=1)
+        m.v3 = Var(initialize=2)
+
+        expr = 0 == (m.v1 - m.v2) / m.v3
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr)
+
+        assert vv == [0, 0]
+        assert len(mm) == 0
+        assert len(cc) == 0
+        assert not k
+        assert not tr
+
+    @pytest.mark.unit
+    def test_mole_fraction_constraint(self):
+        m = ConcreteModel()
+        m.mole_frac_a = Var(initialize=0.5)
+        m.flow_a = Var(initialize=100)
+        m.flow_b = Var(initialize=100)
+
+        expr = m.mole_frac_a * (m.flow_a + m.flow_b) == m.flow_a
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr)
+
+        assert vv == [pytest.approx(-100, rel=1e-5), pytest.approx(100, rel=1e-5)]
+        assert len(mm) == 0
+        assert len(cc) == 0
+        assert not k
+        assert not tr
+
+    @pytest.mark.unit
+    def test_mole_fraction_constraint_trace(self):
+        m = ConcreteModel()
+        m.mole_frac_a = Var(initialize=0.999810)
+        m.flow_a = Var(initialize=122.746)
+        m.flow_b = Var(initialize=0.0233239)
+
+        expr = m.mole_frac_a * (m.flow_a + m.flow_b) == m.flow_a
+        vv, mm, cc, k, tr = ConstraintTermAnalysisVisitor().walk_expression(expr=expr)
+
+        assert vv == [
+            pytest.approx(-122.746, rel=1e-5),
+            pytest.approx(122.746, rel=1e-5),
+        ]
+        assert len(mm) == 0
+        assert len(cc) == 0
+        assert not k
+        assert not tr
diff --git a/setup.py b/setup.py
index d5303b756a..19868f5d1d 100644
--- a/setup.py
+++ b/setup.py
@@ -82,7 +82,7 @@ def __getitem__(self, key):
     # Put abstract (non-versioned) deps here.
     # Concrete dependencies go in requirements[-dev].txt
     install_requires=[
-        "pyomo >= 6.8.0",
+        "pyomo @ https://github.com/IDAES/pyomo/archive/6.8.0.idaes.2024.10.25.zip",
         "pint >= 0.24.1",  # required to use Pyomo units. Pint 0.24.1 needed for Python 3.9 support
         "networkx",  # required to use Pyomo network
         "numpy>=1,<3",

From 7eac46139444f0cad681ede7154151123726f13b Mon Sep 17 00:00:00 2001
From: Miranda Mundt <55767766+mrmundt@users.noreply.github.com>
Date: Mon, 11 Nov 2024 18:36:13 -0700
Subject: [PATCH 2/2] New URL for Downloads Badge (#1526)

---
 README.md | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/README.md b/README.md
index 8309621c00..2d5ee719b9 100644
--- a/README.md
+++ b/README.md
@@ -14,7 +14,7 @@ scale-up, operation and troubleshooting of innovative, advanced energy systems.
 [![GitHub contributors](https://img.shields.io/github/contributors/IDAES/idaes-pse.svg)](https://github.com/IDAES/idaes-pse/graphs/contributors)
 [![Merged PRs](https://img.shields.io/github/issues-pr-closed-raw/IDAES/idaes-pse.svg?label=merged+PRs)](https://github.com/IDAES/idaes-pse/pulls?q=is:pr+is:merged)
 [![Issue stats](http://isitmaintained.com/badge/resolution/IDAES/idaes-pse.svg)](http://isitmaintained.com/project/IDAES/idaes-pse)
-[![Downloads](https://pepy.tech/badge/idaes-pse)](https://pepy.tech/project/idaes-pse)
+[![Downloads](https://static.pepy.tech/badge/idaes-pse)](https://pepy.tech/project/idaes-pse)
 <!-- END Status badges -->
 
 ## Getting Started