diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index c7f5335c26..56f36087ed 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 @@ -67,6 +68,7 @@ ConfigValue, document_kwargs_from_configdict, PositiveInt, + NonNegativeFloat, ) from pyomo.util.check_units import identify_inconsistent_units from pyomo.contrib.incidence_analysis import IncidenceGraphInterface @@ -78,6 +80,9 @@ 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.model_statistics import ( @@ -189,7 +194,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 +203,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 +212,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 +223,31 @@ 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( "variable_large_value_tolerance", ConfigValue( default=1e4, - domain=float, + domain=NonNegativeFloat, description="Absolute tolerance for considering a value to be large.", ), ) @@ -234,7 +255,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 +263,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 +271,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 +279,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 +287,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 +295,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 +311,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 +319,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 +357,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 +365,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 +392,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 +400,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 +408,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 +1073,107 @@ 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): + walker = ConstraintTermAnalysisVisitor( + term_mismatch_tolerance=self.config.constraint_term_mismatch_tolerance, + term_cancellation_tolerance=self.config.constraint_term_cancellation_tolerance, + ) + + 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 + ) + + for i in expr_mismatch: + mismatch.append(f"{c.name}: {i}") + for i in expr_cancellation: + cancellation.append(f"{c.name}: {i}") + 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:", + header="=", + footer="=", + ) + + def display_constraints_with_canceling_terms(self, stream=None): + """ + Display constraints in model which contain additive terms which potentially cancel each other. + + 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:", + 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 +1457,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, @@ -4122,3 +4266,297 @@ 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 + + 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, + ): + super().__init__() + + # Tolerance attributes + self._log_mm_tol = log10(term_mismatch_tolerance) + self._sum_tol = term_cancellation_tolerance + + # Placeholders for collecting results + self.canceling_terms = ComponentSet() + self.mismatched_terms = ComponentSet() + + def _get_value_for_sum_subexpression(self, child_data): + return sum(i for i in child_data[0]) + + def _generate_sum_combinations(self, inputs): + # We want to test all combinations of terms for cancellation + # Single terms cannot cancel, thus we want all combinations of length 2 to all terms + combo_range = range(2, len(inputs) + 1) + for i in chain.from_iterable(combinations(inputs, r) for r in combo_range): + # We want the absolute value of the sum of the combination + yield abs(sum(i)) + + def _check_sum_cancellations(self, values_list): + # 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 i != 0] + + 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 False + + # 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_sum_combinations(stripped): + # We do not need to check all combinations, as we only need determine if any + # combination cancels. Thus, if we find a cancellation, stop and return True + if i <= self._sum_tol * max_value: + return True + + return False + + 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 + if self._check_sum_cancellations(d[0]): + self.canceling_terms.add(node) + + # Expression is not constant if any child is not 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 + # First, need to negate one side of the expression + mdata = [] + + vals = [] + for j in child_data[0][0]: + vals.append(-j) + mdata.append((vals, child_data[0][1])) + 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. + # In this case, we can safely ignore expressions of the form constant = sum() + # We can also ignore any constraint that is already flagged as mismatched + # We will also ignore any constraints where a == b and neither a nor b are sum expressions + if not child_data[0][2] and not child_data[1][2]: + # Linking constraint, skip checks + pass + elif node not in self.mismatched_terms and not any(d[1] for d in mdata): + # No constant terms, 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 + if self._check_sum_cancellations(t): + self.canceling_terms.add(node) + + return vals, const, False + + 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 = [self._get_value_for_sum_subexpression(i) for i in child_data] + + # 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): + 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 + if not d[1]: + const = False + + # Check for mismatched terms + if len(vals) > 1: + absvals = [abs(v) for v in vals] + 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.add(node) + + 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_general_expr, + EXPR.MonomialTermExpression: _check_general_expr, + EXPR.NPV_ProductExpression: _check_general_expr, + EXPR.DivisionExpression: _check_general_expr, + EXPR.NPV_DivisionExpression: _check_general_expr, + EXPR.PowExpression: _check_general_expr, + EXPR.NPV_PowExpression: _check_general_expr, + EXPR.NegationExpression: _check_general_expr, + EXPR.NPV_NegationExpression: _check_general_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: + 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 = ComponentSet() + self.mismatched_terms = ComponentSet() + + # Call parent walk_expression method + vals, const, _ = super().walk_expression(expr) + + # Return results + return vals, self.mismatched_terms, self.canceling_terms, const diff --git a/idaes/core/util/tests/test_model_diagnostics.py b/idaes/core/util/tests/test_model_diagnostics.py index ec5b671e34..f9bce71f55 100644 --- a/idaes/core/util/tests/test_model_diagnostics.py +++ b/idaes/core/util/tests/test_model_diagnostics.py @@ -13,52 +13,68 @@ """ 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.scaling import set_scaling_factor 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 +97,7 @@ _collect_model_statistics, check_parallel_jacobian, compute_ill_conditioning_certificate, + ConstraintTermAnalysisVisitor, ) from idaes.core.util.parameter_sweep import ( SequentialSweepRunner, @@ -90,6 +107,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 +1017,111 @@ 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: v4 + v5"] + assert cancellation == ["c2: v3 == v4 + v5", "c3: v6 == 10 + v3 - v4"] + assert constant == ["c1"] + + @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: v4 + v5 + +==================================================================================== +""" + + 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: v6 == 10 + v3 - v4 + +==================================================================================== +""" + + 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,7 +1274,8 @@ def test_collect_numerical_cautions(self, model): dt = DiagnosticsToolbox(model=model.b) cautions = dt._collect_numerical_cautions() - assert len(cautions) == 5 + + assert len(cautions) == 6 assert ( "Caution: 2 Variables with value close to their bounds (abs=1.0E-04, rel=1.0E-04)" in cautions @@ -1162,6 +1286,7 @@ def test_collect_numerical_cautions(self, model): assert ( "Caution: 1 Variable with extreme value (<1.0E-04 or >1.0E+04)" in cautions ) + assert "Caution: 1 Constraint with potential cancellation of terms" in cautions @pytest.mark.component def test_collect_numerical_cautions_jacobian(self): @@ -1344,9 +1469,9 @@ def test_report_numerical_issues_ok(self): No warnings found! ------------------------------------------------------------------------------------ -0 Cautions +1 Cautions - No cautions found! + Caution: 1 Constraint with potential cancellation of terms ------------------------------------------------------------------------------------ Suggested next steps: @@ -1425,12 +1550,13 @@ def test_report_numerical_issues(self, model): WARNING: 1 Variable at or outside bounds (tol=0.0E+00) ------------------------------------------------------------------------------------ -5 Cautions +6 Cautions Caution: 2 Variables with value close to their bounds (abs=1.0E-04, rel=1.0E-04) Caution: 2 Variables with value close to zero (tol=1.0E-08) Caution: 1 Variable with extreme value (<1.0E-04 or >1.0E+04) Caution: 1 Variable with None value + Caution: 1 Constraint with potential cancellation of terms Caution: 1 extreme Jacobian Entry (<1.0E-04 or >1.0E+04) ------------------------------------------------------------------------------------ @@ -4053,3 +4179,770 @@ 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() + sums = [i for i in visitor._generate_sum_combinations(terms)] + + expected = [ + 3, + 4, + 5, + 6, + 5, + 6, + 7, + 7, + 8, + 9, # 2-term sums + 6, + 7, + 8, + 8, + 9, + 10, # 3-term sums starting with 1 + 9, + 10, + 11, # 3-term sums starting with 2 + 12, # 3-term sum starting with 3 + 10, + 11, + 12, + 13, + 14, # 4-term sums + 15, # 5-term sum + ] + + assert sums == expected + + @pytest.mark.unit + def test_int(self): + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression(expr=7) + + assert vv == [7] + assert len(mm) == 0 + assert len(cc) == 0 + assert k + + @pytest.mark.unit + def test_float(self): + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression(expr=7.7) + + assert vv == [7.7] + assert len(mm) == 0 + assert len(cc) == 0 + assert k + + @pytest.mark.unit + def test_scalar_param(self): + m = ConcreteModel() + m.scalar_param = Param(initialize=1, mutable=True) + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( + expr=m.scalar_param + ) + + assert vv == [1] + assert len(mm) == 0 + assert len(cc) == 0 + assert k + + @pytest.mark.unit + def test_indexed_param(self): + m = ConcreteModel() + m.indexed_param = Param(["a", "b"], initialize=1, mutable=True) + + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( + expr=m.indexed_param["a"] + ) + assert vv == [1] + assert len(mm) == 0 + assert len(cc) == 0 + assert k + + @pytest.mark.unit + def test_scalar_var(self): + m = ConcreteModel() + m.scalar_var = Var(initialize=1) + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( + expr=m.scalar_var + ) + + assert vv == [1] + assert len(mm) == 0 + assert len(cc) == 0 + assert not k + + @pytest.mark.unit + def test_indexed_var(self): + m = ConcreteModel() + m.indexed_var = Var(["a", "b"], initialize=1) + + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( + expr=m.indexed_var["a"] + ) + assert vv == [1] + assert len(mm) == 0 + assert len(cc) == 0 + assert not k + + @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 = ConstraintTermAnalysisVisitor().walk_expression( + expr=m.scalar_var + ) + + assert vv == [1] + assert len(mm) == 0 + assert len(cc) == 0 + assert k + + @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 = ConstraintTermAnalysisVisitor().walk_expression( + expr=m.indexed_var["a"] + ) + assert vv == [1] + assert len(mm) == 0 + assert len(cc) == 0 + assert k + + @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 = 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 + + @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 = 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 + + # Fix v2, now constant + m.v2.fix() + + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) + + assert vv == [-1e-7, 1e7] + assert expr in mm + assert len(mm) == 1 + assert len(cc) == 0 + assert k + + @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 = 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 + + @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 = 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 + + @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 = 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 + + @pytest.mark.unit + def test_product_sum_expr(self, model): + expr = (model.v1 + model.v2) * (model.v3 + model.v4) + vv, mm, cc, k = 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 + + @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 = 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 + + @pytest.mark.unit + def test_division_expr(self, model): + expr = model.v1 / model.v2 / model.v3 + vv, mm, cc, k = 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 + + @pytest.mark.unit + def test_division_sum_expr(self, model): + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( + expr=(model.v1 + model.v2) / (model.v3 + model.v4) + ) + + assert vv == [pytest.approx((2 + 3) / (5.0000001 + 5), rel=1e-8)] + assert len(mm) == 0 + assert len(cc) == 0 + assert not k + + @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 = 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 + + @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 = ConstraintTermAnalysisVisitor().walk_expression( + expr=model.v1**model.v2 + ) + + assert vv == [8] + assert len(mm) == 0 + assert len(cc) == 0 + assert not k + + @pytest.mark.unit + def test_pow_sum_expr(self, model): + vv, mm, cc, k = 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 + + @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 = 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 + + @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 = ConstraintTermAnalysisVisitor().walk_expression( + expr=-func_model.v1 + ) + + assert vv == [-1] + assert len(mm) == 0 + assert len(cc) == 0 + assert not k + + @pytest.mark.unit + def test_negation_sum_expr(self, func_model): + expr = -(func_model.v1 - func_model.v2) + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) + + assert vv == [pytest.approx(-0.00001, rel=1e-8)] + assert len(mm) == 0 + assert expr in cc + assert len(cc) == 1 + assert not k + + # 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 = 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 + + @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 = 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 + + @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 = 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 + + @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 = 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 + + @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 = ConstraintTermAnalysisVisitor().walk_expression(expr=Z) + + assert vv == [pytest.approx(-2.1149075414767577, rel=1e-8)] + assert len(mm) == 0 + assert Z in cc + # We actually get two issues here, as one of the input expressions also cancels + assert len(cc) == 2 + assert not k + + # 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 = 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 + + @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 = 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 + + @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 = 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 + + @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 = 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 + + # Fix v1 and v2 to make first two terms constant + m.v1.fix() + m.v2.fix() + vv, mm, cc, k = 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 + + # Fix v3 to make all terms constant + m.v3.fix() + vv, mm, cc, k = 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 + + @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 = 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 + + @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 = 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 + + @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 = 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 + + # Check for tolerance of sum cancellation + m.v2.set_value(2.00022) + + vv, mm, cc, k = ConstraintTermAnalysisVisitor( + term_cancellation_tolerance=1e-4 + ).walk_expression(expr=expr) + + assert vv == [pytest.approx(3 * -0.00022, rel=1e-8)] + assert len(mm) == 0 + # This should pass as the difference is greater than tol + assert len(cc) == 0 + assert not k + + # Change tolerance so it should identify cancellation + vv, mm, cc, k = ConstraintTermAnalysisVisitor( + term_cancellation_tolerance=1e-3 + ).walk_expression(expr=expr) + + assert vv == [pytest.approx(3 * -0.00022, rel=1e-8)] + assert len(mm) == 0 + assert expr in cc + assert len(cc) == 1 + assert not k + + @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 = 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 + + @pytest.mark.unit + def test_canceling_equality_expr_issue(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 = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) + + assert vv == [0, pytest.approx(0, abs=1e-12)] + assert len(mm) == 0 + assert expr in cc + assert len(cc) == 1 + assert not k + + @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 = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) + + assert vv == [0, pytest.approx(0, abs=1e-12)] + assert len(mm) == 0 + assert expr in cc + assert len(cc) == 1 + assert not k + + # If we fix v3, this should be OK + m.v3.fix() + + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) + + assert vv == [0, pytest.approx(0, abs=1e-12)] + assert len(mm) == 0 + assert len(cc) == 0 + assert not k + + # 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 = 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 + + # Check to make sure simple linking constraints are not flagged as cancelling + @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 = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) + + assert vv == [-1, 1] + assert len(mm) == 0 + assert len(cc) == 0 + assert not k + + @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 = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) + + assert vv == [-1, 1] + assert len(mm) == 0 + assert len(cc) == 0 + assert not k + + @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]) + + # This is triggering an unexpected exception + # Still working to debug this + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( + expr=m.state[0].temperature + ) + + assert vv == [pytest.approx(235.0, rel=1e-8)] + assert len(mm) == 0 + assert len(cc) == 0 + assert not k