Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Diagnostics tool for ill-posed constraints #1454

Merged
merged 62 commits into from
Nov 5, 2024
Merged
Show file tree
Hide file tree
Changes from 1 commit
Commits
Show all changes
62 commits
Select commit Hold shift + click to select a range
37e0fae
Working on cancellation detection
Jul 22, 2024
1664593
Merge branch 'main' of https://github.com/IDAES/idaes-pse into cons_m…
Jul 22, 2024
4d48741
Fixing typo
Jul 22, 2024
a1e91d2
Merge branch 'main' of https://github.com/IDAES/idaes-pse into cons_m…
Aug 6, 2024
33a4435
Improving efficiency
Aug 6, 2024
747eb7d
More efficiency and testing
Aug 6, 2024
707de1e
Reducing duplicated code in tests
Aug 9, 2024
3a6400b
Finishing testing
Aug 9, 2024
c3e315c
Adding additioanl tests and clean up
Aug 9, 2024
2ab55ca
Running pylint
Aug 9, 2024
2af9478
Adding constraint walker to diagnostics toolbox
Aug 9, 2024
b9d12ea
Merge branch 'main' into cons_mag_diagnostics
Aug 23, 2024
162f631
Addressing simple review comments
Aug 23, 2024
90d2cd3
Accumulate node rather than str(node)
Aug 23, 2024
152b044
Clean up some tests
Aug 23, 2024
81c1835
Supporting ranged expressions
Aug 23, 2024
849647c
Merge branch 'main' into cons_mag_diagnostics
Aug 23, 2024
bff77d5
Apply suggestions from code review
andrewlee94 Aug 23, 2024
54e87c4
Merge branch 'main' of https://github.com/IDAES/idaes-pse into cons_m…
Aug 25, 2024
11d3f33
Merge branch 'cons_mag_diagnostics' of https://github.com/andrewlee94…
Aug 25, 2024
85d1956
Regorganising and fixing pylint issues
Aug 25, 2024
d486e85
Fixing false positives for linking constraints
Aug 25, 2024
1048b32
Fixing unused variable warning
Aug 25, 2024
d10a493
Fixing bug related to external functions with string arguments
Sep 6, 2024
00a3434
Merge branch 'main' of https://github.com/IDAES/idaes-pse into cons_m…
Sep 10, 2024
d27de22
Merge branch 'old_scaling_bug' into cons_mag_diagnostics
Sep 10, 2024
6d6fe2c
Addressing comments
Sep 10, 2024
9dbd521
Merge branch 'main' into cons_mag_diagnostics
Sep 11, 2024
678caf8
removing debugging prints
Sep 11, 2024
47a58f4
Merge branch 'main' of https://github.com/IDAES/idaes-pse into cons_m…
Sep 23, 2024
1b5d1d2
Addressing more comments
Sep 23, 2024
18a6ff1
Merge branch 'cons_mag_diagnostics' of https://github.com/andrewlee94…
Sep 23, 2024
ef2e228
Debugging issue with nested external functions
Sep 24, 2024
70636cb
Adding todo
Sep 24, 2024
5d5da12
Relax test tolerance due to Windows failures
Sep 24, 2024
c173984
Adding catches for string arguments/values
Sep 30, 2024
7c03467
Test for nested external functions
Sep 30, 2024
21cca53
Fixing bug in walker logic for string arguments to external funcitons
Oct 8, 2024
ad6dd99
Merge branch 'main' into cons_mag_diagnostics
Oct 9, 2024
a3cde35
Merge branch 'main' into cons_mag_diagnostics
Oct 25, 2024
ddc8fb8
Address double counting of constraints
Oct 25, 2024
bd331e8
Fixing typo
Oct 25, 2024
b3098e1
Bumping Pyomo version
Oct 25, 2024
60e71cb
Adding zero tolerance
Oct 25, 2024
aebb118
Adding test for zero tolerance
Oct 25, 2024
9956286
Fixing typo
Oct 28, 2024
8cf3bbe
Logic to skip scaled equality constraints
Oct 28, 2024
34a2b58
Merge branch 'main' into cons_mag_diagnostics_2
Oct 28, 2024
93ef9d8
Handling mole fraction type constraints
Oct 28, 2024
3284d3e
Merge branch 'main' into cons_mag_diagnostics
Oct 29, 2024
2d5ec93
Catching negations
Oct 29, 2024
5ffec52
Fixing typos
Oct 29, 2024
161b15b
Limit of cancelling combinations
Oct 30, 2024
7deb640
Fixing typos
Oct 30, 2024
86bd61c
Displaying cancelling terms to user
Oct 30, 2024
9b13cb1
More detail of mismatches
Oct 30, 2024
6577481
Adding hooks and warnings about new config arguments
Oct 30, 2024
94fe8c0
Expression aware to string method
Nov 1, 2024
360bb39
Moving expression writer to util/misc
Nov 1, 2024
1fb2743
Typos and comments
Nov 1, 2024
628f445
Using StrEnum
Nov 5, 2024
1a4a1ee
Merge branch 'main' into cons_mag_diagnostics
Nov 5, 2024
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
Next Next commit
Working on cancellation detection
  • Loading branch information
