diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index 6400d219a3..5194ee0064 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -20,7 +20,9 @@ from operator import itemgetter import sys from inspect import signature +import math from math import log +from typing import List import numpy as np from scipy.linalg import svd @@ -29,6 +31,7 @@ from pyomo.environ import ( Binary, + Integers, Block, check_optimal_termination, ConcreteModel, @@ -40,9 +43,21 @@ value, Var, ) +from pyomo.core.expr.numeric_expr import ( + DivisionExpression, + NPV_DivisionExpression, + PowExpression, + NPV_PowExpression, + UnaryFunctionExpression, + NPV_UnaryFunctionExpression, + NumericExpression, +) from pyomo.core.base.block import _BlockData -from pyomo.core.base.var import _VarData +from pyomo.core.base.var import _GeneralVarData, _VarData from pyomo.core.base.constraint import _ConstraintData +from pyomo.repn.standard_repn import ( # pylint: disable=no-name-in-module + generate_standard_repn, +) from pyomo.common.collections import ComponentSet from pyomo.common.config import ( ConfigDict, @@ -52,9 +67,10 @@ ) from pyomo.util.check_units import identify_inconsistent_units from pyomo.contrib.incidence_analysis import IncidenceGraphInterface -from pyomo.core.expr.visitor import identify_variables +from pyomo.core.expr.visitor import identify_variables, StreamBasedExpressionVisitor from pyomo.contrib.pynumero.interfaces.pyomo_nlp import PyomoNLP from pyomo.contrib.pynumero.asl import AmplInterface +from pyomo.contrib.fbbt.fbbt import compute_bounds_on_expr from pyomo.common.deprecation import deprecation_warning from pyomo.common.errors import PyomoException @@ -251,6 +267,14 @@ def svd_sparse(jacobian, number_singular_values): description="Tolerance for raising a warning for small Jacobian values.", ), ) +CONFIG.declare( + "warn_for_evaluation_error_at_bounds", + ConfigValue( + default=True, + domain=bool, + description="If False, warnings will not be generated for things like log(x) with x >= 0", + ), +) SVDCONFIG = ConfigDict() @@ -936,6 +960,13 @@ def _collect_structural_warnings(self): if any(len(x) > 0 for x in [oc_var, oc_con]): next_steps.append(self.display_overconstrained_set.__name__ + "()") + eval_warnings = self._collect_potential_eval_errors() + if len(eval_warnings) > 0: + warnings.append( + f"WARNING: Found {len(eval_warnings)} potential evaluation errors." + ) + next_steps.append(self.display_potential_evaluation_errors.__name__ + "()") + return warnings, next_steps def _collect_structural_cautions(self): @@ -1303,6 +1334,51 @@ def report_numerical_issues(self, stream=None): footer="=", ) + def _collect_potential_eval_errors(self) -> List[str]: + warnings = list() + for con in self.model.component_data_objects( + Constraint, active=True, descend_into=True + ): + walker = _EvalErrorWalker(self.config) + con_warnings = walker.walk_expression(con.body) + for msg in con_warnings: + msg = f"{con.name}: " + msg + warnings.append(msg) + for obj in self.model.component_data_objects( + Objective, active=True, descend_into=True + ): + walker = _EvalErrorWalker(self.config) + obj_warnings = walker.walk_expression(obj.expr) + for msg in obj_warnings: + msg = f"{obj.name}: " + msg + warnings.append(msg) + + return warnings + + def display_potential_evaluation_errors(self, stream=None): + """ + Prints constraints that may be prone to evaluation errors + (e.g., log of a negative number) based on variable bounds. + + Args: + stream: an I/O object to write the output to (default = stdout) + + Returns: + None + """ + if stream is None: + stream = sys.stdout + + warnings = self._collect_potential_eval_errors() + _write_report_section( + stream=stream, + lines_list=warnings, + title=f"{len(warnings)} WARNINGS", + line_if_empty="No warnings found!", + header="=", + footer="=", + ) + @document_kwargs_from_configdict(SVDCONFIG) def prepare_svd_toolbox(self, **kwargs): """ @@ -1623,6 +1699,170 @@ def display_variables_in_constraint(self, constraint, stream=None): ) +def _get_bounds_with_inf(node: NumericExpression): + lb, ub = compute_bounds_on_expr(node) + if lb is None: + lb = -math.inf + if ub is None: + ub = math.inf + return lb, ub + + +def _check_eval_error_division( + node: NumericExpression, warn_list: List[str], config: ConfigDict +): + lb, ub = _get_bounds_with_inf(node.args[1]) + if (config.warn_for_evaluation_error_at_bounds and (lb <= 0 <= ub)) or ( + lb < 0 < ub + ): + msg = f"Potential division by 0 in {node}; Denominator bounds are ({lb}, {ub})" + warn_list.append(msg) + + +def _check_eval_error_pow( + node: NumericExpression, warn_list: List[str], config: ConfigDict +): + arg1, arg2 = node.args + lb1, ub1 = _get_bounds_with_inf(arg1) + lb2, ub2 = _get_bounds_with_inf(arg2) + + integer_domains = ComponentSet([Binary, Integers]) + + integer_exponent = False + # if the exponent is an integer, there should not be any evaluation errors + if isinstance(arg2, _GeneralVarData) and arg2.domain in integer_domains: + # The exponent is an integer variable + # check if the base can be zero + integer_exponent = True + if lb2 == ub2 and lb2 == round(lb2): + # The exponent is fixed to an integer + integer_exponent = True + repn = generate_standard_repn(arg2, quadratic=True) + if ( + repn.nonlinear_expr is None + and repn.constant == round(repn.constant) + and all(i.domain in integer_domains for i in repn.linear_vars) + and all(i[0].domain in integer_domains for i in repn.quadratic_vars) + and all(i[1].domain in integer_domains for i in repn.quadratic_vars) + and all(i == round(i) for i in repn.linear_coefs) + and all(i == round(i) for i in repn.quadratic_coefs) + ): + # The exponent is a linear or quadratic expression containing + # only integer variables with integer coefficients + integer_exponent = True + + if integer_exponent and ( + (lb1 > 0 or ub1 < 0) + or (not config.warn_for_evaluation_error_at_bounds and (lb1 >= 0 or ub1 <= 0)) + ): + # life is good; the exponent is an integer and the base is nonzero + return None + elif integer_exponent and lb2 >= 0: + # life is good; the exponent is a nonnegative integer + return None + + # if the base is positive, there should not be any evaluation errors + if lb1 > 0 or (not config.warn_for_evaluation_error_at_bounds and lb1 >= 0): + return None + if lb1 >= 0 and lb2 >= 0: + return None + + msg = f"Potential evaluation error in {node}; " + msg += f"base bounds are ({lb1}, {ub1}); " + msg += f"exponent bounds are ({lb2}, {ub2})" + warn_list.append(msg) + + +def _check_eval_error_log( + node: NumericExpression, warn_list: List[str], config: ConfigDict +): + lb, ub = _get_bounds_with_inf(node.args[0]) + if (config.warn_for_evaluation_error_at_bounds and lb <= 0) or lb < 0: + msg = f"Potential log of a non-positive number in {node}; Argument bounds are ({lb}, {ub})" + warn_list.append(msg) + + +def _check_eval_error_tan( + node: NumericExpression, warn_list: List[str], config: ConfigDict +): + lb, ub = _get_bounds_with_inf(node) + if not (math.isfinite(lb) and math.isfinite(ub)): + msg = f"{node} may evaluate to -inf or inf; Argument bounds are {_get_bounds_with_inf(node.args[0])}" + warn_list.append(msg) + + +def _check_eval_error_asin( + node: NumericExpression, warn_list: List[str], config: ConfigDict +): + lb, ub = _get_bounds_with_inf(node.args[0]) + if lb < -1 or ub > 1: + msg = f"Potential evaluation of asin outside [-1, 1] in {node}; Argument bounds are ({lb}, {ub})" + warn_list.append(msg) + + +def _check_eval_error_acos( + node: NumericExpression, warn_list: List[str], config: ConfigDict +): + lb, ub = _get_bounds_with_inf(node.args[0]) + if lb < -1 or ub > 1: + msg = f"Potential evaluation of acos outside [-1, 1] in {node}; Argument bounds are ({lb}, {ub})" + warn_list.append(msg) + + +def _check_eval_error_sqrt( + node: NumericExpression, warn_list: List[str], config: ConfigDict +): + lb, ub = _get_bounds_with_inf(node.args[0]) + if lb < 0: + msg = f"Potential square root of a negative number in {node}; Argument bounds are ({lb}, {ub})" + warn_list.append(msg) + + +_unary_eval_err_handler = dict() +_unary_eval_err_handler["log"] = _check_eval_error_log +_unary_eval_err_handler["log10"] = _check_eval_error_log +_unary_eval_err_handler["tan"] = _check_eval_error_tan +_unary_eval_err_handler["asin"] = _check_eval_error_asin +_unary_eval_err_handler["acos"] = _check_eval_error_acos +_unary_eval_err_handler["sqrt"] = _check_eval_error_sqrt + + +def _check_eval_error_unary( + node: NumericExpression, warn_list: List[str], config: ConfigDict +): + if node.getname() in _unary_eval_err_handler: + _unary_eval_err_handler[node.getname()](node, warn_list, config) + + +_eval_err_handler = dict() +_eval_err_handler[DivisionExpression] = _check_eval_error_division +_eval_err_handler[NPV_DivisionExpression] = _check_eval_error_division +_eval_err_handler[PowExpression] = _check_eval_error_pow +_eval_err_handler[NPV_PowExpression] = _check_eval_error_pow +_eval_err_handler[UnaryFunctionExpression] = _check_eval_error_unary +_eval_err_handler[NPV_UnaryFunctionExpression] = _check_eval_error_unary + + +class _EvalErrorWalker(StreamBasedExpressionVisitor): + def __init__(self, config: ConfigDict): + super().__init__() + self._warn_list = list() + self._config = config + + def exitNode(self, node, data): + """ + callback to be called as the visitor moves from the leaf + nodes back to the root node. + + Args: + node: a pyomo expression node + data: not used in this walker + """ + if type(node) in _eval_err_handler: + _eval_err_handler[type(node)](node, self._warn_list, self._config) + return self._warn_list + + # TODO: Rename and redirect once old DegeneracyHunter is removed @document_kwargs_from_configdict(DHCONFIG) class DegeneracyHunter2: diff --git a/idaes/core/util/tests/test_model_diagnostics.py b/idaes/core/util/tests/test_model_diagnostics.py index e91a69160b..fe59206ad9 100644 --- a/idaes/core/util/tests/test_model_diagnostics.py +++ b/idaes/core/util/tests/test_model_diagnostics.py @@ -14,6 +14,7 @@ This module contains model diagnostic utility functions for use in IDAES (Pyomo) models. """ from io import StringIO +import math import numpy as np import pytest @@ -23,6 +24,10 @@ Constraint, Expression, log, + tan, + asin, + acos, + sqrt, Objective, Set, SolverFactory, @@ -31,6 +36,8 @@ units, value, Var, + Param, + Integers, ) from pyomo.common.collections import ComponentSet from pyomo.contrib.pynumero.asl import AmplInterface @@ -41,6 +48,7 @@ from idaes.core.solvers import get_solver from idaes.core import FlowsheetBlock from idaes.core.util.testing import PhysicalParameterTestBlock +from unittest import TestCase # TODO: Add pyomo.dae test case """ @@ -2579,3 +2587,294 @@ def test_ipopt_solve_halt_on_error(capsys): captured = capsys.readouterr() assert "c: can't evaluate log(-5)." in captured.out + + +class TestEvalErrorDetection(TestCase): + @pytest.mark.unit + def test_div(self): + m = ConcreteModel() + m.x = Var(bounds=(1, None)) + m.y = Var() + m.c = Constraint(expr=m.y == 1 / m.x) + dtb = DiagnosticsToolbox(m) + warnings = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 0) + + m.x.setlb(0) + warnings = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 1) + w = warnings[0] + self.assertEqual( + w, "c: Potential division by 0 in 1/x; Denominator bounds are (0, inf)" + ) + + dtb.config.warn_for_evaluation_error_at_bounds = False + warnings = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 0) + + m.x.setlb(-1) + warnings = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 1) + w = warnings[0] + self.assertEqual( + w, "c: Potential division by 0 in 1/x; Denominator bounds are (-1, inf)" + ) + + @pytest.mark.unit + def test_pow1(self): + m = ConcreteModel() + m.x = Var(bounds=(None, None)) + m.y = Var() + m.p = Param(initialize=2, mutable=True) + m.c = Constraint(expr=m.y == m.x**m.p) + dtb = DiagnosticsToolbox(m) + warnings = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 0) + + m.p.value = 2.5 + warnings = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 1) + w = warnings[0] + self.assertEqual( + w, + "c: Potential evaluation error in x**p; base bounds are (-inf, inf); exponent bounds are (2.5, 2.5)", + ) + + m.x.setlb(1) + warnings = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 0) + + @pytest.mark.unit + def test_pow2(self): + m = ConcreteModel() + m.x = Var(bounds=(1, None)) + m.y = Var() + m.p = Var(domain=Integers) + m.c = Constraint(expr=m.y == m.x**m.p) + dtb = DiagnosticsToolbox(m) + warnings = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 0) + + m.x.setlb(None) + warnings = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 1) + w = warnings[0] + self.assertEqual( + w, + "c: Potential evaluation error in x**p; base bounds are (-inf, inf); exponent bounds are (-inf, inf)", + ) + + @pytest.mark.unit + def test_pow3(self): + m = ConcreteModel() + m.x = Var(bounds=(0, None)) + m.y = Var() + m.p = Var(bounds=(0, None)) + m.c = Constraint(expr=m.y == m.x**m.p) + dtb = DiagnosticsToolbox(m) + warnings = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 0) + + @pytest.mark.unit + def test_pow4(self): + m = ConcreteModel() + m.x = Var(bounds=(0, None)) + m.y = Var() + m.c = Constraint(expr=m.y == m.x ** (-2)) + dtb = DiagnosticsToolbox(m) + warnings = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 1) + w = warnings[0] + self.assertEqual( + w, + "c: Potential evaluation error in x**-2; base bounds are (0, inf); exponent bounds are (-2, -2)", + ) + + dtb.config.warn_for_evaluation_error_at_bounds = False + warnings = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 0) + + m.x.setlb(-1) + warnings = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 1) + w = warnings[0] + self.assertEqual( + w, + "c: Potential evaluation error in x**-2; base bounds are (-1, inf); exponent bounds are (-2, -2)", + ) + + @pytest.mark.unit + def test_pow5(self): + m = ConcreteModel() + m.x = Var(bounds=(0, None)) + m.y = Var() + m.c = Constraint(expr=m.y == m.x ** (-2.5)) + dtb = DiagnosticsToolbox(m) + warnings = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 1) + w = warnings[0] + self.assertEqual( + w, + "c: Potential evaluation error in x**-2.5; base bounds are (0, inf); exponent bounds are (-2.5, -2.5)", + ) + + dtb.config.warn_for_evaluation_error_at_bounds = False + warnings = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 0) + + m.x.setlb(-1) + warnings = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 1) + w = warnings[0] + self.assertEqual( + w, + "c: Potential evaluation error in x**-2.5; base bounds are (-1, inf); exponent bounds are (-2.5, -2.5)", + ) + + @pytest.mark.unit + def test_log(self): + m = ConcreteModel() + m.x = Var(bounds=(1, None)) + m.y = Var() + m.c = Constraint(expr=m.y == log(m.x)) + dtb = DiagnosticsToolbox(m) + warnings = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 0) + + m.x.setlb(0) + warnings = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 1) + w = warnings[0] + self.assertEqual( + w, + "c: Potential log of a non-positive number in log(x); Argument bounds are (0, inf)", + ) + + dtb.config.warn_for_evaluation_error_at_bounds = False + warnings = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 0) + + m.x.setlb(-1) + warnings = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 1) + w = warnings[0] + self.assertEqual( + w, + "c: Potential log of a non-positive number in log(x); Argument bounds are (-1, inf)", + ) + + @pytest.mark.unit + def test_tan(self): + m = ConcreteModel() + m.x = Var(bounds=(-math.pi / 4, math.pi / 4)) + m.y = Var() + m.c = Constraint(expr=m.y == tan(m.x)) + dtb = DiagnosticsToolbox(m) + warnings = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 0) + + m.x.setlb(-math.pi) + warnings = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 1) + w = warnings[0] + self.assertEqual( + w, + "c: tan(x) may evaluate to -inf or inf; Argument bounds are (-3.141592653589793, 0.7853981633974483)", + ) + + @pytest.mark.unit + def test_asin(self): + m = ConcreteModel() + m.x = Var(bounds=(-0.5, 0.5)) + m.y = Var() + m.c = Constraint(expr=m.y == asin(m.x)) + dtb = DiagnosticsToolbox(m) + warnings = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 0) + + m.x.setlb(None) + warnings = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 1) + w = warnings[0] + self.assertEqual( + w, + "c: Potential evaluation of asin outside [-1, 1] in asin(x); Argument bounds are (-inf, 0.5)", + ) + + m.x.setlb(-0.5) + m.x.setub(None) + warnings = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 1) + w = warnings[0] + self.assertEqual( + w, + "c: Potential evaluation of asin outside [-1, 1] in asin(x); Argument bounds are (-0.5, inf)", + ) + + @pytest.mark.unit + def test_acos(self): + m = ConcreteModel() + m.x = Var(bounds=(-0.5, 0.5)) + m.y = Var() + m.c = Constraint(expr=m.y == acos(m.x)) + dtb = DiagnosticsToolbox(m) + warnings = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 0) + + m.x.setlb(None) + warnings = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 1) + w = warnings[0] + self.assertEqual( + w, + "c: Potential evaluation of acos outside [-1, 1] in acos(x); Argument bounds are (-inf, 0.5)", + ) + + m.x.setlb(-0.5) + m.x.setub(None) + warnings = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 1) + w = warnings[0] + self.assertEqual( + w, + "c: Potential evaluation of acos outside [-1, 1] in acos(x); Argument bounds are (-0.5, inf)", + ) + + @pytest.mark.unit + def test_sqrt(self): + m = ConcreteModel() + m.x = Var(bounds=(1, None)) + m.y = Var() + m.c = Constraint(expr=m.y == sqrt(m.x)) + dtb = DiagnosticsToolbox(m) + warnings = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 0) + + m.x.setlb(-1) + warnings = dtb._collect_potential_eval_errors() + self.assertEqual(len(warnings), 1) + w = warnings[0] + self.assertEqual( + w, + "c: Potential square root of a negative number in sqrt(x); Argument bounds are (-1, inf)", + ) + + @pytest.mark.unit + def test_display(self): + stream = StringIO() + m = ConcreteModel() + m.x = Var() + m.y = Var() + m.obj = Objective(expr=m.x**2 + m.y**2.5) + m.c1 = Constraint(expr=m.y >= log(m.x)) + m.c2 = Constraint(expr=m.y >= (m.x - 1) ** 2.5) + m.c3 = Constraint(expr=m.x - 1 >= 0) + dtb = DiagnosticsToolbox(m) + dtb.display_potential_evaluation_errors(stream=stream) + expected = "====================================================================================\n3 WARNINGS\n\n c1: Potential log of a non-positive number in log(x); Argument bounds are (-inf, inf)\n c2: Potential evaluation error in (x - 1)**2.5; base bounds are (-inf, inf); exponent bounds are (2.5, 2.5)\n obj: Potential evaluation error in y**2.5; base bounds are (-inf, inf); exponent bounds are (2.5, 2.5)\n\n====================================================================================\n" + got = stream.getvalue() + exp_list = expected.split("\n") + got_list = got.split("\n") + self.assertEqual(len(exp_list), len(got_list)) + for _exp, _got in zip(exp_list, got_list): + self.assertEqual(_exp, _got)