From 37e0faee6f708a507f3d273498a6efb28b4352bd Mon Sep 17 00:00:00 2001 From: Andrew Lee Date: Mon, 22 Jul 2024 16:11:52 -0400 Subject: [PATCH 01/20] Working on cancellation detection --- idaes/core/util/model_diagnostics.py | 167 +++- .../core/util/tests/test_model_diagnostics.py | 875 ++++++++++++++++++ 2 files changed, 1041 insertions(+), 1 deletion(-) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index f4488f4680..3d54c5a091 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 @@ -58,6 +59,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.param import ParamData from pyomo.repn.standard_repn import ( # pylint: disable=no-name-in-module generate_standard_repn, ) @@ -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 ( @@ -4122,3 +4127,163 @@ def _collect_model_statistics(model): ) return stats + + +class ConstraintTermAnalysisVisitor(EXPR.StreamBasedExpressionVisitor): + def __init__(self, term_mismatch_tolerance=1e6, term_cancellation_tolerance=1e-4): + super().__init__() + + self._mm_tol = log10(term_mismatch_tolerance) + self._sum_tol = term_cancellation_tolerance + + def _check_base_type(self, node): + # [value], [mismatched terms], [canceling terms], constant + if isinstance(node, VarData): + const = node.fixed + else: + const = True + return [value(node)], [], [], const + + def _check_sum_expression(self, node, child_data): + # For sums, collect all child values into a list + vals = [] + mismatch = [] + # We will check for cancellation in this node at the next level + # Pyomo is generally good at simplifying compound sums + cancelling = [] + const = True + # Collect data from child nodes + for d in child_data: + vals.append(sum(i for i in d[0])) + mismatch += d[1] + cancelling += d[2] + + # Expression is not constant if any child is not constant + if not d[3]: + 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 and log10(vl / vs) > self._mm_tol: + mismatch.append(str(node)) + + # [value], [mismatched terms], [canceling terms] + return vals, mismatch, cancelling, const + + def _check_other_expression(self, node, child_data): + # First, need to check to see if any child data is a list + # This indicates a sum expression + mismatch = [] + cancelling = [] + const = True + + for d in child_data: + # Collect all warnings about mismatched terms from child nodes + mismatch += d[1] + cancelling += d[2] + + # We will check for cancelling 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 + sums = self._sum_combinations(d[0]) + if any(i <= self._sum_tol * max(d[0]) for i in sums): + cancelling.append(str(node)) + + # Expression is not constant if any child is not constant + if not d[3]: + const = False + + # TODO: See if we can be more efficient about getting the value - might need node specific methods + # [value], [mismatched terms], [canceling terms], constant + return [value(node)], mismatch, cancelling, const + + def _check_equality_expression(self, node, child_data): + # (In)equality expressions are a special case of sum expressions + # We can start by just calling the method to check the sum expression + vals, mismatch, cancelling, const = self._check_sum_expression(node, child_data) + + # 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 + if str(node) not in mismatch and not any(d[3] for d in child_data): + # No constant terms, check for cancellation + # First, collect terms from both sides + t = [] + for d in child_data: + t += d[0] + + # Then check for cancellations + sums = self._sum_combinations(t) + if any(i <= self._sum_tol * max(t) for i in sums): + cancelling.append(str(node)) + + return vals, mismatch, cancelling, const + + def _sum_combinations(self, values_list): + sums = [] + for i in chain.from_iterable( + combinations(values_list, r) for r in range(2, len(values_list) + 1) + ): + sums.append(abs(sum(i))) + return sums + + node_type_method_map = { + EXPR.EqualityExpression: _check_equality_expression, + EXPR.InequalityExpression: _check_equality_expression, + EXPR.RangedExpression: _check_sum_expression, + EXPR.SumExpression: _check_sum_expression, + EXPR.NPV_SumExpression: _check_sum_expression, + EXPR.ProductExpression: _check_other_expression, + EXPR.MonomialTermExpression: _check_other_expression, + EXPR.NPV_ProductExpression: _check_other_expression, + EXPR.DivisionExpression: _check_other_expression, + EXPR.NPV_DivisionExpression: _check_other_expression, + EXPR.PowExpression: _check_other_expression, + EXPR.NPV_PowExpression: _check_other_expression, + EXPR.NegationExpression: _check_other_expression, + EXPR.NPV_NegationExpression: _check_other_expression, + EXPR.AbsExpression: _check_other_expression, + EXPR.NPV_AbsExpression: _check_other_expression, + EXPR.UnaryFunctionExpression: _check_other_expression, + EXPR.NPV_UnaryFunctionExpression: _check_other_expression, + 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): + # Return [node values], [mismatched terms], [cancelling terms], constant + # first check if the node is a leaf + nodetype = type(node) + + if nodetype in native_types: + return [node], [], [], True + + node_func = self.node_type_method_map.get(nodetype, None) + if node_func is not None: + return node_func(self, node, data) + + elif not node.is_expression_type(): + # this is a leaf, but not a native type + if nodetype is _PyomoUnit: + return [1], [], [], True + else: + # 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)}" + ) diff --git a/idaes/core/util/tests/test_model_diagnostics.py b/idaes/core/util/tests/test_model_diagnostics.py index 0a8bbfe4f3..7c1587900b 100644 --- a/idaes/core/util/tests/test_model_diagnostics.py +++ b/idaes/core/util/tests/test_model_diagnostics.py @@ -81,6 +81,7 @@ _collect_model_statistics, check_parallel_jacobian, compute_ill_conditioning_certificate, + ConstraintTermAnalysisVisitor, ) from idaes.core.util.parameter_sweep import ( SequentialSweepRunner, @@ -4053,3 +4054,877 @@ 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 ot generate sums of all combinations of terms + # excludes single term sums + terms = [1, 2, 3, 4, 5] + visitor = ConstraintTermAnalysisVisitor() + sums = visitor._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 mm == [] + assert cc == [] + assert k + + @pytest.mark.unit + def test_float(self): + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression(expr=7.7) + + assert vv == [7.7] + assert mm == [] + assert cc == [] + 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 mm == [] + assert cc == [] + 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 mm == [] + assert cc == [] + 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 mm == [] + assert cc == [] + 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 mm == [] + assert cc == [] + 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 mm == [] + assert cc == [] + 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 mm == [] + assert cc == [] + assert k + + @pytest.mark.unit + def test_equality_expr(self): + m = ConcreteModel() + m.v1 = Var(initialize=1e-7) + m.v2 = Var(initialize=1e7) + + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( + expr=m.v1 == m.v2 + ) + + assert vv == [1e-7, 1e7] + assert mm == ["v1 == v2"] + assert cc == [] + 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() + + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( + expr=m.v1 == m.v2 + ) + + assert vv == [1e-7, 1e7] + assert mm == ["v1 == v2"] + assert cc == [] + assert not k + + # Fix v2, now constant + m.v2.fix() + + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( + expr=m.v1 == m.v2 + ) + + assert vv == [1e-7, 1e7] + assert mm == ["v1 == v2"] + assert cc == [] + 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) + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( + expr=sum(m.v1[i] for i in m.v1) + ) + + assert vv == [1e-7, 1e7, 1e7] + assert mm == ["v1[a] + v1[b] + v1[c]"] + assert cc == [] + 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() + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( + expr=sum(m.v1[i] for i in m.v1) + ) + + assert vv == [1e-7, 1e7, 1e7] + assert mm == ["v1[a] + v1[b] + v1[c]"] + assert cc == [] + assert k + + # @pytest.mark.unit + # def test_inequality_expr(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=m.scalar_var <= m.indexed_var["a"] + # ) == [21, 22] + # + # @pytest.mark.unit + # def test_sum_expr(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=sum(m.indexed_var[i] for i in m.set) + # ) == [22, 23, 24] + # + # @pytest.mark.unit + # def test_additive_expr(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=m.scalar_var + m.indexed_var["a"] + m.scalar_param + # ) == [21, 22, 12] + # + # @pytest.mark.unit + # def test_additive_expr_w_negation(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=m.scalar_var + m.indexed_var["a"] - m.scalar_param + # ) == [21, 22, -12] + # + # @pytest.mark.unit + # def test_product_expr(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=m.scalar_var * m.indexed_var["a"] * m.scalar_param + # ) == [21 * 22 * 12] + # + # @pytest.mark.unit + # def test_product_sum_expr(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=(m.scalar_var + m.indexed_var["a"]) + # * (m.scalar_param + m.indexed_var["b"]) + # ) == [21 * 12, 21 * 23, 22 * 12, 22 * 23] + # + # @pytest.mark.unit + # def test_product_sum_expr_w_negation(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=(m.scalar_var - m.indexed_var["a"]) + # * (m.scalar_param - m.indexed_var["b"]) + # ) == [21 * 12, -21 * 23, -22 * 12, 22 * 23] + # + # @pytest.mark.unit + # def test_division_expr(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=m.scalar_var / m.indexed_var["a"] / m.scalar_param + # ) == [21 / 22 / 12] + # + # @pytest.mark.unit + # def test_division_sum_expr(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=(m.scalar_var + m.indexed_var["a"]) + # / (m.scalar_param + m.indexed_var["b"]) + # ) == [(21 + 22) / (12 + 23)] + # + # @pytest.mark.unit + # def test_division_sum_expr_w_negation(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=(m.scalar_var - m.indexed_var["a"]) + # / (m.scalar_param - m.indexed_var["b"]) + # ) == [(21 - 22) / (12 - 23)] + # + # @pytest.mark.unit + # def test_division_expr_error(self, m, caplog): + # caplog.set_level(logging.DEBUG, logger="idaes.core.util.scaling") + # sc.NominalValueExtractionVisitor().walk_expression(expr=1 / (m.scalar_var - 21)) + # + # expected = "Nominal value of 0 found in denominator of division expression. " + # "Assigning a value of 1. You should check you scaling factors and models to " + # "ensure there are no values of 0 that can appear in these functions." + # + # assert expected in caplog.text + # + # @pytest.mark.unit + # def test_pow_expr(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=m.scalar_var ** m.indexed_var["a"] + # ) == pytest.approx([21 ** 22], rel=1e-12) + # + # @pytest.mark.unit + # def test_pow_sum_expr(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=(m.scalar_var + m.indexed_var["a"]) + # ** (m.scalar_param + m.indexed_var["b"]) + # ) == [ + # pytest.approx((21 + 22) ** (12 + 23), rel=1e-12), + # ] + # + # @pytest.mark.unit + # def test_pow_sum_expr_w_negation(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=(m.scalar_var - m.indexed_var["a"]) + # ** (m.scalar_param - m.indexed_var["b"]) + # ) == [ + # pytest.approx(abs(21 - 22) ** (12 - 23), rel=1e-12), + # ] + # + # @pytest.mark.unit + # def test_negation_expr(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=-m.scalar_var + # ) == [-21] + # + # @pytest.mark.unit + # def test_negation_sum_expr(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=-(m.scalar_var + m.indexed_var["a"]) + # ) == [-21, -22] + # + # @pytest.mark.unit + # def test_log_expr(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.log(m.scalar_var) + # ) == [pytest.approx(math.log(21), rel=1e-12)] + # + # @pytest.mark.unit + # def test_log_expr_error(self, m): + # with pytest.raises( + # ValueError, + # match="Evaluation error occurred when getting nominal value in log expression " + # "with input 0.0. You should check you scaling factors and model to " + # "address any numerical issues or scale this constraint manually.", + # ): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.log(m.scalar_var - 21) + # ) + # + # @pytest.mark.unit + # def test_log_sum_expr(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.log(m.scalar_var + m.indexed_var["a"]) + # ) == [pytest.approx(math.log(21 + 22), rel=1e-12)] + # + # @pytest.mark.unit + # def test_log_sum_expr_w_negation(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.log(-m.scalar_var + m.indexed_var["a"]) + # ) == [pytest.approx(math.log(-21 + 22), rel=1e-12)] + # + # @pytest.mark.unit + # def test_log10_expr(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.log10(m.scalar_var) + # ) == [pytest.approx(math.log10(21), rel=1e-12)] + # + # @pytest.mark.unit + # def test_log10_sum_expr(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.log10(m.scalar_var + m.indexed_var["a"]) + # ) == [pytest.approx(math.log10(21 + 22), rel=1e-12)] + # + # @pytest.mark.unit + # def test_log10_sum_expr_w_negation(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.log10(-m.scalar_var + m.indexed_var["a"]) + # ) == [pytest.approx(math.log10(-21 + 22), rel=1e-12)] + # + # @pytest.mark.unit + # def test_log10_expr_error(self, m): + # with pytest.raises( + # ValueError, + # match="Evaluation error occurred when getting nominal value in log10 expression " + # "with input 0.0. You should check you scaling factors and model to " + # "address any numerical issues or scale this constraint manually.", + # ): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.log10(m.scalar_var - 21) + # ) + # + # @pytest.mark.unit + # def test_sqrt_expr(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.sqrt(m.scalar_var) + # ) == [pytest.approx(21 ** 0.5, rel=1e-12)] + # + # @pytest.mark.unit + # def test_sqrt_sum_expr(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.sqrt(m.scalar_var + m.indexed_var["a"]) + # ) == [pytest.approx((21 + 22) ** 0.5, rel=1e-12)] + # + # @pytest.mark.unit + # def test_sqrt_sum_expr_w_negation(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.sqrt(-m.scalar_var + m.indexed_var["a"]) + # ) == [pytest.approx((-21 + 22) ** 0.5, rel=1e-12)] + # + # @pytest.mark.unit + # def test_sqrt_expr_error(self, m): + # with pytest.raises( + # ValueError, + # match="Evaluation error occurred when getting nominal value in sqrt expression " + # "with input -21.0. You should check you scaling factors and model to " + # "address any numerical issues or scale this constraint manually.", + # ): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.sqrt(-m.scalar_var) + # ) + # + # @pytest.mark.unit + # def test_sin_expr(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.sin(m.scalar_var) + # ) == [pytest.approx(math.sin(21), rel=1e-12)] + # + # @pytest.mark.unit + # def test_sin_sum_expr(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.sin(m.scalar_var + m.indexed_var["a"]) + # ) == [pytest.approx(math.sin(21 + 22), rel=1e-12)] + # + # @pytest.mark.unit + # def test_sin_sum_expr_w_negation(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.sin(m.scalar_var - m.indexed_var["a"]) + # ) == [pytest.approx(math.sin(21 - 22), rel=1e-12)] + # + # @pytest.mark.unit + # def test_cos_expr(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.cos(m.scalar_var) + # ) == [pytest.approx(math.cos(21), rel=1e-12)] + # + # @pytest.mark.unit + # def test_cos_sum_expr(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.cos(m.scalar_var + m.indexed_var["a"]) + # ) == [pytest.approx(math.cos(21 + 22), rel=1e-12)] + # + # @pytest.mark.unit + # def test_cos_sum_expr_w_negation(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.cos(m.scalar_var - m.indexed_var["a"]) + # ) == [pytest.approx(math.cos(21 - 22), rel=1e-12)] + # + # @pytest.mark.unit + # def test_tan_expr(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.tan(m.scalar_var) + # ) == [pytest.approx(math.tan(21), rel=1e-12)] + # + # @pytest.mark.unit + # def test_tan_sum_expr(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.tan(m.scalar_var + m.indexed_var["a"]) + # ) == [pytest.approx(math.tan(21 + 22), rel=1e-12)] + # + # @pytest.mark.unit + # def test_tan_sum_expr_w_negation(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.tan(m.scalar_var - m.indexed_var["a"]) + # ) == [pytest.approx(math.tan(21 - 22), rel=1e-12)] + # + # @pytest.mark.unit + # def test_sinh_expr(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.sinh(m.scalar_var) + # ) == [pytest.approx(math.sinh(21), rel=1e-12)] + # + # @pytest.mark.unit + # def test_sinh_sum_expr(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.sinh(m.scalar_var + m.indexed_var["a"]) + # ) == [pytest.approx(math.sinh(21 + 22), rel=1e-12)] + # + # @pytest.mark.unit + # def test_sinh_sum_expr_w_negation(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.sinh(m.scalar_var - m.indexed_var["a"]) + # ) == [pytest.approx(math.sinh(21 - 22), rel=1e-12)] + # + # @pytest.mark.unit + # def test_cosh_expr(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.cosh(m.scalar_var) + # ) == [pytest.approx(math.cosh(21), rel=1e-12)] + # + # @pytest.mark.unit + # def test_cosh_sum_expr(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.cosh(m.scalar_var + m.indexed_var["a"]) + # ) == [pytest.approx(math.cosh(21 + 22), rel=1e-12)] + # + # @pytest.mark.unit + # def test_cosh_sum_expr_w_negation(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.cosh(m.scalar_var - m.indexed_var["a"]) + # ) == [pytest.approx(math.cosh(21 - 22), rel=1e-12)] + # + # @pytest.mark.unit + # def test_tanh_expr(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.tanh(m.scalar_var) + # ) == [pytest.approx(math.tanh(21), rel=1e-12)] + # + # @pytest.mark.unit + # def test_tanh_sum_expr(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.tanh(m.scalar_var + m.indexed_var["a"]) + # ) == [pytest.approx(math.tanh(21 + 22), rel=1e-12)] + # + # @pytest.mark.unit + # def test_tanh_sum_expr_w_negation(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.tanh(m.scalar_var - m.indexed_var["a"]) + # ) == [pytest.approx(math.tanh(21 - 22), rel=1e-12)] + # + # @pytest.mark.unit + # def test_asin_expr(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression(expr=pyo.asin(1)) == [ + # pytest.approx(math.asin(1), rel=1e-12) + # ] + # + # @pytest.mark.unit + # def test_asin_sum_expr(self, m): + # sc.set_scaling_factor(m.scalar_param, 2) + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.asin(0.5 + m.scalar_param) + # ) == [pytest.approx(math.asin(1), rel=1e-12)] + # + # @pytest.mark.unit + # def test_asin_sum_expr_negation(self, m): + # sc.set_scaling_factor(m.scalar_param, 2) + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.asin(0.5 - m.scalar_param) + # ) == [pytest.approx(math.asin(0), rel=1e-12)] + # + # @pytest.mark.unit + # def test_acos_expr(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.acos(m.scalar_param) + # ) == [pytest.approx(math.acos(0.5), rel=1e-12)] + # + # @pytest.mark.unit + # def test_acos_sum_expr(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.acos(0.5 + m.scalar_param) + # ) == [pytest.approx(math.acos(1), rel=1e-12)] + # + # @pytest.mark.unit + # def test_acos_sum_expr_w_negation(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.acos(0.5 - m.scalar_param) + # ) == [pytest.approx(math.acos(0), rel=1e-12)] + # + # @pytest.mark.unit + # def test_asinh_expr(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.asinh(m.scalar_var) + # ) == [pytest.approx(math.asinh(21), rel=1e-12)] + # + # @pytest.mark.unit + # def test_asinh_sum_expr(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.asinh(m.scalar_var + m.indexed_var["a"]) + # ) == [pytest.approx(math.asinh(21 + 22), rel=1e-12)] + # + # @pytest.mark.unit + # def test_asinh_sum_expr_w_negation(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.asinh(m.scalar_var - m.indexed_var["a"]) + # ) == [pytest.approx(math.asinh(21 - 22), rel=1e-12)] + # + # @pytest.mark.unit + # def test_acosh_expr(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.acosh(m.scalar_var) + # ) == [pytest.approx(math.acosh(21), rel=1e-12)] + # + # @pytest.mark.unit + # def test_acosh_sum_expr(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.acosh(m.scalar_var + m.indexed_var["a"]) + # ) == [pytest.approx(math.acosh(21 + 22), rel=1e-12)] + # + # @pytest.mark.unit + # def test_acosh_sum_expr_w_negation(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.acosh(-m.scalar_var + m.indexed_var["a"]) + # ) == [pytest.approx(math.acosh(-21 + 22), rel=1e-12)] + # + # @pytest.mark.unit + # def test_atanh_expr(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.atanh(m.scalar_param) + # ) == [pytest.approx(math.atanh(0.5), rel=1e-12)] + # + # @pytest.mark.unit + # def test_atanh_sum_expr(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.atanh(0.4 + m.scalar_param) + # ) == [pytest.approx(math.atanh(0.9), rel=1e-12)] + # + # @pytest.mark.unit + # def test_atanh_sum_expr_w_negation(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.atanh(0.4 - m.scalar_param) + # ) == [pytest.approx(math.atanh(-0.1), rel=1e-12)] + # + # @pytest.mark.unit + # def test_exp_expr(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.exp(m.scalar_param) + # ) == [pytest.approx(math.exp(0.5), rel=1e-12)] + # + # @pytest.mark.unit + # def test_exp_sum_expr(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.exp(0.4 + m.scalar_param) + # ) == [pytest.approx(math.exp(0.9), rel=1e-12)] + # + # @pytest.mark.unit + # def test_exp_sum_expr_w_negation(self, m): + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=pyo.exp(-0.4 + m.scalar_param) + # ) == [pytest.approx(math.exp(0.1), rel=1e-12)] + # + # @pytest.mark.unit + # def test_expr_if(self, m): + # m.exprif = pyo.Expr_if( + # IF=m.scalar_param, + # THEN=m.indexed_var["a"], + # ELSE=m.indexed_var["b"] + m.indexed_var["c"], + # ) + # + # assert sc.NominalValueExtractionVisitor().walk_expression(expr=m.exprif) == [ + # 22, + # 23, + # 24, + # ] + # + # @pytest.mark.unit + # def test_expr_if_w_negation(self, m): + # m.exprif = pyo.Expr_if( + # IF=m.scalar_param, + # THEN=m.indexed_var["a"], + # ELSE=m.indexed_var["b"] - m.indexed_var["c"], + # ) + # + # assert sc.NominalValueExtractionVisitor().walk_expression(expr=m.exprif) == [ + # 22, + # 23, + # -24, + # ] + # + # @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 = pyo.ConcreteModel() + # m.a = pyo.Var(initialize=1) + # m.b = pyo.Var(initialize=1) + # + # sc.set_scaling_factor(m.a, 1 / 2) + # sc.set_scaling_factor(m.b, 1 / 4) + # + # m.expr_write = CubicThermoExpressions(m) + # Z = m.expr_write.z_liq(eos=CubicEoS.PR, A=m.a, B=m.b) + # + # expected_mag = -9.489811292072448 + # assert sc.NominalValueExtractionVisitor().walk_expression(expr=Z) == [ + # pytest.approx(expected_mag, rel=1e-8) + # ] + # + # # Check that model state did not change + # assert pyo.value(m.a) == 1 + # assert pyo.value(m.b) == 1 + # assert pyo.value(Z) == pytest.approx(-2.1149075414767577, rel=1e-8) + # + # # Now, change the actual state to the expected magnitudes and confirm result + # m.a.set_value(2) + # m.b.set_value(4) + # assert pyo.value(Z) == pytest.approx(expected_mag, rel=1e-8) + # + # @pytest.mark.unit + # def test_Expression(self, m): + # m.expression = pyo.Expression( + # expr=m.scalar_param ** (sum(m.indexed_var[i] for i in m.set)) + # ) + # + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=m.expression + # ) == [0.5 ** (22 + 23 + 24)] + # + # @pytest.mark.unit + # def test_constraint(self, m): + # m.constraint = pyo.Constraint(expr=m.scalar_var == m.expression) + # + # assert sc.NominalValueExtractionVisitor().walk_expression( + # expr=m.constraint.expr + # ) == [21, 0.5 ** (22 + 23 + 24)] + + @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) + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( + expr=m.v2 == sum(m.v1[i] for i in m.v1) + ) + + assert vv == [1e-7, 3e7] + assert mm == ["v2 == v1[a] + v1[b] + v1[c]"] + assert cc == [] + 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) + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( + expr=6 * m.v2 == 8 * sum(m.v1[i] for i in m.v1) + ) + + assert vv == [6e-7, 2.4e8] + assert mm == ["6*v2 == 8*(v1[a] + v1[b] + v1[c])"] + assert cc == [] + assert not 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) + + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( + expr=6 * m.v2 == 8 * sum(m.v1[i] for i in m.v1) + m.v3 + ) + + assert vv == [6e-7, pytest.approx(8 * (2e7 + 1e-7) + 1000, rel=1e-8)] + assert mm == [ + "v1[a] + v1[b] + v1[c]", + "6*v2 == 8*(v1[a] + v1[b] + v1[c]) + v3", + ] + assert cc == [] + assert not k + + @pytest.mark.unit + def test_cancelling_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 mm == [] + # We do not check cancellation at the sum, so this should be empty + assert cc == [] + assert not k + + @pytest.mark.unit + def test_expr_w_cancelling_sum(self): + m = ConcreteModel() + m.v1 = Var(initialize=2) + m.v2 = Var(initialize=2) + m.v3 = Var(initialize=3) + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( + expr=m.v3 * (m.v1 - m.v2) + ) + + assert vv == [0] + assert mm == [] + # We should get a warning about cancelling sums here + assert cc == ["v3*(v1 - v2)"] + 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=m.v3 * (m.v1 - m.v2)) + + assert vv == [pytest.approx(3 * -0.00022, rel=1e-8)] + assert mm == [] + # This should pass as the difference is greater than tol + assert cc == [] + assert not k + + # Change tolerance so it should identify cancellation + vv, mm, cc, k = ConstraintTermAnalysisVisitor( + term_cancellation_tolerance=1e-3 + ).walk_expression(expr=m.v3 * (m.v1 - m.v2)) + + assert vv == [pytest.approx(3 * -0.00022, rel=1e-8)] + assert mm == [] + assert cc == ["v3*(v1 - v2)"] + assert not k + + @pytest.mark.unit + def test_cancelling_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 mm == [] + assert cc == [] + assert not k + + @pytest.mark.unit + def test_cancelling_equality_expr_issue(self): + m = ConcreteModel() + m.v1 = Var(initialize=1e7) + m.v2 = Var(initialize=1e7) + m.v3 = Var(initialize=0) + + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( + expr=m.v3 == m.v1 - m.v2 + ) + + assert vv == [0, pytest.approx(0, abs=1e-12)] + assert mm == [] + assert cc == ["v3 == v1 - v2"] + assert not k + + @pytest.mark.unit + def test_cancelling_equality_expr_compound(self): + m = ConcreteModel() + m.v1 = Var(["a", "b"], initialize=5e6) + m.v2 = Var(initialize=1e7) + m.v3 = Var(initialize=0) + + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( + expr=m.v3 == sum(v for v in m.v1.values()) - m.v2 + ) + + assert vv == [0, pytest.approx(0, abs=1e-12)] + assert mm == [] + assert cc == ["v3 == v1[a] + v1[b] - v2"] + assert not k + + # If we fix v3, this should be OK + m.v3.fix() + + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( + expr=m.v3 == sum(v for v in m.v1.values()) - m.v2 + ) + + assert vv == [0, pytest.approx(0, abs=1e-12)] + assert mm == [] + assert cc == [] + assert not k From 4d48741acc7e8a6d2d875591044347e2c2e9f30e Mon Sep 17 00:00:00 2001 From: Andrew Lee Date: Mon, 22 Jul 2024 16:34:19 -0400 Subject: [PATCH 02/20] Fixing typo --- idaes/core/util/tests/test_model_diagnostics.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/idaes/core/util/tests/test_model_diagnostics.py b/idaes/core/util/tests/test_model_diagnostics.py index 4f504f3d8f..7ebe1e400c 100644 --- a/idaes/core/util/tests/test_model_diagnostics.py +++ b/idaes/core/util/tests/test_model_diagnostics.py @@ -4059,7 +4059,7 @@ def test_output(self, model): class TestConstraintTermAnalysisVisitor: @pytest.mark.unit def test_sum_combinations(self): - # Check method ot generate sums of all combinations of terms + # Check method to generate sums of all combinations of terms # excludes single term sums terms = [1, 2, 3, 4, 5] visitor = ConstraintTermAnalysisVisitor() From 33a4435bd4225b1e3d75f2ac50cb95876a493d0c Mon Sep 17 00:00:00 2001 From: Andrew Lee Date: Tue, 6 Aug 2024 14:30:17 -0400 Subject: [PATCH 03/20] Improving efficiency --- idaes/core/util/model_diagnostics.py | 140 +++++++++- .../core/util/tests/test_model_diagnostics.py | 239 ++++++++++++------ 2 files changed, 289 insertions(+), 90 deletions(-) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index 528f832468..bff5f8e966 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -59,7 +59,6 @@ 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.param import ParamData from pyomo.repn.standard_repn import ( # pylint: disable=no-name-in-module generate_standard_repn, ) @@ -4130,7 +4129,30 @@ def _collect_model_statistics(model): class ConstraintTermAnalysisVisitor(EXPR.StreamBasedExpressionVisitor): - def __init__(self, term_mismatch_tolerance=1e6, term_cancellation_tolerance=1e-4): + """ + 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__() self._mm_tol = log10(term_mismatch_tolerance) @@ -4144,7 +4166,11 @@ def _check_base_type(self, node): const = True return [value(node)], [], [], const + def _get_value_for_sum_subexpression(self, child_data): + return sum(i for i in child_data[0]) + def _check_sum_expression(self, node, child_data): + # Sum expressions need special handling # For sums, collect all child values into a list vals = [] mismatch = [] @@ -4154,7 +4180,7 @@ def _check_sum_expression(self, node, child_data): const = True # Collect data from child nodes for d in child_data: - vals.append(sum(i for i in d[0])) + vals.append(self._get_value_for_sum_subexpression(d)) mismatch += d[1] cancelling += d[2] @@ -4170,10 +4196,94 @@ def _check_sum_expression(self, node, child_data): if vl != vs and log10(vl / vs) > self._mm_tol: mismatch.append(str(node)) - # [value], [mismatched terms], [canceling terms] return vals, mismatch, cancelling, const + def _check_product(self, node, child_data): + mismatch, cancelling, const = self._perform_checks(node, child_data) + + val = self._get_value_for_sum_subexpression( + child_data[0] + ) * self._get_value_for_sum_subexpression(child_data[1]) + + return [val], mismatch, cancelling, const + + def _check_division(self, node, child_data): + mismatch, cancelling, const = self._perform_checks(node, child_data) + + numerator = self._get_value_for_sum_subexpression(child_data[0]) + denominator = self._get_value_for_sum_subexpression(child_data[1]) + # TODO: should we look at mismatched magnitudes in num and denom? + if denominator == 0: + raise ValueError( + f"Error in ConstraintTermAnalysisVisitor: found division with denominator of 0 " + f"({str(node)})." + ) + + return [numerator / denominator], mismatch, cancelling, const + + def _check_power(self, node, child_data): + mismatch, cancelling, const = self._perform_checks(node, child_data) + + base = self._get_value_for_sum_subexpression(child_data[0]) + exponent = self._get_value_for_sum_subexpression(child_data[1]) + + return [base**exponent], mismatch, cancelling, const + + # def _get_nominal_value_single_child(self, node, child_nominal_values): + # assert len(child_nominal_values) == 1 + # return child_nominal_values[0] + # + # def _get_nominal_value_abs(self, node, child_nominal_values): + # assert len(child_nominal_values) == 1 + # return [abs(i) for i in child_nominal_values[0]] + # + # def _get_nominal_value_negation(self, node, child_nominal_values): + # assert len(child_nominal_values) == 1 + # return [-i for i in child_nominal_values[0]] + # + # def _get_nominal_value_for_unary_function(self, node, child_nominal_values): + # assert len(child_nominal_values) == 1 + # func_name = node.getname() + # # TODO: Some of these need the absolute value of the nominal value (e.g. sqrt) + # func_nominal = self._get_nominal_value_for_sum_subexpression( + # child_nominal_values[0] + # ) + # func = getattr(math, func_name) + # try: + # return [func(func_nominal)] + # except ValueError: + # raise ValueError( + # f"Evaluation error occurred when getting nominal value in {func_name} " + # f"expression with input {func_nominal}. You should check you scaling factors " + # f"and model to address any numerical issues or scale this constraint manually." + # ) + # + # def _get_nominal_value_expr_if(self, node, child_nominal_values): + # assert len(child_nominal_values) == 3 + # return child_nominal_values[1] + child_nominal_values[2] + # + # def _get_nominal_value_external_function(self, node, child_nominal_values): + # # First, need to get expected magnitudes of input terms, which may be sub-expressions + # input_mag = [ + # self._get_nominal_value_for_sum_subexpression(i) + # for i in child_nominal_values + # ] + # + # # Next, create a copy of the external function with expected magnitudes as inputs + # newfunc = node.create_node_with_local_data(input_mag) + # + # # Evaluate new function and return the absolute value + # return [pyo.value(newfunc)] + def _check_other_expression(self, node, child_data): + mismatch, cancelling, const = self._perform_checks(node, child_data) + + # TODO: See if we can be more efficient about getting the value - might need node specific methods + # [value], [mismatched terms], [canceling terms], constant + return [value(node)], mismatch, cancelling, const + + 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 mismatch = [] @@ -4196,9 +4306,8 @@ def _check_other_expression(self, node, child_data): if not d[3]: const = False - # TODO: See if we can be more efficient about getting the value - might need node specific methods - # [value], [mismatched terms], [canceling terms], constant - return [value(node)], mismatch, cancelling, const + # Return any problematic terms found + return mismatch, cancelling, const def _check_equality_expression(self, node, child_data): # (In)equality expressions are a special case of sum expressions @@ -4236,13 +4345,13 @@ def _sum_combinations(self, values_list): EXPR.RangedExpression: _check_sum_expression, EXPR.SumExpression: _check_sum_expression, EXPR.NPV_SumExpression: _check_sum_expression, - EXPR.ProductExpression: _check_other_expression, - EXPR.MonomialTermExpression: _check_other_expression, - EXPR.NPV_ProductExpression: _check_other_expression, - EXPR.DivisionExpression: _check_other_expression, - EXPR.NPV_DivisionExpression: _check_other_expression, - EXPR.PowExpression: _check_other_expression, - EXPR.NPV_PowExpression: _check_other_expression, + EXPR.ProductExpression: _check_product, + EXPR.MonomialTermExpression: _check_product, + EXPR.NPV_ProductExpression: _check_product, + EXPR.DivisionExpression: _check_division, + EXPR.NPV_DivisionExpression: _check_division, + EXPR.PowExpression: _check_power, + EXPR.NPV_PowExpression: _check_power, EXPR.NegationExpression: _check_other_expression, EXPR.NPV_NegationExpression: _check_other_expression, EXPR.AbsExpression: _check_other_expression, @@ -4256,6 +4365,9 @@ def _sum_combinations(self, values_list): } def exitNode(self, node, data): + """ + Method to call when exiting node to check for potential issues. + """ # Return [node values], [mismatched terms], [cancelling terms], constant # first check if the node is a leaf nodetype = type(node) diff --git a/idaes/core/util/tests/test_model_diagnostics.py b/idaes/core/util/tests/test_model_diagnostics.py index 7ebe1e400c..1fe0c37e97 100644 --- a/idaes/core/util/tests/test_model_diagnostics.py +++ b/idaes/core/util/tests/test_model_diagnostics.py @@ -4291,82 +4291,169 @@ def test_sum_expr_constant(self): # assert sc.NominalValueExtractionVisitor().walk_expression( # expr=m.scalar_var + m.indexed_var["a"] - m.scalar_param # ) == [21, 22, -12] - # - # @pytest.mark.unit - # def test_product_expr(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=m.scalar_var * m.indexed_var["a"] * m.scalar_param - # ) == [21 * 22 * 12] - # - # @pytest.mark.unit - # def test_product_sum_expr(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=(m.scalar_var + m.indexed_var["a"]) - # * (m.scalar_param + m.indexed_var["b"]) - # ) == [21 * 12, 21 * 23, 22 * 12, 22 * 23] - # - # @pytest.mark.unit - # def test_product_sum_expr_w_negation(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=(m.scalar_var - m.indexed_var["a"]) - # * (m.scalar_param - m.indexed_var["b"]) - # ) == [21 * 12, -21 * 23, -22 * 12, 22 * 23] - # - # @pytest.mark.unit - # def test_division_expr(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=m.scalar_var / m.indexed_var["a"] / m.scalar_param - # ) == [21 / 22 / 12] - # - # @pytest.mark.unit - # def test_division_sum_expr(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=(m.scalar_var + m.indexed_var["a"]) - # / (m.scalar_param + m.indexed_var["b"]) - # ) == [(21 + 22) / (12 + 23)] - # - # @pytest.mark.unit - # def test_division_sum_expr_w_negation(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=(m.scalar_var - m.indexed_var["a"]) - # / (m.scalar_param - m.indexed_var["b"]) - # ) == [(21 - 22) / (12 - 23)] - # - # @pytest.mark.unit - # def test_division_expr_error(self, m, caplog): - # caplog.set_level(logging.DEBUG, logger="idaes.core.util.scaling") - # sc.NominalValueExtractionVisitor().walk_expression(expr=1 / (m.scalar_var - 21)) - # - # expected = "Nominal value of 0 found in denominator of division expression. " - # "Assigning a value of 1. You should check you scaling factors and models to " - # "ensure there are no values of 0 that can appear in these functions." - # - # assert expected in caplog.text - # - # @pytest.mark.unit - # def test_pow_expr(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=m.scalar_var ** m.indexed_var["a"] - # ) == pytest.approx([21 ** 22], rel=1e-12) - # - # @pytest.mark.unit - # def test_pow_sum_expr(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=(m.scalar_var + m.indexed_var["a"]) - # ** (m.scalar_param + m.indexed_var["b"]) - # ) == [ - # pytest.approx((21 + 22) ** (12 + 23), rel=1e-12), - # ] - # - # @pytest.mark.unit - # def test_pow_sum_expr_w_negation(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=(m.scalar_var - m.indexed_var["a"]) - # ** (m.scalar_param - m.indexed_var["b"]) - # ) == [ - # pytest.approx(abs(21 - 22) ** (12 - 23), rel=1e-12), - # ] - # + + @pytest.mark.unit + def test_product_expr(self): + m = ConcreteModel() + m.v1 = Var(initialize=2) + m.v2 = Var(initialize=3) + m.v3 = Var(initialize=5) + + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( + expr=m.v1 * m.v2 * m.v3 + ) + + assert vv == [30] + assert mm == [] + assert cc == [] + assert not k + + @pytest.mark.unit + def test_product_sum_expr(self): + m = ConcreteModel() + m.v1 = Var(initialize=2) + m.v2 = Var(initialize=3) + m.v3 = Var(initialize=5) + m.v4 = Var(initialize=6) + + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( + expr=(m.v1 + m.v2) * (m.v3 + m.v4) + ) + + assert vv == [(2 + 3) * (5 + 6)] + assert mm == [] + assert cc == [] + assert not k + + @pytest.mark.unit + def test_product_sum_expr_w_negation(self): + m = ConcreteModel() + m.v1 = Var(initialize=2) + m.v2 = Var(initialize=3) + m.v3 = Var(initialize=5) + m.v4 = Var(initialize=5) + + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( + expr=(m.v1 + m.v2) * (m.v3 - m.v4) + ) + + assert vv == [0] + assert mm == [] + assert cc == ["(v1 + v2)*(v3 - v4)"] + assert not k + + @pytest.mark.unit + def test_division_expr(self): + m = ConcreteModel() + m.v1 = Var(initialize=2) + m.v2 = Var(initialize=3) + m.v3 = Var(initialize=5) + + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( + expr=m.v1 / m.v2 / m.v3 + ) + + assert vv == [pytest.approx(2 / 3 / 5, rel=1e-8)] + assert mm == [] + assert cc == [] + assert not k + + @pytest.mark.unit + def test_division_sum_expr(self): + m = ConcreteModel() + m.v1 = Var(initialize=2) + m.v2 = Var(initialize=3) + m.v3 = Var(initialize=5) + m.v4 = Var(initialize=6) + + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( + expr=(m.v1 + m.v2) / (m.v3 + m.v4) + ) + + assert vv == [pytest.approx((2 + 3) / (5 + 6), rel=1e-8)] + assert mm == [] + assert cc == [] + assert not k + + @pytest.mark.unit + def test_division_sum_expr_w_negation(self): + m = ConcreteModel() + m.v1 = Var(initialize=2) + m.v2 = Var(initialize=3) + m.v3 = Var(initialize=5.0000001) + m.v4 = Var(initialize=5) + + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( + expr=(m.v1 + m.v2) / (m.v3 - m.v4) + ) + + assert vv == [pytest.approx((2 + 3) / (0.0000001), rel=1e-8)] + assert mm == [] + assert cc == ["(v1 + v2)/(v3 - v4)"] + 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( + ValueError, + 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): + m = ConcreteModel() + m.v1 = Var(initialize=2) + m.v2 = Var(initialize=3) + + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression(expr=m.v1**m.v2) + + assert vv == [8] + assert mm == [] + assert cc == [] + assert not k + + @pytest.mark.unit + def test_pow_sum_expr(self): + m = ConcreteModel() + m.v1 = Var(initialize=2) + m.v2 = Var(initialize=3) + m.v3 = Var(initialize=5) + m.v4 = Var(initialize=6) + + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( + expr=(m.v1 + m.v2) ** (m.v3 + m.v4) + ) + + assert vv == [5**11] + assert mm == [] + assert cc == [] + assert not k + + @pytest.mark.unit + def test_pow_sum_expr_w_negation(self): + m = ConcreteModel() + m.v1 = Var(initialize=2) + m.v2 = Var(initialize=3) + m.v3 = Var(initialize=5.0000001) + m.v4 = Var(initialize=5) + + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( + expr=(m.v1 + m.v2) ** (m.v3 - m.v4) + ) + + assert vv == [pytest.approx((2 + 3) ** (0.0000001), rel=1e-8)] + assert mm == [] + assert cc == ["(v1 + v2)**(v3 - v4)"] + assert not k + # @pytest.mark.unit # def test_negation_expr(self, m): # assert sc.NominalValueExtractionVisitor().walk_expression( From 747eb7df272febdb4f03b1df7b2430924cf5a802 Mon Sep 17 00:00:00 2001 From: Andrew Lee Date: Tue, 6 Aug 2024 14:58:05 -0400 Subject: [PATCH 04/20] More efficiency and testing --- idaes/core/util/model_diagnostics.py | 70 ++--- .../core/util/tests/test_model_diagnostics.py | 271 ++++++++++-------- 2 files changed, 187 insertions(+), 154 deletions(-) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index bff5f8e966..ab3e1252ec 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -17,6 +17,7 @@ __author__ = "Alexander Dowling, Douglas Allan, Andrew Lee, Robby Parker, Ben Knueven" +import math from operator import itemgetter import sys from inspect import signature @@ -4229,35 +4230,34 @@ def _check_power(self, node, child_data): return [base**exponent], mismatch, cancelling, const - # def _get_nominal_value_single_child(self, node, child_nominal_values): - # assert len(child_nominal_values) == 1 - # return child_nominal_values[0] - # - # def _get_nominal_value_abs(self, node, child_nominal_values): - # assert len(child_nominal_values) == 1 - # return [abs(i) for i in child_nominal_values[0]] - # - # def _get_nominal_value_negation(self, node, child_nominal_values): - # assert len(child_nominal_values) == 1 - # return [-i for i in child_nominal_values[0]] - # - # def _get_nominal_value_for_unary_function(self, node, child_nominal_values): - # assert len(child_nominal_values) == 1 - # func_name = node.getname() - # # TODO: Some of these need the absolute value of the nominal value (e.g. sqrt) - # func_nominal = self._get_nominal_value_for_sum_subexpression( - # child_nominal_values[0] - # ) - # func = getattr(math, func_name) - # try: - # return [func(func_nominal)] - # except ValueError: - # raise ValueError( - # f"Evaluation error occurred when getting nominal value in {func_name} " - # f"expression with input {func_nominal}. You should check you scaling factors " - # f"and model to address any numerical issues or scale this constraint manually." - # ) - # + def _check_negation(self, node, child_data): + mismatch, cancelling, const = self._perform_checks(node, child_data) + val = -self._get_value_for_sum_subexpression(child_data[0]) + + return [val], mismatch, cancelling, const + + def _check_abs(self, node, child_data): + mismatch, cancelling, const = self._perform_checks(node, child_data) + val = abs(self._get_value_for_sum_subexpression(child_data[0])) + + return [val], mismatch, cancelling, const + + def _check_unary_function(self, node, child_data): + mismatch, cancelling, const = self._perform_checks(node, child_data) + + func_name = node.getname() + func = getattr(math, func_name) + func_val = self._get_value_for_sum_subexpression(child_data[0]) + + try: + val = func(func_val) + except ValueError: + raise ValueError( + f"Error in ConstraintTermAnalysisVisitor: error evaluating {str(node)}." + ) + + return [val], mismatch, cancelling, const + # def _get_nominal_value_expr_if(self, node, child_nominal_values): # assert len(child_nominal_values) == 3 # return child_nominal_values[1] + child_nominal_values[2] @@ -4352,12 +4352,12 @@ def _sum_combinations(self, values_list): EXPR.NPV_DivisionExpression: _check_division, EXPR.PowExpression: _check_power, EXPR.NPV_PowExpression: _check_power, - EXPR.NegationExpression: _check_other_expression, - EXPR.NPV_NegationExpression: _check_other_expression, - EXPR.AbsExpression: _check_other_expression, - EXPR.NPV_AbsExpression: _check_other_expression, - EXPR.UnaryFunctionExpression: _check_other_expression, - EXPR.NPV_UnaryFunctionExpression: _check_other_expression, + EXPR.NegationExpression: _check_negation, + EXPR.NPV_NegationExpression: _check_negation, + EXPR.AbsExpression: _check_abs, + EXPR.NPV_AbsExpression: _check_abs, + EXPR.UnaryFunctionExpression: _check_unary_function, + EXPR.NPV_UnaryFunctionExpression: _check_unary_function, EXPR.Expr_ifExpression: _check_other_expression, EXPR.ExternalFunctionExpression: _check_other_expression, EXPR.NPV_ExternalFunctionExpression: _check_other_expression, diff --git a/idaes/core/util/tests/test_model_diagnostics.py b/idaes/core/util/tests/test_model_diagnostics.py index 1fe0c37e97..387258fefb 100644 --- a/idaes/core/util/tests/test_model_diagnostics.py +++ b/idaes/core/util/tests/test_model_diagnostics.py @@ -29,6 +29,7 @@ Constraint, Expression, log, + log10, tan, asin, acos, @@ -4454,107 +4455,155 @@ def test_pow_sum_expr_w_negation(self): assert cc == ["(v1 + v2)**(v3 - v4)"] assert not k - # @pytest.mark.unit - # def test_negation_expr(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=-m.scalar_var - # ) == [-21] - # - # @pytest.mark.unit - # def test_negation_sum_expr(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=-(m.scalar_var + m.indexed_var["a"]) - # ) == [-21, -22] - # - # @pytest.mark.unit - # def test_log_expr(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.log(m.scalar_var) - # ) == [pytest.approx(math.log(21), rel=1e-12)] - # - # @pytest.mark.unit - # def test_log_expr_error(self, m): - # with pytest.raises( - # ValueError, - # match="Evaluation error occurred when getting nominal value in log expression " - # "with input 0.0. You should check you scaling factors and model to " - # "address any numerical issues or scale this constraint manually.", - # ): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.log(m.scalar_var - 21) - # ) - # - # @pytest.mark.unit - # def test_log_sum_expr(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.log(m.scalar_var + m.indexed_var["a"]) - # ) == [pytest.approx(math.log(21 + 22), rel=1e-12)] - # - # @pytest.mark.unit - # def test_log_sum_expr_w_negation(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.log(-m.scalar_var + m.indexed_var["a"]) - # ) == [pytest.approx(math.log(-21 + 22), rel=1e-12)] - # - # @pytest.mark.unit - # def test_log10_expr(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.log10(m.scalar_var) - # ) == [pytest.approx(math.log10(21), rel=1e-12)] - # - # @pytest.mark.unit - # def test_log10_sum_expr(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.log10(m.scalar_var + m.indexed_var["a"]) - # ) == [pytest.approx(math.log10(21 + 22), rel=1e-12)] - # - # @pytest.mark.unit - # def test_log10_sum_expr_w_negation(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.log10(-m.scalar_var + m.indexed_var["a"]) - # ) == [pytest.approx(math.log10(-21 + 22), rel=1e-12)] - # - # @pytest.mark.unit - # def test_log10_expr_error(self, m): - # with pytest.raises( - # ValueError, - # match="Evaluation error occurred when getting nominal value in log10 expression " - # "with input 0.0. You should check you scaling factors and model to " - # "address any numerical issues or scale this constraint manually.", - # ): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.log10(m.scalar_var - 21) - # ) - # - # @pytest.mark.unit - # def test_sqrt_expr(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.sqrt(m.scalar_var) - # ) == [pytest.approx(21 ** 0.5, rel=1e-12)] - # - # @pytest.mark.unit - # def test_sqrt_sum_expr(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.sqrt(m.scalar_var + m.indexed_var["a"]) - # ) == [pytest.approx((21 + 22) ** 0.5, rel=1e-12)] - # - # @pytest.mark.unit - # def test_sqrt_sum_expr_w_negation(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.sqrt(-m.scalar_var + m.indexed_var["a"]) - # ) == [pytest.approx((-21 + 22) ** 0.5, rel=1e-12)] - # - # @pytest.mark.unit - # def test_sqrt_expr_error(self, m): - # with pytest.raises( - # ValueError, - # match="Evaluation error occurred when getting nominal value in sqrt expression " - # "with input -21.0. You should check you scaling factors and model to " - # "address any numerical issues or scale this constraint manually.", - # ): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.sqrt(-m.scalar_var) - # ) + @pytest.mark.unit + def test_negation_expr(self): + m = ConcreteModel() + m.v1 = Var(initialize=2) + + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression(expr=-m.v1) + + assert vv == [-2] + assert mm == [] + assert cc == [] + assert not k + + @pytest.mark.unit + def test_negation_sum_expr(self): + m = ConcreteModel() + m.v1 = Var(initialize=2.00001) + m.v2 = Var(initialize=2) + + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( + expr=-(m.v1 - m.v2) + ) + + assert vv == [pytest.approx(-0.00001, rel=1e-8)] + assert mm == [] + assert cc == ["- (v1 - v2)"] + assert not k + + @pytest.mark.unit + def test_log_expr(self): + m = ConcreteModel() + m.v1 = Var(initialize=2) + + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression(expr=log(m.v1)) + + assert vv == [pytest.approx(log(2), rel=1e-8)] + assert mm == [] + assert cc == [] + assert not k + + @pytest.mark.unit + def test_log_expr_error(self): + m = ConcreteModel() + m.v1 = Var(initialize=-2) + + with pytest.raises( + ValueError, + match=re.escape( + "Error in ConstraintTermAnalysisVisitor: error evaluating log(v1)" + ), + ): + ConstraintTermAnalysisVisitor().walk_expression(expr=log(m.v1)) + + @pytest.mark.unit + def test_log_sum_expr(self): + m = ConcreteModel() + m.v1 = Var(initialize=2.00001) + m.v2 = Var(initialize=2) + + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( + expr=log(m.v1 - m.v2) + ) + + assert vv == [pytest.approx(log(0.00001), rel=1e-8)] + assert mm == [] + assert cc == ["log(v1 - v2)"] + assert not k + + @pytest.mark.unit + def test_log10_expr(self): + m = ConcreteModel() + m.v1 = Var(initialize=2) + + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( + expr=log10(m.v1) + ) + + assert vv == [pytest.approx(log10(2), rel=1e-8)] + assert mm == [] + assert cc == [] + assert not k + + @pytest.mark.unit + def test_log10_expr_error(self): + m = ConcreteModel() + m.v1 = Var(initialize=-2) + + with pytest.raises( + ValueError, + match=re.escape( + "Error in ConstraintTermAnalysisVisitor: error evaluating log10(v1)" + ), + ): + ConstraintTermAnalysisVisitor().walk_expression(expr=log10(m.v1)) + + @pytest.mark.unit + def test_log10_sum_expr(self): + m = ConcreteModel() + m.v1 = Var(initialize=2.00001) + m.v2 = Var(initialize=2) + + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( + expr=log10(m.v1 - m.v2) + ) + + assert vv == [pytest.approx(log10(0.00001), rel=1e-8)] + assert mm == [] + assert cc == ["log10(v1 - v2)"] + assert not k + + @pytest.mark.unit + def test_sqrt_expr(self): + m = ConcreteModel() + m.v1 = Var(initialize=2) + + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression(expr=sqrt(m.v1)) + + assert vv == [pytest.approx(sqrt(2), rel=1e-8)] + assert mm == [] + assert cc == [] + assert not k + + @pytest.mark.unit + def test_sqrt_expr_error(self): + m = ConcreteModel() + m.v1 = Var(initialize=-2) + + with pytest.raises( + ValueError, + match=re.escape( + "Error in ConstraintTermAnalysisVisitor: error evaluating sqrt(v1)" + ), + ): + ConstraintTermAnalysisVisitor().walk_expression(expr=sqrt(m.v1)) + + @pytest.mark.unit + def test_sqrt_sum_expr(self): + m = ConcreteModel() + m.v1 = Var(initialize=2.00001) + m.v2 = Var(initialize=2) + + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( + expr=sqrt(m.v1 - m.v2) + ) + + assert vv == [pytest.approx(sqrt(0.00001), rel=1e-8)] + assert mm == [] + assert cc == ["sqrt(v1 - v2)"] + assert not k + # # @pytest.mark.unit # def test_sin_expr(self, m): @@ -4833,24 +4882,6 @@ def test_pow_sum_expr_w_negation(self): # m.a.set_value(2) # m.b.set_value(4) # assert pyo.value(Z) == pytest.approx(expected_mag, rel=1e-8) - # - # @pytest.mark.unit - # def test_Expression(self, m): - # m.expression = pyo.Expression( - # expr=m.scalar_param ** (sum(m.indexed_var[i] for i in m.set)) - # ) - # - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=m.expression - # ) == [0.5 ** (22 + 23 + 24)] - # - # @pytest.mark.unit - # def test_constraint(self, m): - # m.constraint = pyo.Constraint(expr=m.scalar_var == m.expression) - # - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=m.constraint.expr - # ) == [21, 0.5 ** (22 + 23 + 24)] @pytest.mark.unit def test_equality_sum_expr(self): @@ -5015,3 +5046,5 @@ def test_cancelling_equality_expr_compound(self): assert mm == [] assert cc == [] assert not k + + # TODO: Check for a+eps=c forms From 707de1e89f3d3ada0c7d078d6cfb7afa88c1724f Mon Sep 17 00:00:00 2001 From: Andrew Lee Date: Fri, 9 Aug 2024 09:58:50 -0400 Subject: [PATCH 05/20] Reducing duplicated code in tests --- .../core/util/tests/test_model_diagnostics.py | 204 +++++++----------- 1 file changed, 72 insertions(+), 132 deletions(-) diff --git a/idaes/core/util/tests/test_model_diagnostics.py b/idaes/core/util/tests/test_model_diagnostics.py index 387258fefb..3e44015e46 100644 --- a/idaes/core/util/tests/test_model_diagnostics.py +++ b/idaes/core/util/tests/test_model_diagnostics.py @@ -4293,99 +4293,76 @@ def test_sum_expr_constant(self): # expr=m.scalar_var + m.indexed_var["a"] - m.scalar_param # ) == [21, 22, -12] - @pytest.mark.unit - def test_product_expr(self): + @pytest.fixture(scope="class") + def model(self): m = ConcreteModel() m.v1 = Var(initialize=2) m.v2 = Var(initialize=3) - m.v3 = Var(initialize=5) + m.v3 = Var(initialize=5.0000001) + m.v4 = Var(initialize=5) + + return m + @pytest.mark.unit + def test_product_expr(self, model): + m = ConcreteModel() vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( - expr=m.v1 * m.v2 * m.v3 + expr=model.v1 * model.v2 * model.v3 ) - assert vv == [30] + assert vv == [pytest.approx(30.0000006, rel=1e-8)] assert mm == [] assert cc == [] assert not k @pytest.mark.unit - def test_product_sum_expr(self): - m = ConcreteModel() - m.v1 = Var(initialize=2) - m.v2 = Var(initialize=3) - m.v3 = Var(initialize=5) - m.v4 = Var(initialize=6) - + def test_product_sum_expr(self, model): vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( - expr=(m.v1 + m.v2) * (m.v3 + m.v4) + expr=(model.v1 + model.v2) * (model.v3 + model.v4) ) - assert vv == [(2 + 3) * (5 + 6)] + assert vv == [pytest.approx((2 + 3) * (5.0000001 + 5), rel=1e-8)] assert mm == [] assert cc == [] assert not k @pytest.mark.unit - def test_product_sum_expr_w_negation(self): - m = ConcreteModel() - m.v1 = Var(initialize=2) - m.v2 = Var(initialize=3) - m.v3 = Var(initialize=5) - m.v4 = Var(initialize=5) - + def test_product_sum_expr_w_negation(self, model): vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( - expr=(m.v1 + m.v2) * (m.v3 - m.v4) + expr=(model.v1 + model.v2) * (model.v3 - model.v4) ) - assert vv == [0] + assert vv == [pytest.approx(0.0000005, rel=1e-8)] assert mm == [] assert cc == ["(v1 + v2)*(v3 - v4)"] assert not k @pytest.mark.unit - def test_division_expr(self): - m = ConcreteModel() - m.v1 = Var(initialize=2) - m.v2 = Var(initialize=3) - m.v3 = Var(initialize=5) - + def test_division_expr(self, model): vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( - expr=m.v1 / m.v2 / m.v3 + expr=model.v1 / model.v2 / model.v3 ) - assert vv == [pytest.approx(2 / 3 / 5, rel=1e-8)] + assert vv == [pytest.approx(2 / 3 / 5.0000001, rel=1e-8)] assert mm == [] assert cc == [] assert not k @pytest.mark.unit - def test_division_sum_expr(self): - m = ConcreteModel() - m.v1 = Var(initialize=2) - m.v2 = Var(initialize=3) - m.v3 = Var(initialize=5) - m.v4 = Var(initialize=6) - + def test_division_sum_expr(self, model): vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( - expr=(m.v1 + m.v2) / (m.v3 + m.v4) + expr=(model.v1 + model.v2) / (model.v3 + model.v4) ) - assert vv == [pytest.approx((2 + 3) / (5 + 6), rel=1e-8)] + assert vv == [pytest.approx((2 + 3) / (5.0000001 + 5), rel=1e-8)] assert mm == [] assert cc == [] assert not k @pytest.mark.unit - def test_division_sum_expr_w_negation(self): - m = ConcreteModel() - m.v1 = Var(initialize=2) - m.v2 = Var(initialize=3) - m.v3 = Var(initialize=5.0000001) - m.v4 = Var(initialize=5) - + def test_division_sum_expr_w_negation(self, model): vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( - expr=(m.v1 + m.v2) / (m.v3 - m.v4) + expr=(model.v1 + model.v2) / (model.v3 - model.v4) ) assert vv == [pytest.approx((2 + 3) / (0.0000001), rel=1e-8)] @@ -4409,12 +4386,10 @@ def test_division_expr_error(self): ConstraintTermAnalysisVisitor().walk_expression(expr=m.v1 / m.v2) @pytest.mark.unit - def test_pow_expr(self): - m = ConcreteModel() - m.v1 = Var(initialize=2) - m.v2 = Var(initialize=3) - - vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression(expr=m.v1**m.v2) + def test_pow_expr(self, model): + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( + expr=model.v1**model.v2 + ) assert vv == [8] assert mm == [] @@ -4422,32 +4397,20 @@ def test_pow_expr(self): assert not k @pytest.mark.unit - def test_pow_sum_expr(self): - m = ConcreteModel() - m.v1 = Var(initialize=2) - m.v2 = Var(initialize=3) - m.v3 = Var(initialize=5) - m.v4 = Var(initialize=6) - + def test_pow_sum_expr(self, model): vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( - expr=(m.v1 + m.v2) ** (m.v3 + m.v4) + expr=(model.v1 + model.v2) ** (model.v3 + model.v4) ) - assert vv == [5**11] + assert vv == [pytest.approx(5**10.0000001, rel=1e-8)] assert mm == [] assert cc == [] assert not k @pytest.mark.unit - def test_pow_sum_expr_w_negation(self): - m = ConcreteModel() - m.v1 = Var(initialize=2) - m.v2 = Var(initialize=3) - m.v3 = Var(initialize=5.0000001) - m.v4 = Var(initialize=5) - + def test_pow_sum_expr_w_negation(self, model): vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( - expr=(m.v1 + m.v2) ** (m.v3 - m.v4) + expr=(model.v1 + model.v2) ** (model.v3 - model.v4) ) assert vv == [pytest.approx((2 + 3) ** (0.0000001), rel=1e-8)] @@ -4455,12 +4418,19 @@ def test_pow_sum_expr_w_negation(self): assert cc == ["(v1 + v2)**(v3 - v4)"] assert not k - @pytest.mark.unit - def test_negation_expr(self): + @pytest.fixture(scope="class") + def func_model(self): m = ConcreteModel() - m.v1 = Var(initialize=2) + m.v1 = Var(initialize=2.00001) + m.v2 = Var(initialize=2) + + return m - vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression(expr=-m.v1) + @pytest.mark.unit + def test_negation_expr(self, func_model): + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( + expr=-func_model.v2 + ) assert vv == [-2] assert mm == [] @@ -4468,13 +4438,9 @@ def test_negation_expr(self): assert not k @pytest.mark.unit - def test_negation_sum_expr(self): - m = ConcreteModel() - m.v1 = Var(initialize=2.00001) - m.v2 = Var(initialize=2) - + def test_negation_sum_expr(self, func_model): vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( - expr=-(m.v1 - m.v2) + expr=-(func_model.v1 - func_model.v2) ) assert vv == [pytest.approx(-0.00001, rel=1e-8)] @@ -4483,11 +4449,10 @@ def test_negation_sum_expr(self): assert not k @pytest.mark.unit - def test_log_expr(self): - m = ConcreteModel() - m.v1 = Var(initialize=2) - - vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression(expr=log(m.v1)) + def test_log_expr(self, func_model): + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( + expr=log(func_model.v2) + ) assert vv == [pytest.approx(log(2), rel=1e-8)] assert mm == [] @@ -4495,26 +4460,19 @@ def test_log_expr(self): assert not k @pytest.mark.unit - def test_log_expr_error(self): - m = ConcreteModel() - m.v1 = Var(initialize=-2) - + def test_log_expr_error(self, func_model): with pytest.raises( ValueError, match=re.escape( - "Error in ConstraintTermAnalysisVisitor: error evaluating log(v1)" + "Error in ConstraintTermAnalysisVisitor: error evaluating log(- v1)" ), ): - ConstraintTermAnalysisVisitor().walk_expression(expr=log(m.v1)) + ConstraintTermAnalysisVisitor().walk_expression(expr=log(-func_model.v1)) @pytest.mark.unit - def test_log_sum_expr(self): - m = ConcreteModel() - m.v1 = Var(initialize=2.00001) - m.v2 = Var(initialize=2) - + def test_log_sum_expr(self, func_model): vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( - expr=log(m.v1 - m.v2) + expr=log(func_model.v1 - func_model.v2) ) assert vv == [pytest.approx(log(0.00001), rel=1e-8)] @@ -4523,12 +4481,9 @@ def test_log_sum_expr(self): assert not k @pytest.mark.unit - def test_log10_expr(self): - m = ConcreteModel() - m.v1 = Var(initialize=2) - + def test_log10_expr(self, func_model): vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( - expr=log10(m.v1) + expr=log10(func_model.v2) ) assert vv == [pytest.approx(log10(2), rel=1e-8)] @@ -4537,26 +4492,19 @@ def test_log10_expr(self): assert not k @pytest.mark.unit - def test_log10_expr_error(self): - m = ConcreteModel() - m.v1 = Var(initialize=-2) - + def test_log10_expr_error(self, func_model): with pytest.raises( ValueError, match=re.escape( - "Error in ConstraintTermAnalysisVisitor: error evaluating log10(v1)" + "Error in ConstraintTermAnalysisVisitor: error evaluating log10(- v1)" ), ): - ConstraintTermAnalysisVisitor().walk_expression(expr=log10(m.v1)) + ConstraintTermAnalysisVisitor().walk_expression(expr=log10(-func_model.v1)) @pytest.mark.unit - def test_log10_sum_expr(self): - m = ConcreteModel() - m.v1 = Var(initialize=2.00001) - m.v2 = Var(initialize=2) - + def test_log10_sum_expr(self, func_model): vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( - expr=log10(m.v1 - m.v2) + expr=log10(func_model.v1 - func_model.v2) ) assert vv == [pytest.approx(log10(0.00001), rel=1e-8)] @@ -4565,11 +4513,10 @@ def test_log10_sum_expr(self): assert not k @pytest.mark.unit - def test_sqrt_expr(self): - m = ConcreteModel() - m.v1 = Var(initialize=2) - - vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression(expr=sqrt(m.v1)) + def test_sqrt_expr(self, func_model): + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( + expr=sqrt(func_model.v2) + ) assert vv == [pytest.approx(sqrt(2), rel=1e-8)] assert mm == [] @@ -4577,26 +4524,19 @@ def test_sqrt_expr(self): assert not k @pytest.mark.unit - def test_sqrt_expr_error(self): - m = ConcreteModel() - m.v1 = Var(initialize=-2) - + def test_sqrt_expr_error(self, func_model): with pytest.raises( ValueError, match=re.escape( - "Error in ConstraintTermAnalysisVisitor: error evaluating sqrt(v1)" + "Error in ConstraintTermAnalysisVisitor: error evaluating sqrt(- v1)" ), ): - ConstraintTermAnalysisVisitor().walk_expression(expr=sqrt(m.v1)) + ConstraintTermAnalysisVisitor().walk_expression(expr=sqrt(-func_model.v1)) @pytest.mark.unit - def test_sqrt_sum_expr(self): - m = ConcreteModel() - m.v1 = Var(initialize=2.00001) - m.v2 = Var(initialize=2) - + def test_sqrt_sum_expr(self, func_model): vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( - expr=sqrt(m.v1 - m.v2) + expr=sqrt(func_model.v1 - func_model.v2) ) assert vv == [pytest.approx(sqrt(0.00001), rel=1e-8)] From 3a6400b8e879610f83c2bfdc29e77b169e075f7a Mon Sep 17 00:00:00 2001 From: Andrew Lee Date: Fri, 9 Aug 2024 10:56:45 -0400 Subject: [PATCH 06/20] Finishing testing --- idaes/core/util/model_diagnostics.py | 36 +- .../core/util/tests/test_model_diagnostics.py | 445 ++++-------------- 2 files changed, 112 insertions(+), 369 deletions(-) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index ab3e1252ec..cefff5484f 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -4258,29 +4258,17 @@ def _check_unary_function(self, node, child_data): return [val], mismatch, cancelling, const - # def _get_nominal_value_expr_if(self, node, child_nominal_values): - # assert len(child_nominal_values) == 3 - # return child_nominal_values[1] + child_nominal_values[2] - # - # def _get_nominal_value_external_function(self, node, child_nominal_values): - # # First, need to get expected magnitudes of input terms, which may be sub-expressions - # input_mag = [ - # self._get_nominal_value_for_sum_subexpression(i) - # for i in child_nominal_values - # ] - # - # # Next, create a copy of the external function with expected magnitudes as inputs - # newfunc = node.create_node_with_local_data(input_mag) - # - # # Evaluate new function and return the absolute value - # return [pyo.value(newfunc)] - - def _check_other_expression(self, node, child_data): + def _check_other_function(self, node, child_data): mismatch, cancelling, const = self._perform_checks(node, child_data) - # TODO: See if we can be more efficient about getting the value - might need node specific methods - # [value], [mismatched terms], [canceling terms], constant - return [value(node)], mismatch, cancelling, const + # 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 external function with expected magnitudes as inputs + newfunc = node.create_node_with_local_data(input_mag) + + # Evaluate new function and return the absolute value + return [value(newfunc)], mismatch, cancelling, const def _perform_checks(self, node, child_data): # Perform checks for problematic expressions @@ -4358,9 +4346,9 @@ def _sum_combinations(self, values_list): EXPR.NPV_AbsExpression: _check_abs, EXPR.UnaryFunctionExpression: _check_unary_function, EXPR.NPV_UnaryFunctionExpression: _check_unary_function, - EXPR.Expr_ifExpression: _check_other_expression, - EXPR.ExternalFunctionExpression: _check_other_expression, - EXPR.NPV_ExternalFunctionExpression: _check_other_expression, + EXPR.Expr_ifExpression: _check_other_function, + EXPR.ExternalFunctionExpression: _check_other_function, + EXPR.NPV_ExternalFunctionExpression: _check_other_function, EXPR.LinearExpression: _check_sum_expression, } diff --git a/idaes/core/util/tests/test_model_diagnostics.py b/idaes/core/util/tests/test_model_diagnostics.py index 3e44015e46..9b54360451 100644 --- a/idaes/core/util/tests/test_model_diagnostics.py +++ b/idaes/core/util/tests/test_model_diagnostics.py @@ -21,6 +21,7 @@ import os from copy import deepcopy +from unittest import TestCase from pandas import DataFrame from pyomo.environ import ( @@ -28,24 +29,33 @@ ConcreteModel, Constraint, Expression, - log, - log10, - 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 @@ -58,8 +68,11 @@ 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, @@ -4421,18 +4434,21 @@ def test_pow_sum_expr_w_negation(self, model): @pytest.fixture(scope="class") def func_model(self): m = ConcreteModel() - m.v1 = Var(initialize=2.00001) - m.v2 = Var(initialize=2) + 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.v2 + expr=-func_model.v1 ) - assert vv == [-2] + assert vv == [-1] assert mm == [] assert cc == [] assert not k @@ -4448,380 +4464,119 @@ def test_negation_sum_expr(self, func_model): assert cc == ["- (v1 - v2)"] 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 - def test_log_expr(self, func_model): + @pytest.mark.parametrize("func", func_list) + def test_func_expr(self, func_model, func): vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( - expr=log(func_model.v2) + expr=func(func_model.v2) ) - assert vv == [pytest.approx(log(2), rel=1e-8)] + assert vv == [pytest.approx(value(func(0.99999)), rel=1e-8)] assert mm == [] assert cc == [] - assert not k - - @pytest.mark.unit - def test_log_expr_error(self, func_model): - with pytest.raises( - ValueError, - match=re.escape( - "Error in ConstraintTermAnalysisVisitor: error evaluating log(- v1)" - ), - ): - ConstraintTermAnalysisVisitor().walk_expression(expr=log(-func_model.v1)) + assert k @pytest.mark.unit - def test_log_sum_expr(self, func_model): + @pytest.mark.parametrize("func", func_list) + def test_func_sum_expr(self, func_model, func): vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( - expr=log(func_model.v1 - func_model.v2) + expr=func(func_model.v1 - func_model.v2) ) - assert vv == [pytest.approx(log(0.00001), rel=1e-8)] + assert vv == [pytest.approx(value(func(0.00001)), rel=1e-8)] assert mm == [] - assert cc == ["log(v1 - v2)"] + assert cc == [str(func(func_model.v1 - func_model.v2))] assert not k @pytest.mark.unit - def test_log10_expr(self, func_model): + @pytest.mark.parametrize("func", func_list) + def test_func_sum_expr_w_negation(self, func_model, func): vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( - expr=log10(func_model.v2) + expr=func(func_model.v1 - func_model.v2) ) - assert vv == [pytest.approx(log10(2), rel=1e-8)] + assert vv == [pytest.approx(value(func(0.00001)), rel=1e-8)] assert mm == [] - assert cc == [] + assert cc == [str(func(func_model.v1 - func_model.v2))] assert not k @pytest.mark.unit - def test_log10_expr_error(self, func_model): + @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 log10(- v1)" + "Error in ConstraintTermAnalysisVisitor: error evaluating " ), ): - ConstraintTermAnalysisVisitor().walk_expression(expr=log10(-func_model.v1)) + ConstraintTermAnalysisVisitor().walk_expression(expr=func(func_model.v3)) @pytest.mark.unit - def test_log10_sum_expr(self, func_model): - vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( - expr=log10(func_model.v1 - func_model.v2) + def test_expr_if(self, model): + model.exprif = Expr_if( + IF=model.v1, + THEN=model.v1 + model.v2, + ELSE=model.v3 - model.v4, ) - assert vv == [pytest.approx(log10(0.00001), rel=1e-8)] - assert mm == [] - assert cc == ["log10(v1 - v2)"] - assert not k - - @pytest.mark.unit - def test_sqrt_expr(self, func_model): vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( - expr=sqrt(func_model.v2) + expr=model.exprif ) - assert vv == [pytest.approx(sqrt(2), rel=1e-8)] + assert vv == [pytest.approx(5, rel=1e-8)] assert mm == [] - assert cc == [] + assert cc == ["Expr_if( ( v1 ), then=( v1 + v2 ), else=( v3 - v4 ) )"] assert not k @pytest.mark.unit - def test_sqrt_expr_error(self, func_model): - with pytest.raises( - ValueError, - match=re.escape( - "Error in ConstraintTermAnalysisVisitor: error evaluating sqrt(- v1)" - ), - ): - ConstraintTermAnalysisVisitor().walk_expression(expr=sqrt(-func_model.v1)) + @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) - @pytest.mark.unit - def test_sqrt_sum_expr(self, func_model): - vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( - expr=sqrt(func_model.v1 - func_model.v2) - ) + m.expr_write = CubicThermoExpressions(m) + Z = m.expr_write.z_liq(eos=CubicEoS.PR, A=m.a, B=m.b) - assert vv == [pytest.approx(sqrt(0.00001), rel=1e-8)] + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression(expr=Z) + + assert vv == [pytest.approx(-2.1149075414767577, rel=1e-8)] assert mm == [] - assert cc == ["sqrt(v1 - v2)"] + assert cc == [ + "- (1.0 + b - 2*b)", + "compress_fact_liq_func(- (1.0 + b - 2*b), a - b**2 - 2*b - 2*b**2, - a*b + b**2 + b**3)", + "compress_fact_liq_func(- (1.0 + b - 2*b), a - b**2 - 2*b - 2*b**2, - a*b + b**2 + b**3)", + ] assert not k - # - # @pytest.mark.unit - # def test_sin_expr(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.sin(m.scalar_var) - # ) == [pytest.approx(math.sin(21), rel=1e-12)] - # - # @pytest.mark.unit - # def test_sin_sum_expr(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.sin(m.scalar_var + m.indexed_var["a"]) - # ) == [pytest.approx(math.sin(21 + 22), rel=1e-12)] - # - # @pytest.mark.unit - # def test_sin_sum_expr_w_negation(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.sin(m.scalar_var - m.indexed_var["a"]) - # ) == [pytest.approx(math.sin(21 - 22), rel=1e-12)] - # - # @pytest.mark.unit - # def test_cos_expr(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.cos(m.scalar_var) - # ) == [pytest.approx(math.cos(21), rel=1e-12)] - # - # @pytest.mark.unit - # def test_cos_sum_expr(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.cos(m.scalar_var + m.indexed_var["a"]) - # ) == [pytest.approx(math.cos(21 + 22), rel=1e-12)] - # - # @pytest.mark.unit - # def test_cos_sum_expr_w_negation(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.cos(m.scalar_var - m.indexed_var["a"]) - # ) == [pytest.approx(math.cos(21 - 22), rel=1e-12)] - # - # @pytest.mark.unit - # def test_tan_expr(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.tan(m.scalar_var) - # ) == [pytest.approx(math.tan(21), rel=1e-12)] - # - # @pytest.mark.unit - # def test_tan_sum_expr(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.tan(m.scalar_var + m.indexed_var["a"]) - # ) == [pytest.approx(math.tan(21 + 22), rel=1e-12)] - # - # @pytest.mark.unit - # def test_tan_sum_expr_w_negation(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.tan(m.scalar_var - m.indexed_var["a"]) - # ) == [pytest.approx(math.tan(21 - 22), rel=1e-12)] - # - # @pytest.mark.unit - # def test_sinh_expr(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.sinh(m.scalar_var) - # ) == [pytest.approx(math.sinh(21), rel=1e-12)] - # - # @pytest.mark.unit - # def test_sinh_sum_expr(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.sinh(m.scalar_var + m.indexed_var["a"]) - # ) == [pytest.approx(math.sinh(21 + 22), rel=1e-12)] - # - # @pytest.mark.unit - # def test_sinh_sum_expr_w_negation(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.sinh(m.scalar_var - m.indexed_var["a"]) - # ) == [pytest.approx(math.sinh(21 - 22), rel=1e-12)] - # - # @pytest.mark.unit - # def test_cosh_expr(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.cosh(m.scalar_var) - # ) == [pytest.approx(math.cosh(21), rel=1e-12)] - # - # @pytest.mark.unit - # def test_cosh_sum_expr(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.cosh(m.scalar_var + m.indexed_var["a"]) - # ) == [pytest.approx(math.cosh(21 + 22), rel=1e-12)] - # - # @pytest.mark.unit - # def test_cosh_sum_expr_w_negation(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.cosh(m.scalar_var - m.indexed_var["a"]) - # ) == [pytest.approx(math.cosh(21 - 22), rel=1e-12)] - # - # @pytest.mark.unit - # def test_tanh_expr(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.tanh(m.scalar_var) - # ) == [pytest.approx(math.tanh(21), rel=1e-12)] - # - # @pytest.mark.unit - # def test_tanh_sum_expr(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.tanh(m.scalar_var + m.indexed_var["a"]) - # ) == [pytest.approx(math.tanh(21 + 22), rel=1e-12)] - # - # @pytest.mark.unit - # def test_tanh_sum_expr_w_negation(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.tanh(m.scalar_var - m.indexed_var["a"]) - # ) == [pytest.approx(math.tanh(21 - 22), rel=1e-12)] - # - # @pytest.mark.unit - # def test_asin_expr(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression(expr=pyo.asin(1)) == [ - # pytest.approx(math.asin(1), rel=1e-12) - # ] - # - # @pytest.mark.unit - # def test_asin_sum_expr(self, m): - # sc.set_scaling_factor(m.scalar_param, 2) - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.asin(0.5 + m.scalar_param) - # ) == [pytest.approx(math.asin(1), rel=1e-12)] - # - # @pytest.mark.unit - # def test_asin_sum_expr_negation(self, m): - # sc.set_scaling_factor(m.scalar_param, 2) - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.asin(0.5 - m.scalar_param) - # ) == [pytest.approx(math.asin(0), rel=1e-12)] - # - # @pytest.mark.unit - # def test_acos_expr(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.acos(m.scalar_param) - # ) == [pytest.approx(math.acos(0.5), rel=1e-12)] - # - # @pytest.mark.unit - # def test_acos_sum_expr(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.acos(0.5 + m.scalar_param) - # ) == [pytest.approx(math.acos(1), rel=1e-12)] - # - # @pytest.mark.unit - # def test_acos_sum_expr_w_negation(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.acos(0.5 - m.scalar_param) - # ) == [pytest.approx(math.acos(0), rel=1e-12)] - # - # @pytest.mark.unit - # def test_asinh_expr(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.asinh(m.scalar_var) - # ) == [pytest.approx(math.asinh(21), rel=1e-12)] - # - # @pytest.mark.unit - # def test_asinh_sum_expr(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.asinh(m.scalar_var + m.indexed_var["a"]) - # ) == [pytest.approx(math.asinh(21 + 22), rel=1e-12)] - # - # @pytest.mark.unit - # def test_asinh_sum_expr_w_negation(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.asinh(m.scalar_var - m.indexed_var["a"]) - # ) == [pytest.approx(math.asinh(21 - 22), rel=1e-12)] - # - # @pytest.mark.unit - # def test_acosh_expr(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.acosh(m.scalar_var) - # ) == [pytest.approx(math.acosh(21), rel=1e-12)] - # - # @pytest.mark.unit - # def test_acosh_sum_expr(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.acosh(m.scalar_var + m.indexed_var["a"]) - # ) == [pytest.approx(math.acosh(21 + 22), rel=1e-12)] - # - # @pytest.mark.unit - # def test_acosh_sum_expr_w_negation(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.acosh(-m.scalar_var + m.indexed_var["a"]) - # ) == [pytest.approx(math.acosh(-21 + 22), rel=1e-12)] - # - # @pytest.mark.unit - # def test_atanh_expr(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.atanh(m.scalar_param) - # ) == [pytest.approx(math.atanh(0.5), rel=1e-12)] - # - # @pytest.mark.unit - # def test_atanh_sum_expr(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.atanh(0.4 + m.scalar_param) - # ) == [pytest.approx(math.atanh(0.9), rel=1e-12)] - # - # @pytest.mark.unit - # def test_atanh_sum_expr_w_negation(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.atanh(0.4 - m.scalar_param) - # ) == [pytest.approx(math.atanh(-0.1), rel=1e-12)] - # - # @pytest.mark.unit - # def test_exp_expr(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.exp(m.scalar_param) - # ) == [pytest.approx(math.exp(0.5), rel=1e-12)] - # - # @pytest.mark.unit - # def test_exp_sum_expr(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.exp(0.4 + m.scalar_param) - # ) == [pytest.approx(math.exp(0.9), rel=1e-12)] - # - # @pytest.mark.unit - # def test_exp_sum_expr_w_negation(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=pyo.exp(-0.4 + m.scalar_param) - # ) == [pytest.approx(math.exp(0.1), rel=1e-12)] - # - # @pytest.mark.unit - # def test_expr_if(self, m): - # m.exprif = pyo.Expr_if( - # IF=m.scalar_param, - # THEN=m.indexed_var["a"], - # ELSE=m.indexed_var["b"] + m.indexed_var["c"], - # ) - # - # assert sc.NominalValueExtractionVisitor().walk_expression(expr=m.exprif) == [ - # 22, - # 23, - # 24, - # ] - # - # @pytest.mark.unit - # def test_expr_if_w_negation(self, m): - # m.exprif = pyo.Expr_if( - # IF=m.scalar_param, - # THEN=m.indexed_var["a"], - # ELSE=m.indexed_var["b"] - m.indexed_var["c"], - # ) - # - # assert sc.NominalValueExtractionVisitor().walk_expression(expr=m.exprif) == [ - # 22, - # 23, - # -24, - # ] - # - # @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 = pyo.ConcreteModel() - # m.a = pyo.Var(initialize=1) - # m.b = pyo.Var(initialize=1) - # - # sc.set_scaling_factor(m.a, 1 / 2) - # sc.set_scaling_factor(m.b, 1 / 4) - # - # m.expr_write = CubicThermoExpressions(m) - # Z = m.expr_write.z_liq(eos=CubicEoS.PR, A=m.a, B=m.b) - # - # expected_mag = -9.489811292072448 - # assert sc.NominalValueExtractionVisitor().walk_expression(expr=Z) == [ - # pytest.approx(expected_mag, rel=1e-8) - # ] - # - # # Check that model state did not change - # assert pyo.value(m.a) == 1 - # assert pyo.value(m.b) == 1 - # assert pyo.value(Z) == pytest.approx(-2.1149075414767577, rel=1e-8) - # - # # Now, change the actual state to the expected magnitudes and confirm result - # m.a.set_value(2) - # m.b.set_value(4) - # assert pyo.value(Z) == pytest.approx(expected_mag, rel=1e-8) + # 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): From c3e315c3b2248cc7f11674d6f47e4c6061359cc9 Mon Sep 17 00:00:00 2001 From: Andrew Lee Date: Fri, 9 Aug 2024 11:02:38 -0400 Subject: [PATCH 07/20] Adding additioanl tests and clean up --- .../core/util/tests/test_model_diagnostics.py | 55 ++++++++++--------- 1 file changed, 30 insertions(+), 25 deletions(-) diff --git a/idaes/core/util/tests/test_model_diagnostics.py b/idaes/core/util/tests/test_model_diagnostics.py index 9b54360451..bb89a86dc8 100644 --- a/idaes/core/util/tests/test_model_diagnostics.py +++ b/idaes/core/util/tests/test_model_diagnostics.py @@ -4282,30 +4282,6 @@ def test_sum_expr_constant(self): assert cc == [] assert k - # @pytest.mark.unit - # def test_inequality_expr(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=m.scalar_var <= m.indexed_var["a"] - # ) == [21, 22] - # - # @pytest.mark.unit - # def test_sum_expr(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=sum(m.indexed_var[i] for i in m.set) - # ) == [22, 23, 24] - # - # @pytest.mark.unit - # def test_additive_expr(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=m.scalar_var + m.indexed_var["a"] + m.scalar_param - # ) == [21, 22, 12] - # - # @pytest.mark.unit - # def test_additive_expr_w_negation(self, m): - # assert sc.NominalValueExtractionVisitor().walk_expression( - # expr=m.scalar_var + m.indexed_var["a"] - m.scalar_param - # ) == [21, 22, -12] - @pytest.fixture(scope="class") def model(self): m = ConcreteModel() @@ -4592,6 +4568,20 @@ def test_equality_sum_expr(self): assert cc == [] 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) + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( + expr=m.v2 <= sum(m.v1[i] for i in m.v1) + ) + + assert vv == [1e-7, 3e7] + assert mm == ["v2 <= v1[a] + v1[b] + v1[c]"] + assert cc == [] + assert not k + @pytest.mark.unit def test_compound_equality_expr_1(self): m = ConcreteModel() @@ -4742,4 +4732,19 @@ def test_cancelling_equality_expr_compound(self): assert cc == [] assert not k - # TODO: Check for a+eps=c forms + # Double check for a+eps=c form gets flagged in some way + @pytest.mark.unit + def test_cancelling_equality_expr_canceling_sides(self): + m = ConcreteModel() + m.v1 = Var(initialize=1) + m.v2 = Var(initialize=1) + m.v3 = Var(initialize=1e-8) + + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( + expr=m.v2 == m.v1 + m.v3 + ) + + assert vv == [1, pytest.approx(1 + 1e-8, abs=1e-8)] + assert mm == ["v1 + v3"] + assert cc == [] + assert not k From 2ab55ca209d23e9a222246793a25152d0dde07a0 Mon Sep 17 00:00:00 2001 From: Andrew Lee Date: Fri, 9 Aug 2024 11:09:05 -0400 Subject: [PATCH 08/20] Running pylint --- idaes/core/util/model_diagnostics.py | 18 +++++++++--------- .../core/util/tests/test_model_diagnostics.py | 10 +++++----- 2 files changed, 14 insertions(+), 14 deletions(-) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index cefff5484f..bc9a37502d 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -4258,7 +4258,7 @@ def _check_unary_function(self, node, child_data): return [val], mismatch, cancelling, const - def _check_other_function(self, node, child_data): + def _check_other_expression(self, node, child_data): mismatch, cancelling, const = self._perform_checks(node, child_data) # First, need to get value of input terms, which may be sub-expressions @@ -4346,9 +4346,9 @@ def _sum_combinations(self, values_list): EXPR.NPV_AbsExpression: _check_abs, EXPR.UnaryFunctionExpression: _check_unary_function, EXPR.NPV_UnaryFunctionExpression: _check_unary_function, - EXPR.Expr_ifExpression: _check_other_function, - EXPR.ExternalFunctionExpression: _check_other_function, - EXPR.NPV_ExternalFunctionExpression: _check_other_function, + EXPR.Expr_ifExpression: _check_other_expression, + EXPR.ExternalFunctionExpression: _check_other_expression, + EXPR.NPV_ExternalFunctionExpression: _check_other_expression, EXPR.LinearExpression: _check_sum_expression, } @@ -4367,14 +4367,14 @@ def exitNode(self, node, data): if node_func is not None: return node_func(self, node, data) - elif not node.is_expression_type(): + if not node.is_expression_type(): # this is a leaf, but not a native type if nodetype is _PyomoUnit: return [1], [], [], True - else: - # Var or Param - return self._check_base_type(node) - # might want to add other common types here + + # 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 ( diff --git a/idaes/core/util/tests/test_model_diagnostics.py b/idaes/core/util/tests/test_model_diagnostics.py index bb89a86dc8..fc182259e3 100644 --- a/idaes/core/util/tests/test_model_diagnostics.py +++ b/idaes/core/util/tests/test_model_diagnostics.py @@ -13,16 +13,16 @@ """ 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, From 2af9478ac9f0f495bbf20b574b0287a1e3b6ef52 Mon Sep 17 00:00:00 2001 From: Andrew Lee Date: Fri, 9 Aug 2024 14:47:59 -0400 Subject: [PATCH 09/20] Adding constraint walker to diagnostics toolbox --- idaes/core/util/model_diagnostics.py | 147 +++++++++++++++++- .../core/util/tests/test_model_diagnostics.py | 117 +++++++++++++- 2 files changed, 258 insertions(+), 6 deletions(-) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index bc9a37502d..ec7b191b23 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -227,6 +227,22 @@ def svd_sparse(jacobian, number_singular_values): description="Absolute tolerance to use when checking constraint residuals.", ), ) +CONFIG.declare( + "constraint_term_mismatch_tolerance", + ConfigValue( + default=1e6, + domain=float, + description="Magnitude difference to use when checking for mismatched additive terms in constraints.", + ), +) +CONFIG.declare( + "constraint_term_cancellation_tolerance", + ConfigValue( + default=1e-4, + domain=float, + description="Absolute tolerance to use when checking for cancelling additive terms in constraints.", + ), +) CONFIG.declare( "variable_large_value_tolerance", ConfigValue( @@ -1057,6 +1073,105 @@ 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 + ): + _, mm, cc, k = walker.walk_expression(c.expr) + + for i in mm: + mismatch.append(f"{c.name}: {i}") + for i in cc: + cancellation.append(f"{c.name}: {i}") + if k: + 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_cancelling_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 + + _, cancelling, _ = self._collect_constraint_mismatches() + + # Write the output + _write_report_section( + stream=stream, + lines_list=cancelling, + title="The following constraints have cancelling 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 ): @@ -1340,6 +1455,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, @@ -4194,8 +4331,14 @@ def _check_sum_expression(self, node, child_data): absvals = [abs(v) for v in vals] vl = max(absvals) vs = min(absvals) - if vl != vs and log10(vl / vs) > self._mm_tol: - mismatch.append(str(node)) + if vl != vs: + if vs == 0: + diff = log10(vl) + else: + diff = log10(vl / vs) + + if diff >= self._mm_tol: + mismatch.append(str(node)) return vals, mismatch, cancelling, const diff --git a/idaes/core/util/tests/test_model_diagnostics.py b/idaes/core/util/tests/test_model_diagnostics.py index fc182259e3..51bd7569f8 100644 --- a/idaes/core/util/tests/test_model_diagnostics.py +++ b/idaes/core/util/tests/test_model_diagnostics.py @@ -1014,6 +1014,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 == ["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_cancelling_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_cancelling_terms(stream=stream) + + expected = """==================================================================================== +The following constraints have cancelling 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 + ==================================================================================== """ @@ -1166,7 +1271,7 @@ 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 @@ -1177,6 +1282,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): @@ -1193,7 +1299,7 @@ def test_collect_numerical_cautions_jacobian(self): cautions = dt._collect_numerical_cautions() - assert len(cautions) == 4 + assert len(cautions) == 5 assert "Caution: 3 Variables with value close to zero (tol=1.0E-08)" in cautions assert ( "Caution: 3 Variables with extreme Jacobian values (<1.0E-04 or >1.0E+04)" @@ -1204,6 +1310,7 @@ def test_collect_numerical_cautions_jacobian(self): in cautions ) assert "Caution: 4 extreme Jacobian Entries (<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_assert_no_structural_warnings(self, model): @@ -1440,12 +1547,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) ------------------------------------------------------------------------------------ @@ -1490,9 +1598,10 @@ def test_report_numerical_issues_jacobian(self): WARNING: 3 pairs of variables are parallel (to tolerance 1.0E-08) ------------------------------------------------------------------------------------ -4 Cautions +5 Cautions Caution: 3 Variables with value close to zero (tol=1.0E-08) + Caution: 1 Constraint with potential cancellation of terms Caution: 3 Variables with extreme Jacobian values (<1.0E-04 or >1.0E+04) Caution: 1 Constraint with extreme Jacobian values (<1.0E-04 or >1.0E+04) Caution: 4 extreme Jacobian Entries (<1.0E-04 or >1.0E+04) From 162f6311c0232d068a224ee502fe65f3fdddebd8 Mon Sep 17 00:00:00 2001 From: Andrew Lee Date: Fri, 23 Aug 2024 10:20:49 -0400 Subject: [PATCH 10/20] Addressing simple review comments --- idaes/core/util/model_diagnostics.py | 133 +++++++----------- .../core/util/tests/test_model_diagnostics.py | 2 +- 2 files changed, 48 insertions(+), 87 deletions(-) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index ec7b191b23..e03574a018 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -69,6 +69,7 @@ ConfigValue, document_kwargs_from_configdict, PositiveInt, + NonNegativeFloat, ) from pyomo.util.check_units import identify_inconsistent_units from pyomo.contrib.incidence_analysis import IncidenceGraphInterface @@ -194,7 +195,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.", ), @@ -203,7 +204,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.", ), @@ -212,7 +213,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 " @@ -223,7 +224,7 @@ 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.", ), ) @@ -231,7 +232,7 @@ def svd_sparse(jacobian, number_singular_values): "constraint_term_mismatch_tolerance", ConfigValue( default=1e6, - domain=float, + domain=NonNegativeFloat, description="Magnitude difference to use when checking for mismatched additive terms in constraints.", ), ) @@ -239,7 +240,7 @@ def svd_sparse(jacobian, number_singular_values): "constraint_term_cancellation_tolerance", ConfigValue( default=1e-4, - domain=float, + domain=NonNegativeFloat, description="Absolute tolerance to use when checking for cancelling additive terms in constraints.", ), ) @@ -247,7 +248,7 @@ def svd_sparse(jacobian, number_singular_values): "variable_large_value_tolerance", ConfigValue( default=1e4, - domain=float, + domain=NonNegativeFloat, description="Absolute tolerance for considering a value to be large.", ), ) @@ -255,7 +256,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.", ), ) @@ -263,7 +264,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.", ), ) @@ -271,7 +272,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.", ), ) @@ -279,7 +280,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.", ), ) @@ -287,7 +288,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.", ), ) @@ -295,7 +296,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.", ), ) @@ -311,7 +312,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", ), ) @@ -319,7 +320,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", ), ) @@ -357,7 +358,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", ), ) @@ -365,7 +366,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", ), @@ -392,7 +393,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.", ), ) @@ -400,7 +401,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.", ), ) @@ -408,7 +409,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.", ), ) @@ -4293,7 +4294,7 @@ def __init__( ): super().__init__() - self._mm_tol = log10(term_mismatch_tolerance) + self._log_mm_tol = log10(term_mismatch_tolerance) self._sum_tol = term_cancellation_tolerance def _check_base_type(self, node): @@ -4337,67 +4338,27 @@ def _check_sum_expression(self, node, child_data): else: diff = log10(vl / vs) - if diff >= self._mm_tol: + if diff >= self._log_mm_tol: mismatch.append(str(node)) return vals, mismatch, cancelling, const - def _check_product(self, node, child_data): + def _check_general_expr(self, node, child_data): mismatch, cancelling, const = self._perform_checks(node, child_data) - val = self._get_value_for_sum_subexpression( - child_data[0] - ) * self._get_value_for_sum_subexpression(child_data[1]) - - return [val], mismatch, cancelling, const - - def _check_division(self, node, child_data): - mismatch, cancelling, const = self._perform_checks(node, child_data) - - numerator = self._get_value_for_sum_subexpression(child_data[0]) - denominator = self._get_value_for_sum_subexpression(child_data[1]) - # TODO: should we look at mismatched magnitudes in num and denom? - if denominator == 0: - raise ValueError( - f"Error in ConstraintTermAnalysisVisitor: found division with denominator of 0 " - f"({str(node)})." - ) - - return [numerator / denominator], mismatch, cancelling, const - - def _check_power(self, node, child_data): - mismatch, cancelling, const = self._perform_checks(node, child_data) - - base = self._get_value_for_sum_subexpression(child_data[0]) - exponent = self._get_value_for_sum_subexpression(child_data[1]) - - return [base**exponent], mismatch, cancelling, const - - def _check_negation(self, node, child_data): - mismatch, cancelling, const = self._perform_checks(node, child_data) - val = -self._get_value_for_sum_subexpression(child_data[0]) - - return [val], mismatch, cancelling, const - - def _check_abs(self, node, child_data): - mismatch, cancelling, const = self._perform_checks(node, child_data) - val = abs(self._get_value_for_sum_subexpression(child_data[0])) - - return [val], mismatch, cancelling, const - - def _check_unary_function(self, node, child_data): - mismatch, cancelling, const = self._perform_checks(node, child_data) - - func_name = node.getname() - func = getattr(math, func_name) - func_val = self._get_value_for_sum_subexpression(child_data[0]) - try: - val = func(func_val) + 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)})." + ) return [val], mismatch, cancelling, const @@ -4407,10 +4368,10 @@ def _check_other_expression(self, 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 external function with expected magnitudes as inputs + # 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 absolute value + # Evaluate new function and return the value along with check results return [value(newfunc)], mismatch, cancelling, const def _perform_checks(self, node, child_data): @@ -4476,19 +4437,19 @@ def _sum_combinations(self, values_list): EXPR.RangedExpression: _check_sum_expression, EXPR.SumExpression: _check_sum_expression, EXPR.NPV_SumExpression: _check_sum_expression, - EXPR.ProductExpression: _check_product, - EXPR.MonomialTermExpression: _check_product, - EXPR.NPV_ProductExpression: _check_product, - EXPR.DivisionExpression: _check_division, - EXPR.NPV_DivisionExpression: _check_division, - EXPR.PowExpression: _check_power, - EXPR.NPV_PowExpression: _check_power, - EXPR.NegationExpression: _check_negation, - EXPR.NPV_NegationExpression: _check_negation, - EXPR.AbsExpression: _check_abs, - EXPR.NPV_AbsExpression: _check_abs, - EXPR.UnaryFunctionExpression: _check_unary_function, - EXPR.NPV_UnaryFunctionExpression: _check_unary_function, + 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, diff --git a/idaes/core/util/tests/test_model_diagnostics.py b/idaes/core/util/tests/test_model_diagnostics.py index 51bd7569f8..3fae854220 100644 --- a/idaes/core/util/tests/test_model_diagnostics.py +++ b/idaes/core/util/tests/test_model_diagnostics.py @@ -4475,7 +4475,7 @@ def test_division_expr_error(self): m.v2 = Var(initialize=0) with pytest.raises( - ValueError, + ZeroDivisionError, match=re.escape( "Error in ConstraintTermAnalysisVisitor: found division with " "denominator of 0 (v1/v2)." From 90d2cd3a9de8f912b809be6fbca04d9deda6a213 Mon Sep 17 00:00:00 2001 From: Andrew Lee Date: Fri, 23 Aug 2024 12:56:31 -0400 Subject: [PATCH 11/20] Accumulate node rather than str(node) --- idaes/core/util/model_diagnostics.py | 89 +++-- .../core/util/tests/test_model_diagnostics.py | 325 +++++++++--------- 2 files changed, 217 insertions(+), 197 deletions(-) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index e03574a018..3f692874a8 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -241,7 +241,7 @@ def svd_sparse(jacobian, number_singular_values): ConfigValue( default=1e-4, domain=NonNegativeFloat, - description="Absolute tolerance to use when checking for cancelling additive terms in constraints.", + description="Absolute tolerance to use when checking for canceling additive terms in constraints.", ), ) CONFIG.declare( @@ -1123,7 +1123,7 @@ def display_constraints_with_mismatched_terms(self, stream=None): footer="=", ) - def display_constraints_with_cancelling_terms(self, stream=None): + def display_constraints_with_canceling_terms(self, stream=None): """ Display constraints in model which contain additive terms which potentially cancel each other. @@ -1137,13 +1137,13 @@ def display_constraints_with_cancelling_terms(self, stream=None): if stream is None: stream = sys.stdout - _, cancelling, _ = self._collect_constraint_mismatches() + _, canceling, _ = self._collect_constraint_mismatches() # Write the output _write_report_section( stream=stream, - lines_list=cancelling, - title="The following constraints have cancelling terms:", + lines_list=canceling, + title="The following constraints have canceling terms:", header="=", footer="=", ) @@ -4294,16 +4294,21 @@ def __init__( ): 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 _check_base_type(self, node): - # [value], [mismatched terms], [canceling terms], constant + # [value], constant if isinstance(node, VarData): const = node.fixed else: const = True - return [value(node)], [], [], const + return [value(node)], const def _get_value_for_sum_subexpression(self, child_data): return sum(i for i in child_data[0]) @@ -4312,19 +4317,15 @@ def _check_sum_expression(self, node, child_data): # Sum expressions need special handling # For sums, collect all child values into a list vals = [] - mismatch = [] # We will check for cancellation in this node at the next level # Pyomo is generally good at simplifying compound sums - cancelling = [] const = True # Collect data from child nodes for d in child_data: vals.append(self._get_value_for_sum_subexpression(d)) - mismatch += d[1] - cancelling += d[2] # Expression is not constant if any child is not constant - if not d[3]: + if not d[1]: const = False # Check for mismatched terms @@ -4339,12 +4340,12 @@ def _check_sum_expression(self, node, child_data): diff = log10(vl / vs) if diff >= self._log_mm_tol: - mismatch.append(str(node)) + self.mismatched_terms.add(node) - return vals, mismatch, cancelling, const + return vals, const def _check_general_expr(self, node, child_data): - mismatch, cancelling, const = self._perform_checks(node, child_data) + const = self._perform_checks(node, child_data) try: val = node._apply_operation( @@ -4360,10 +4361,10 @@ def _check_general_expr(self, node, child_data): f"({str(node)})." ) - return [val], mismatch, cancelling, const + return [val], const def _check_other_expression(self, node, child_data): - mismatch, cancelling, const = self._perform_checks(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] @@ -4372,44 +4373,38 @@ def _check_other_expression(self, node, child_data): newfunc = node.create_node_with_local_data(input_mag) # Evaluate new function and return the value along with check results - return [value(newfunc)], mismatch, cancelling, const + return [value(newfunc)], const 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 - mismatch = [] - cancelling = [] const = True for d in child_data: - # Collect all warnings about mismatched terms from child nodes - mismatch += d[1] - cancelling += d[2] - - # We will check for cancelling terms here, rather than the sum itself, to handle special cases + # 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 sums = self._sum_combinations(d[0]) if any(i <= self._sum_tol * max(d[0]) for i in sums): - cancelling.append(str(node)) + self.canceling_terms.add(node) # Expression is not constant if any child is not constant - if not d[3]: + if not d[1]: const = False # Return any problematic terms found - return mismatch, cancelling, const + return const def _check_equality_expression(self, node, child_data): # (In)equality expressions are a special case of sum expressions # We can start by just calling the method to check the sum expression - vals, mismatch, cancelling, const = self._check_sum_expression(node, child_data) + vals, const = self._check_sum_expression(node, child_data) # 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 - if str(node) not in mismatch and not any(d[3] for d in child_data): + if node not in self.mismatched_terms and not any(d[1] for d in child_data): # No constant terms, check for cancellation # First, collect terms from both sides t = [] @@ -4419,9 +4414,9 @@ def _check_equality_expression(self, node, child_data): # Then check for cancellations sums = self._sum_combinations(t) if any(i <= self._sum_tol * max(t) for i in sums): - cancelling.append(str(node)) + self.canceling_terms.add(node) - return vals, mismatch, cancelling, const + return vals, const def _sum_combinations(self, values_list): sums = [] @@ -4429,6 +4424,7 @@ def _sum_combinations(self, values_list): combinations(values_list, r) for r in range(2, len(values_list) + 1) ): sums.append(abs(sum(i))) + return sums node_type_method_map = { @@ -4460,12 +4456,12 @@ def exitNode(self, node, data): """ Method to call when exiting node to check for potential issues. """ - # Return [node values], [mismatched terms], [cancelling terms], constant + # Return [node values], [constant # first check if the node is a leaf nodetype = type(node) if nodetype in native_types: - return [node], [], [], True + return [node], True node_func = self.node_type_method_map.get(nodetype, None) if node_func is not None: @@ -4474,7 +4470,7 @@ def exitNode(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 + return [1], True # Var or Param return self._check_base_type(node) @@ -4491,3 +4487,26 @@ def exitNode(self, node, data): f"An unhandled expression node type: {str(nodetype)} was encountered while " f"analyzing constraint terms {str(node)}" ) + + def walk_expression(self, expr): + """ + Main method ot call to walk an expression and return analysis. + + Args: + expr - expression ot be analysed + + 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 3fae854220..80ba143bb4 100644 --- a/idaes/core/util/tests/test_model_diagnostics.py +++ b/idaes/core/util/tests/test_model_diagnostics.py @@ -4224,8 +4224,8 @@ def test_int(self): vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression(expr=7) assert vv == [7] - assert mm == [] - assert cc == [] + assert len(mm) == 0 + assert len(cc) == 0 assert k @pytest.mark.unit @@ -4233,8 +4233,8 @@ def test_float(self): vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression(expr=7.7) assert vv == [7.7] - assert mm == [] - assert cc == [] + assert len(mm) == 0 + assert len(cc) == 0 assert k @pytest.mark.unit @@ -4246,8 +4246,8 @@ def test_scalar_param(self): ) assert vv == [1] - assert mm == [] - assert cc == [] + assert len(mm) == 0 + assert len(cc) == 0 assert k @pytest.mark.unit @@ -4259,8 +4259,8 @@ def test_indexed_param(self): expr=m.indexed_param["a"] ) assert vv == [1] - assert mm == [] - assert cc == [] + assert len(mm) == 0 + assert len(cc) == 0 assert k @pytest.mark.unit @@ -4272,8 +4272,8 @@ def test_scalar_var(self): ) assert vv == [1] - assert mm == [] - assert cc == [] + assert len(mm) == 0 + assert len(cc) == 0 assert not k @pytest.mark.unit @@ -4285,8 +4285,8 @@ def test_indexed_var(self): expr=m.indexed_var["a"] ) assert vv == [1] - assert mm == [] - assert cc == [] + assert len(mm) == 0 + assert len(cc) == 0 assert not k @pytest.mark.unit @@ -4299,8 +4299,8 @@ def test_scalar_var_fixed(self): ) assert vv == [1] - assert mm == [] - assert cc == [] + assert len(mm) == 0 + assert len(cc) == 0 assert k @pytest.mark.unit @@ -4313,8 +4313,8 @@ def test_indexed_var_fixed(self): expr=m.indexed_var["a"] ) assert vv == [1] - assert mm == [] - assert cc == [] + assert len(mm) == 0 + assert len(cc) == 0 assert k @pytest.mark.unit @@ -4323,13 +4323,13 @@ def test_equality_expr(self): m.v1 = Var(initialize=1e-7) m.v2 = Var(initialize=1e7) - vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( - expr=m.v1 == m.v2 - ) + expr = m.v1 == m.v2 + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) assert vv == [1e-7, 1e7] - assert mm == ["v1 == v2"] - assert cc == [] + assert expr in mm + assert len(mm) == 1 + assert len(cc) == 0 assert not k @pytest.mark.unit @@ -4341,25 +4341,24 @@ def test_equality_expr_constant(self): # Fix v1, not constant yet as v2 still free m.v1.fix() - vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( - expr=m.v1 == m.v2 - ) + expr = m.v1 == m.v2 + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) assert vv == [1e-7, 1e7] - assert mm == ["v1 == v2"] - assert cc == [] + 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=m.v1 == m.v2 - ) + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) assert vv == [1e-7, 1e7] - assert mm == ["v1 == v2"] - assert cc == [] + assert expr in mm + assert len(mm) == 1 + assert len(cc) == 0 assert k @pytest.mark.unit @@ -4367,13 +4366,14 @@ def test_sum_expr(self): m = ConcreteModel() m.v1 = Var(["a", "b", "c"], initialize=1e7) m.v1["a"].set_value(1e-7) - vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( - expr=sum(m.v1[i] for i in m.v1) - ) + + 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 mm == ["v1[a] + v1[b] + v1[c]"] - assert cc == [] + assert expr in mm + assert len(mm) == 1 + assert len(cc) == 0 assert not k @pytest.mark.unit @@ -4382,13 +4382,14 @@ def test_sum_expr_constant(self): m.v1 = Var(["a", "b", "c"], initialize=1e7) m.v1["a"].set_value(1e-7) m.v1.fix() - vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( - expr=sum(m.v1[i] for i in m.v1) - ) + + 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 mm == ["v1[a] + v1[b] + v1[c]"] - assert cc == [] + assert expr in mm + assert len(mm) == 1 + assert len(cc) == 0 assert k @pytest.fixture(scope="class") @@ -4404,46 +4405,43 @@ def model(self): @pytest.mark.unit def test_product_expr(self, model): m = ConcreteModel() - vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( - expr=model.v1 * model.v2 * model.v3 - ) + 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 mm == [] - assert cc == [] + assert len(mm) == 0 + assert len(cc) == 0 assert not k @pytest.mark.unit def test_product_sum_expr(self, model): - vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( - expr=(model.v1 + model.v2) * (model.v3 + model.v4) - ) + 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 mm == [] - assert cc == [] + assert len(mm) == 0 + assert len(cc) == 0 assert not k @pytest.mark.unit def test_product_sum_expr_w_negation(self, model): - vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( - expr=(model.v1 + model.v2) * (model.v3 - model.v4) - ) + 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 mm == [] - assert cc == ["(v1 + v2)*(v3 - v4)"] + assert len(mm) == 0 + assert expr in cc + assert len(cc) == 1 assert not k @pytest.mark.unit def test_division_expr(self, model): - vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( - expr=model.v1 / model.v2 / model.v3 - ) + 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 mm == [] - assert cc == [] + assert len(mm) == 0 + assert len(cc) == 0 assert not k @pytest.mark.unit @@ -4453,19 +4451,19 @@ def test_division_sum_expr(self, model): ) assert vv == [pytest.approx((2 + 3) / (5.0000001 + 5), rel=1e-8)] - assert mm == [] - assert cc == [] + assert len(mm) == 0 + assert len(cc) == 0 assert not k @pytest.mark.unit def test_division_sum_expr_w_negation(self, model): - vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( - expr=(model.v1 + model.v2) / (model.v3 - model.v4) - ) + 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 mm == [] - assert cc == ["(v1 + v2)/(v3 - v4)"] + assert len(mm) == 0 + assert expr in cc + assert len(cc) == 1 assert not k @pytest.mark.unit @@ -4490,8 +4488,8 @@ def test_pow_expr(self, model): ) assert vv == [8] - assert mm == [] - assert cc == [] + assert len(mm) == 0 + assert len(cc) == 0 assert not k @pytest.mark.unit @@ -4501,19 +4499,19 @@ def test_pow_sum_expr(self, model): ) assert vv == [pytest.approx(5**10.0000001, rel=1e-8)] - assert mm == [] - assert cc == [] + assert len(mm) == 0 + assert len(cc) == 0 assert not k @pytest.mark.unit def test_pow_sum_expr_w_negation(self, model): - vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( - expr=(model.v1 + model.v2) ** (model.v3 - model.v4) - ) + 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 mm == [] - assert cc == ["(v1 + v2)**(v3 - v4)"] + assert len(mm) == 0 + assert expr in cc + assert len(cc) == 1 assert not k @pytest.fixture(scope="class") @@ -4534,19 +4532,19 @@ def test_negation_expr(self, func_model): ) assert vv == [-1] - assert mm == [] - assert cc == [] + assert len(mm) == 0 + assert len(cc) == 0 assert not k @pytest.mark.unit def test_negation_sum_expr(self, func_model): - vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( - expr=-(func_model.v1 - func_model.v2) - ) + 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 mm == [] - assert cc == ["- (v1 - v2)"] + 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 @@ -4577,32 +4575,32 @@ def test_func_expr(self, func_model, func): ) assert vv == [pytest.approx(value(func(0.99999)), rel=1e-8)] - assert mm == [] - assert cc == [] + 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): - vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( - expr=func(func_model.v1 - func_model.v2) - ) + 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 mm == [] - assert cc == [str(func(func_model.v1 - func_model.v2))] + 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): - vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( - expr=func(func_model.v1 - func_model.v2) - ) + 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 mm == [] - assert cc == [str(func(func_model.v1 - func_model.v2))] + assert len(mm) == 0 + assert expr in cc + assert len(cc) == 1 assert not k @pytest.mark.unit @@ -4629,8 +4627,9 @@ def test_expr_if(self, model): ) assert vv == [pytest.approx(5, rel=1e-8)] - assert mm == [] - assert cc == ["Expr_if( ( v1 ), then=( v1 + v2 ), else=( v3 - v4 ) )"] + assert len(mm) == 0 + assert model.exprif in cc + assert len(cc) == 1 assert not k @pytest.mark.unit @@ -4650,12 +4649,10 @@ def test_ext_func(self): vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression(expr=Z) assert vv == [pytest.approx(-2.1149075414767577, rel=1e-8)] - assert mm == [] - assert cc == [ - "- (1.0 + b - 2*b)", - "compress_fact_liq_func(- (1.0 + b - 2*b), a - b**2 - 2*b - 2*b**2, - a*b + b**2 + b**3)", - "compress_fact_liq_func(- (1.0 + b - 2*b), a - b**2 - 2*b - 2*b**2, - a*b + b**2 + b**3)", - ] + 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 @@ -4668,13 +4665,14 @@ def test_equality_sum_expr(self): m = ConcreteModel() m.v1 = Var(["a", "b", "c"], initialize=1e7) m.v2 = Var(initialize=1e-7) - vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( - expr=m.v2 == sum(m.v1[i] for i in m.v1) - ) + + 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 mm == ["v2 == v1[a] + v1[b] + v1[c]"] - assert cc == [] + assert expr in mm + assert len(mm) == 1 + assert len(cc) == 0 assert not k @pytest.mark.unit @@ -4682,13 +4680,14 @@ def test_inequality_sum_expr(self): m = ConcreteModel() m.v1 = Var(["a", "b", "c"], initialize=1e7) m.v2 = Var(initialize=1e-7) - vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( - expr=m.v2 <= sum(m.v1[i] for i in m.v1) - ) + + 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 mm == ["v2 <= v1[a] + v1[b] + v1[c]"] - assert cc == [] + assert expr in mm + assert len(mm) == 1 + assert len(cc) == 0 assert not k @pytest.mark.unit @@ -4696,13 +4695,14 @@ def test_compound_equality_expr_1(self): m = ConcreteModel() m.v1 = Var(["a", "b", "c"], initialize=1e7) m.v2 = Var(initialize=1e-7) - vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( - expr=6 * m.v2 == 8 * sum(m.v1[i] for i in m.v1) - ) + + 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 mm == ["6*v2 == 8*(v1[a] + v1[b] + v1[c])"] - assert cc == [] + assert expr in mm + assert len(mm) == 1 + assert len(cc) == 0 assert not k @pytest.mark.unit @@ -4715,16 +4715,15 @@ def test_compound_equality_expr_2(self): # Set this small so we get two mismatched warnings m.v1["a"].set_value(1e-7) - vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( - expr=6 * m.v2 == 8 * sum(m.v1[i] for i in m.v1) + m.v3 - ) + 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 mm == [ - "v1[a] + v1[b] + v1[c]", - "6*v2 == 8*(v1[a] + v1[b] + v1[c]) + v3", - ] - assert cc == [] + assert expr in mm + assert expr1 in mm + assert len(mm) == 2 + assert len(cc) == 0 assert not k @pytest.mark.unit @@ -4737,9 +4736,9 @@ def test_cancelling_sum_expr(self): ) assert vv == [2, -2] - assert mm == [] + assert len(mm) == 0 # We do not check cancellation at the sum, so this should be empty - assert cc == [] + assert len(cc) == 0 assert not k @pytest.mark.unit @@ -4748,14 +4747,15 @@ def test_expr_w_cancelling_sum(self): m.v1 = Var(initialize=2) m.v2 = Var(initialize=2) m.v3 = Var(initialize=3) - vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( - expr=m.v3 * (m.v1 - m.v2) - ) + + expr = m.v3 * (m.v1 - m.v2) + vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) assert vv == [0] - assert mm == [] + assert len(mm) == 0 # We should get a warning about cancelling sums here - assert cc == ["v3*(v1 - v2)"] + assert expr in cc + assert len(cc) == 1 assert not k # Check for tolerance of sum cancellation @@ -4763,22 +4763,23 @@ def test_expr_w_cancelling_sum(self): vv, mm, cc, k = ConstraintTermAnalysisVisitor( term_cancellation_tolerance=1e-4 - ).walk_expression(expr=m.v3 * (m.v1 - m.v2)) + ).walk_expression(expr=expr) assert vv == [pytest.approx(3 * -0.00022, rel=1e-8)] - assert mm == [] + assert len(mm) == 0 # This should pass as the difference is greater than tol - assert cc == [] + 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=m.v3 * (m.v1 - m.v2)) + ).walk_expression(expr=expr) assert vv == [pytest.approx(3 * -0.00022, rel=1e-8)] - assert mm == [] - assert cc == ["v3*(v1 - v2)"] + assert len(mm) == 0 + assert expr in cc + assert len(cc) == 1 assert not k @pytest.mark.unit @@ -4793,8 +4794,8 @@ def test_cancelling_equality_expr_safe(self): ) assert vv == [0, pytest.approx(0, abs=1e-12)] - assert mm == [] - assert cc == [] + assert len(mm) == 0 + assert len(cc) == 0 assert not k @pytest.mark.unit @@ -4804,13 +4805,13 @@ def test_cancelling_equality_expr_issue(self): m.v2 = Var(initialize=1e7) m.v3 = Var(initialize=0) - vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( - expr=m.v3 == m.v1 - m.v2 - ) + 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 mm == [] - assert cc == ["v3 == v1 - v2"] + assert len(mm) == 0 + assert expr in cc + assert len(cc) == 1 assert not k @pytest.mark.unit @@ -4820,25 +4821,24 @@ def test_cancelling_equality_expr_compound(self): m.v2 = Var(initialize=1e7) m.v3 = Var(initialize=0) - vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( - expr=m.v3 == sum(v for v in m.v1.values()) - m.v2 - ) + 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 mm == [] - assert cc == ["v3 == v1[a] + v1[b] - v2"] + assert len(mm) == 0 + print(cc) + 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=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 mm == [] - assert cc == [] + assert len(mm) == 0 + assert len(cc) == 0 assert not k # Double check for a+eps=c form gets flagged in some way @@ -4849,11 +4849,12 @@ def test_cancelling_equality_expr_canceling_sides(self): m.v2 = Var(initialize=1) m.v3 = Var(initialize=1e-8) - vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression( - expr=m.v2 == m.v1 + m.v3 - ) + 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 mm == ["v1 + v3"] - assert cc == [] + assert expr1 in mm + assert len(mm) == 1 + assert len(cc) == 0 assert not k From 152b044d385b50bf6d6686690d11cf4f602e025f Mon Sep 17 00:00:00 2001 From: Andrew Lee Date: Fri, 23 Aug 2024 14:13:30 -0400 Subject: [PATCH 12/20] Clean up some tests --- .../core/util/tests/test_model_diagnostics.py | 20 +++++++++---------- 1 file changed, 10 insertions(+), 10 deletions(-) diff --git a/idaes/core/util/tests/test_model_diagnostics.py b/idaes/core/util/tests/test_model_diagnostics.py index 80ba143bb4..4496ed3350 100644 --- a/idaes/core/util/tests/test_model_diagnostics.py +++ b/idaes/core/util/tests/test_model_diagnostics.py @@ -1073,7 +1073,7 @@ def test_display_constraints_with_mismatched_terms(self): assert stream.getvalue() == expected @pytest.mark.component - def test_display_constraints_with_cancelling_terms(self): + def test_display_constraints_with_canceling_terms(self): m = ConcreteModel() # Constraint with mismatched terms m.v3 = Var(initialize=10) @@ -1086,10 +1086,10 @@ def test_display_constraints_with_cancelling_terms(self): dt = DiagnosticsToolbox(model=m) stream = StringIO() - dt.display_constraints_with_cancelling_terms(stream=stream) + dt.display_constraints_with_canceling_terms(stream=stream) expected = """==================================================================================== -The following constraints have cancelling terms: +The following constraints have canceling terms: c3: v6 == 10 + v3 - v4 @@ -4727,7 +4727,7 @@ def test_compound_equality_expr_2(self): assert not k @pytest.mark.unit - def test_cancelling_sum_expr(self): + def test_canceling_sum_expr(self): m = ConcreteModel() m.v1 = Var(initialize=2) m.v2 = Var(initialize=2) @@ -4742,7 +4742,7 @@ def test_cancelling_sum_expr(self): assert not k @pytest.mark.unit - def test_expr_w_cancelling_sum(self): + def test_expr_w_canceling_sum(self): m = ConcreteModel() m.v1 = Var(initialize=2) m.v2 = Var(initialize=2) @@ -4753,7 +4753,7 @@ def test_expr_w_cancelling_sum(self): assert vv == [0] assert len(mm) == 0 - # We should get a warning about cancelling sums here + # We should get a warning about canceling sums here assert expr in cc assert len(cc) == 1 assert not k @@ -4783,7 +4783,7 @@ def test_expr_w_cancelling_sum(self): assert not k @pytest.mark.unit - def test_cancelling_equality_expr_safe(self): + def test_canceling_equality_expr_safe(self): m = ConcreteModel() m.v1 = Var(initialize=1e7) m.v2 = Var(initialize=1e7) @@ -4799,7 +4799,7 @@ def test_cancelling_equality_expr_safe(self): assert not k @pytest.mark.unit - def test_cancelling_equality_expr_issue(self): + def test_canceling_equality_expr_issue(self): m = ConcreteModel() m.v1 = Var(initialize=1e7) m.v2 = Var(initialize=1e7) @@ -4815,7 +4815,7 @@ def test_cancelling_equality_expr_issue(self): assert not k @pytest.mark.unit - def test_cancelling_equality_expr_compound(self): + def test_canceling_equality_expr_compound(self): m = ConcreteModel() m.v1 = Var(["a", "b"], initialize=5e6) m.v2 = Var(initialize=1e7) @@ -4843,7 +4843,7 @@ def test_cancelling_equality_expr_compound(self): # Double check for a+eps=c form gets flagged in some way @pytest.mark.unit - def test_cancelling_equality_expr_canceling_sides(self): + def test_canceling_equality_expr_canceling_sides(self): m = ConcreteModel() m.v1 = Var(initialize=1) m.v2 = Var(initialize=1) From 81c18352359a2deb28890fe782e06e1bc3bd2397 Mon Sep 17 00:00:00 2001 From: Andrew Lee Date: Fri, 23 Aug 2024 15:07:14 -0400 Subject: [PATCH 13/20] Supporting ranged expressions --- idaes/core/util/model_diagnostics.py | 35 ++++++++-- .../core/util/tests/test_model_diagnostics.py | 70 +++++++++++++++---- 2 files changed, 86 insertions(+), 19 deletions(-) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index 3f692874a8..c9592fe787 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -4398,17 +4398,28 @@ def _perform_checks(self, node, child_data): def _check_equality_expression(self, node, child_data): # (In)equality expressions are a special case of sum expressions - # We can start by just calling the method to check the sum expression - vals, const = self._check_sum_expression(node, child_data) + # First, need to negate one side of the expression + mdata = [] + for i in [0, 1]: + if i == 0: + vals = [] + for j in child_data[i][0]: + vals.append(-j) + mdata.append((vals, child_data[i][1])) + else: + mdata.append(child_data[i]) + + # 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 - if node not in self.mismatched_terms and not any(d[1] for d in child_data): + if 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 t = [] - for d in child_data: + for d in mdata: t += d[0] # Then check for cancellations @@ -4418,6 +4429,20 @@ def _check_equality_expression(self, node, child_data): return vals, const + def _check_ranged(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 + def _sum_combinations(self, values_list): sums = [] for i in chain.from_iterable( @@ -4430,7 +4455,7 @@ def _sum_combinations(self, values_list): node_type_method_map = { EXPR.EqualityExpression: _check_equality_expression, EXPR.InequalityExpression: _check_equality_expression, - EXPR.RangedExpression: _check_sum_expression, + EXPR.RangedExpression: _check_ranged, EXPR.SumExpression: _check_sum_expression, EXPR.NPV_SumExpression: _check_sum_expression, EXPR.ProductExpression: _check_general_expr, diff --git a/idaes/core/util/tests/test_model_diagnostics.py b/idaes/core/util/tests/test_model_diagnostics.py index 4496ed3350..86bb4fd64e 100644 --- a/idaes/core/util/tests/test_model_diagnostics.py +++ b/idaes/core/util/tests/test_model_diagnostics.py @@ -62,6 +62,7 @@ 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 @@ -1045,7 +1046,7 @@ def test_collect_constraint_mismatches(self): mismatch, cancellation, constant = dt._collect_constraint_mismatches() assert mismatch == ["c2: v4 + v5"] - assert cancellation == ["c3: v6 == 10 + v3 - v4"] + assert cancellation == ["c2: v3 == v4 + v5", "c3: v6 == 10 + v3 - v4"] assert constant == ["c1"] @pytest.mark.component @@ -1271,6 +1272,7 @@ def test_collect_numerical_cautions(self, model): dt = DiagnosticsToolbox(model=model.b) cautions = dt._collect_numerical_cautions() + assert len(cautions) == 6 assert ( "Caution: 2 Variables with value close to their bounds (abs=1.0E-04, rel=1.0E-04)" @@ -1282,7 +1284,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 + assert "Caution: 2 Constraints with potential cancellation of terms" in cautions @pytest.mark.component def test_collect_numerical_cautions_jacobian(self): @@ -1466,9 +1468,9 @@ def test_report_numerical_issues_ok(self): No warnings found! ------------------------------------------------------------------------------------ -0 Cautions +1 Cautions - No cautions found! + Caution: 2 Constraints with potential cancellation of terms ------------------------------------------------------------------------------------ Suggested next steps: @@ -1553,7 +1555,7 @@ def test_report_numerical_issues(self, model): 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: 2 Constraints with potential cancellation of terms Caution: 1 extreme Jacobian Entry (<1.0E-04 or >1.0E+04) ------------------------------------------------------------------------------------ @@ -4326,7 +4328,7 @@ def test_equality_expr(self): expr = m.v1 == m.v2 vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) - assert vv == [1e-7, 1e7] + assert vv == [-1e-7, 1e7] assert expr in mm assert len(mm) == 1 assert len(cc) == 0 @@ -4344,7 +4346,7 @@ def test_equality_expr_constant(self): expr = m.v1 == m.v2 vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) - assert vv == [1e-7, 1e7] + assert vv == [-1e-7, 1e7] assert expr in mm assert len(mm) == 1 assert len(cc) == 0 @@ -4355,7 +4357,7 @@ def test_equality_expr_constant(self): vv, mm, cc, k = ConstraintTermAnalysisVisitor().walk_expression(expr=expr) - assert vv == [1e-7, 1e7] + assert vv == [-1e-7, 1e7] assert expr in mm assert len(mm) == 1 assert len(cc) == 0 @@ -4669,7 +4671,7 @@ def test_equality_sum_expr(self): 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 vv == [-1e-7, 3e7] assert expr in mm assert len(mm) == 1 assert len(cc) == 0 @@ -4684,7 +4686,7 @@ def test_inequality_sum_expr(self): 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 vv == [-1e-7, 3e7] assert expr in mm assert len(mm) == 1 assert len(cc) == 0 @@ -4699,12 +4701,51 @@ def test_compound_equality_expr_1(self): 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 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() @@ -4719,7 +4760,7 @@ def test_compound_equality_expr_2(self): 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 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 @@ -4853,8 +4894,9 @@ def test_canceling_equality_expr_canceling_sides(self): 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 vv == [-1, pytest.approx(1 + 1e-8, abs=1e-8)] assert expr1 in mm assert len(mm) == 1 - assert len(cc) == 0 + assert expr in cc + assert len(cc) == 1 assert not k From bff77d50c7e4e5b0c0d3afffc319122034366394 Mon Sep 17 00:00:00 2001 From: Andrew Lee Date: Fri, 23 Aug 2024 15:27:21 -0400 Subject: [PATCH 14/20] Apply suggestions from code review --- idaes/core/util/model_diagnostics.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index c9592fe787..6b808e7c54 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -4515,10 +4515,10 @@ def exitNode(self, node, data): def walk_expression(self, expr): """ - Main method ot call to walk an expression and return analysis. + Main method to call to walk an expression and return analysis. Args: - expr - expression ot be analysed + expr - expression to be analyzed Returns: list of values of top-level additive terms From 85d1956d3d948bee8c1175c85979a9d8d95eb2b0 Mon Sep 17 00:00:00 2001 From: Andrew Lee Date: Sun, 25 Aug 2024 15:30:48 -0400 Subject: [PATCH 15/20] Regorganising and fixing pylint issues --- idaes/core/util/model_diagnostics.py | 142 +++++++++++++-------------- 1 file changed, 71 insertions(+), 71 deletions(-) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index 6b808e7c54..08099b34ce 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -17,7 +17,6 @@ __author__ = "Alexander Dowling, Douglas Allan, Andrew Lee, Robby Parker, Ben Knueven" -import math from operator import itemgetter import sys from inspect import signature @@ -4313,67 +4312,14 @@ def _check_base_type(self, node): def _get_value_for_sum_subexpression(self, child_data): return sum(i for i in child_data[0]) - 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 - - def _check_general_expr(self, node, child_data): - const = self._perform_checks(node, child_data) - - try: - 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)})." - ) - - return [val], const - - 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) + def _sum_combinations(self, values_list): + sums = [] + for i in chain.from_iterable( + combinations(values_list, r) for r in range(2, len(values_list) + 1) + ): + sums.append(abs(sum(i))) - # Evaluate new function and return the value along with check results - return [value(newfunc)], const + return sums def _perform_checks(self, node, child_data): # Perform checks for problematic expressions @@ -4429,7 +4375,39 @@ def _check_equality_expression(self, node, child_data): return vals, const - def _check_ranged(self, node, child_data): + 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)})." + ) + + return [val], const + + 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 + + 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:]) @@ -4443,19 +4421,41 @@ def _check_ranged(self, node, child_data): return vals, const - def _sum_combinations(self, values_list): - sums = [] - for i in chain.from_iterable( - combinations(values_list, r) for r in range(2, len(values_list) + 1) - ): - sums.append(abs(sum(i))) + 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)) - return sums + # 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 node_type_method_map = { EXPR.EqualityExpression: _check_equality_expression, EXPR.InequalityExpression: _check_equality_expression, - EXPR.RangedExpression: _check_ranged, + EXPR.RangedExpression: _check_ranged_expression, EXPR.SumExpression: _check_sum_expression, EXPR.NPV_SumExpression: _check_sum_expression, EXPR.ProductExpression: _check_general_expr, @@ -4481,7 +4481,7 @@ def exitNode(self, node, data): """ Method to call when exiting node to check for potential issues. """ - # Return [node values], [constant + # Return [node values], constant # first check if the node is a leaf nodetype = type(node) From d486e852b6788956be401afee8d6061309b7dd0a Mon Sep 17 00:00:00 2001 From: Andrew Lee Date: Sun, 25 Aug 2024 16:03:31 -0400 Subject: [PATCH 16/20] Fixing false positives for linking constraints --- idaes/core/util/model_diagnostics.py | 45 ++++++++++--------- .../core/util/tests/test_model_diagnostics.py | 38 ++++++++++++++-- 2 files changed, 58 insertions(+), 25 deletions(-) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index 08099b34ce..bfc1a2565f 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -4301,14 +4301,6 @@ def __init__( self.canceling_terms = ComponentSet() self.mismatched_terms = ComponentSet() - def _check_base_type(self, node): - # [value], constant - if isinstance(node, VarData): - const = node.fixed - else: - const = True - return [value(node)], const - def _get_value_for_sum_subexpression(self, child_data): return sum(i for i in child_data[0]) @@ -4342,6 +4334,13 @@ def _perform_checks(self, node, child_data): # 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 @@ -4356,12 +4355,16 @@ def _check_equality_expression(self, node, child_data): mdata.append(child_data[i]) # Next, call the method to check the sum expression - vals, const = self._check_sum_expression(node, mdata) + vals, const, sum_expr = 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 - if node not in self.mismatched_terms and not any(d[1] for d in mdata): + # 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 t = [] @@ -4373,7 +4376,7 @@ def _check_equality_expression(self, node, child_data): if any(i <= self._sum_tol * max(t) for i in sums): self.canceling_terms.add(node) - return vals, const + return vals, const, False def _check_general_expr(self, node, child_data): const = self._perform_checks(node, child_data) @@ -4393,7 +4396,7 @@ def _check_general_expr(self, node, child_data): f"({str(node)})." ) - return [val], const + return [val], const, False def _check_other_expression(self, node, child_data): const = self._perform_checks(node, child_data) @@ -4405,11 +4408,11 @@ def _check_other_expression(self, node, child_data): newfunc = node.create_node_with_local_data(input_mag) # Evaluate new function and return the value along with check results - return [value(newfunc)], const + 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:]) + 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 @@ -4419,7 +4422,7 @@ def _check_ranged_expression(self, node, child_data): # Also for sign convention, we will negate the outer terms vals = lhs_vals + [-rhs_vals[1]] - return vals, const + return vals, const, False def _check_sum_expression(self, node, child_data): # Sum expressions need special handling @@ -4450,7 +4453,7 @@ def _check_sum_expression(self, node, child_data): if diff >= self._log_mm_tol: self.mismatched_terms.add(node) - return vals, const + return vals, const, True node_type_method_map = { EXPR.EqualityExpression: _check_equality_expression, @@ -4481,12 +4484,12 @@ def exitNode(self, node, data): """ Method to call when exiting node to check for potential issues. """ - # Return [node values], constant + # 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 + return [node], True, False node_func = self.node_type_method_map.get(nodetype, None) if node_func is not None: @@ -4495,7 +4498,7 @@ def exitNode(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 + return [1], True, False # Var or Param return self._check_base_type(node) @@ -4531,7 +4534,7 @@ def walk_expression(self, expr): self.mismatched_terms = ComponentSet() # Call parent walk_expression method - vals, const = super().walk_expression(expr) + 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 86bb4fd64e..0869d6bd66 100644 --- a/idaes/core/util/tests/test_model_diagnostics.py +++ b/idaes/core/util/tests/test_model_diagnostics.py @@ -1284,7 +1284,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: 2 Constraints with potential cancellation of terms" in cautions + assert "Caution: 1 Constraint with potential cancellation of terms" in cautions @pytest.mark.component def test_collect_numerical_cautions_jacobian(self): @@ -1470,7 +1470,7 @@ def test_report_numerical_issues_ok(self): ------------------------------------------------------------------------------------ 1 Cautions - Caution: 2 Constraints with potential cancellation of terms + Caution: 1 Constraint with potential cancellation of terms ------------------------------------------------------------------------------------ Suggested next steps: @@ -1482,7 +1482,7 @@ def test_report_numerical_issues_ok(self): ==================================================================================== """ - + print(stream.getvalue()) assert stream.getvalue() == expected @pytest.mark.component @@ -1555,7 +1555,7 @@ def test_report_numerical_issues(self, model): 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: 2 Constraints with potential cancellation of terms + Caution: 1 Constraint with potential cancellation of terms Caution: 1 extreme Jacobian Entry (<1.0E-04 or >1.0E+04) ------------------------------------------------------------------------------------ @@ -4900,3 +4900,33 @@ def test_canceling_equality_expr_canceling_sides(self): 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 From 1048b324bb2c422ac3a8fa9d2590d0f7977c3dcc Mon Sep 17 00:00:00 2001 From: Andrew Lee Date: Sun, 25 Aug 2024 16:10:47 -0400 Subject: [PATCH 17/20] Fixing unused variable warning --- idaes/core/util/model_diagnostics.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index bfc1a2565f..97bc305e07 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -4355,7 +4355,7 @@ def _check_equality_expression(self, node, child_data): mdata.append(child_data[i]) # Next, call the method to check the sum expression - vals, const, sum_expr = self._check_sum_expression(node, mdata) + 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() From d10a493c75b34093605f174ba94629647ad01c76 Mon Sep 17 00:00:00 2001 From: Andrew Lee Date: Fri, 6 Sep 2024 17:12:58 -0400 Subject: [PATCH 18/20] Fixing bug related to external functions with string arguments --- idaes/core/util/scaling.py | 12 ++++++++---- idaes/core/util/tests/test_scaling.py | 11 +++++++++++ 2 files changed, 19 insertions(+), 4 deletions(-) diff --git a/idaes/core/util/scaling.py b/idaes/core/util/scaling.py index 5e03a2a9b8..ccb305c041 100644 --- a/idaes/core/util/scaling.py +++ b/idaes/core/util/scaling.py @@ -1472,10 +1472,14 @@ def _get_nominal_value_expr_if(self, node, child_nominal_values): def _get_nominal_value_external_function(self, node, child_nominal_values): # First, need to get expected magnitudes of input terms, which may be sub-expressions - input_mag = [ - self._get_nominal_value_for_sum_subexpression(i) - for i in child_nominal_values - ] + input_mag = [] + for i in child_nominal_values: + if isinstance(i[0], str): + # Sometimes external functions might have string arguments + # Check here, and return the string if true + input_mag.append(i[0]) + else: + input_mag.append(self._get_nominal_value_for_sum_subexpression(i)) # Next, create a copy of the external function with expected magnitudes as inputs newfunc = node.create_node_with_local_data(input_mag) diff --git a/idaes/core/util/tests/test_scaling.py b/idaes/core/util/tests/test_scaling.py index e5ed63844d..f07ad129cd 100644 --- a/idaes/core/util/tests/test_scaling.py +++ b/idaes/core/util/tests/test_scaling.py @@ -39,6 +39,7 @@ CubicThermoExpressions, CubicType as CubicEoS, ) +from idaes.models.properties import iapws95 __author__ = "John Eslick, Tim Bartholomew" @@ -2230,6 +2231,16 @@ def test_constraint(self, m): expr=m.constraint.expr ) == [21, 0.5 ** (22 + 23 + 24)] + @pytest.mark.component + def test_external_function_w_string_argument(self): + m = pyo.ConcreteModel() + m.properties = iapws95.Iapws95ParameterBlock() + m.state = m.properties.build_state_block([0]) + + assert sc.NominalValueExtractionVisitor().walk_expression( + expr=m.state[0].temperature + ) == [pytest.approx(235.0, rel=1e-8)] + @pytest.fixture(scope="function") def m(): From 6d6fe2c930cc21f5d3dcbaff380f6e05e0bad2fc Mon Sep 17 00:00:00 2001 From: Andrew Lee Date: Tue, 10 Sep 2024 15:54:04 -0400 Subject: [PATCH 19/20] Addressing comments --- idaes/core/util/model_diagnostics.py | 74 ++++++++++++------- .../core/util/tests/test_model_diagnostics.py | 29 ++++++-- 2 files changed, 71 insertions(+), 32 deletions(-) diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index 97bc305e07..56f36087ed 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -1086,13 +1086,15 @@ def _collect_constraint_mismatches(self, descend_into=True): for c in self._model.component_data_objects( Constraint, descend_into=descend_into ): - _, mm, cc, k = walker.walk_expression(c.expr) + _, expr_mismatch, expr_cancellation, expr_constant = walker.walk_expression( + c.expr + ) - for i in mm: + for i in expr_mismatch: mismatch.append(f"{c.name}: {i}") - for i in cc: + for i in expr_cancellation: cancellation.append(f"{c.name}: {i}") - if k: + if expr_constant: constant.append(c.name) return mismatch, cancellation, constant @@ -4304,14 +4306,35 @@ def __init__( def _get_value_for_sum_subexpression(self, child_data): return sum(i for i in child_data[0]) - def _sum_combinations(self, values_list): - sums = [] - for i in chain.from_iterable( - combinations(values_list, r) for r in range(2, len(values_list) + 1) - ): - sums.append(abs(sum(i))) + 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)) - return sums + 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 @@ -4323,8 +4346,7 @@ def _perform_checks(self, node, 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 - sums = self._sum_combinations(d[0]) - if any(i <= self._sum_tol * max(d[0]) for i in sums): + if self._check_sum_cancellations(d[0]): self.canceling_terms.add(node) # Expression is not constant if any child is not constant @@ -4345,14 +4367,12 @@ 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 = [] - for i in [0, 1]: - if i == 0: - vals = [] - for j in child_data[i][0]: - vals.append(-j) - mdata.append((vals, child_data[i][1])) - else: - mdata.append(child_data[i]) + + 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) @@ -4367,13 +4387,11 @@ def _check_equality_expression(self, node, child_data): 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 - t = [] - for d in mdata: - t += d[0] + # Note: outer loop comes first in list comprehension + t = [i for d in mdata for i in d[0]] # Then check for cancellations - sums = self._sum_combinations(t) - if any(i <= self._sum_tol * max(t) for i in sums): + if self._check_sum_cancellations(t): self.canceling_terms.add(node) return vals, const, False @@ -4395,6 +4413,10 @@ def _check_general_expr(self, node, child_data): 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 diff --git a/idaes/core/util/tests/test_model_diagnostics.py b/idaes/core/util/tests/test_model_diagnostics.py index 0869d6bd66..d67380b780 100644 --- a/idaes/core/util/tests/test_model_diagnostics.py +++ b/idaes/core/util/tests/test_model_diagnostics.py @@ -68,6 +68,7 @@ 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 idaes.models.properties.modular_properties.eos.ceos_common import ( cubic_roots_available, @@ -106,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" @@ -1301,7 +1303,7 @@ def test_collect_numerical_cautions_jacobian(self): cautions = dt._collect_numerical_cautions() - assert len(cautions) == 5 + assert len(cautions) == 4 assert "Caution: 3 Variables with value close to zero (tol=1.0E-08)" in cautions assert ( "Caution: 3 Variables with extreme Jacobian values (<1.0E-04 or >1.0E+04)" @@ -1312,7 +1314,6 @@ def test_collect_numerical_cautions_jacobian(self): in cautions ) assert "Caution: 4 extreme Jacobian Entries (<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_assert_no_structural_warnings(self, model): @@ -1600,10 +1601,9 @@ def test_report_numerical_issues_jacobian(self): WARNING: 3 pairs of variables are parallel (to tolerance 1.0E-08) ------------------------------------------------------------------------------------ -5 Cautions +4 Cautions Caution: 3 Variables with value close to zero (tol=1.0E-08) - Caution: 1 Constraint with potential cancellation of terms Caution: 3 Variables with extreme Jacobian values (<1.0E-04 or >1.0E+04) Caution: 1 Constraint with extreme Jacobian values (<1.0E-04 or >1.0E+04) Caution: 4 extreme Jacobian Entries (<1.0E-04 or >1.0E+04) @@ -4188,7 +4188,7 @@ def test_sum_combinations(self): # excludes single term sums terms = [1, 2, 3, 4, 5] visitor = ConstraintTermAnalysisVisitor() - sums = visitor._sum_combinations(terms) + sums = [i for i in visitor._generate_sum_combinations(terms)] expected = [ 3, @@ -4218,7 +4218,7 @@ def test_sum_combinations(self): 14, # 4-term sums 15, # 5-term sum ] - + print(sums) assert sums == expected @pytest.mark.unit @@ -4930,3 +4930,20 @@ def test_linking_equality_expr_compound(self): 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 From 678caf874b88c21f626c4e4ded8bb37887ef97fa Mon Sep 17 00:00:00 2001 From: Andrew Lee Date: Wed, 11 Sep 2024 11:12:04 -0400 Subject: [PATCH 20/20] removing debugging prints --- idaes/core/util/tests/test_model_diagnostics.py | 5 ++--- 1 file changed, 2 insertions(+), 3 deletions(-) diff --git a/idaes/core/util/tests/test_model_diagnostics.py b/idaes/core/util/tests/test_model_diagnostics.py index d67380b780..f9bce71f55 100644 --- a/idaes/core/util/tests/test_model_diagnostics.py +++ b/idaes/core/util/tests/test_model_diagnostics.py @@ -1483,7 +1483,7 @@ def test_report_numerical_issues_ok(self): ==================================================================================== """ - print(stream.getvalue()) + assert stream.getvalue() == expected @pytest.mark.component @@ -4218,7 +4218,7 @@ def test_sum_combinations(self): 14, # 4-term sums 15, # 5-term sum ] - print(sums) + assert sums == expected @pytest.mark.unit @@ -4867,7 +4867,6 @@ def test_canceling_equality_expr_compound(self): assert vv == [0, pytest.approx(0, abs=1e-12)] assert len(mm) == 0 - print(cc) assert expr in cc assert len(cc) == 1 assert not k