Andrew Lee committed Jul 22, 2024
commit 37e0faee6f708a507f3d273498a6efb28b4352bd
167 changes: 166 additions & 1 deletion idaes/core/util/model_diagnostics.py
Original file line number Diff line number Diff line change
Expand Up @@ -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
Expand Down Expand Up @@ -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,
)
Expand All @@ -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 (
Expand Down Expand Up @@ -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
andrewlee94 marked this conversation as resolved.
Show resolved Hide resolved

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))
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think I understand what you are doing here, but wouldn't this complain loudly about the objective for every parameter estimation problem (MSE) if the problem actually solved?

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

One of my main concerns with this tool is its inability to distinguish between true cases of catastrophic cancellation and benign situations, where, for example, we have a heat of adsorption term in a column which will be close to except near the breakthrough point.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is there an easy way to test for that however? Note that this tool is only used for Cautions in the main toolbox, with the implication that these might be issues (but not guaranteeing that).

I.e., this is intended to be a simple check to find potential issues, but the user will have to look into them all further to decide if they are critical or not.

Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I'm not sure there's an easy way to do it without assuming a well-scaled system. In a well-scaled system, it would show up in the SVD with extremely large or extremely small singular values.

However, once you get to a high enough noise/signal ratio, you start to question about whether including the tool in report_numerical_issues is worthwhile or whether it is distracting the user from more fruitful avenues of diagnostics. I suppose we can just pull it from report_numerical_issues without a deprecation cycle if it proves to not be useful, so we can try it out and see how users find it.

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

@jsiirola This walker was written specifically with Constraints in mind, so I had not considered that. The check for mismatched terms would make sense for an Objective however, so this could be extended to handle them as well. However, how would the expression walker know if it was dealing with an Objective or a Constraint (the input argument is the expression, not the component)?


# 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)
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Don't you need to negate the values for the second argument in the relational expression before treating it like a sum?

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I had forgotten to make that correction - thank you.


# 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,
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Don't you need special handling for Ranged (like Equality / Inequality)? I would probably define a

def _check_ranged(self, node, child_data):
    lhs_vals, lhs_mismatch, lhs_cancelling, lhs_const = self._check_equality(node, child_data[:2])
    rhs_vals, rhs_mismatch, rhs_cancelling, rhs_const = self._check_equality(node, child_data[1:])
    # Then merge the results and return

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think I've fixed this - it would be good if you could double check my logic however.

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()
):
Comment on lines +4833 to +4836
Copy link
Contributor

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I have stumbled across a syntax that is a bit more concise. Not sure if you want to use it, though:

if getattr(node, "is_named_expression_type", bool)():

Copy link
Member Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think I'll keep the current form as it is a little easier to understand what it is doing.

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)}"
)
Loading