Diagnostics tool for ill-posed constraints

[andrewlee94]

Fixes None

Summary/Motivation:

As part of working on the new scaling tools, I realised there are some simple checks we can do for detecting poorly-posed constraints that could cause scaling issues. This PR adds a new expression walker that looks for the following signs of poor scaling in constraints:

• sum expressions with terms which have significant differences in magnitude.
• sum expressions where terms cancel out (with a catch for cases of constant == sum()).

Changes proposed in this PR:

Legal Acknowledgement

By contributing to this software project, I agree to the following terms and conditions for my contribution:

1. I agree my contributions are submitted under the license terms described in the LICENSE.txt file at the top level of this directory.
2. I represent I am authorized to make the contributions and grant the license. If my employer has rights to intellectual property that includes these contributions, I represent that I have received permission to make contributions and grant the required license on behalf of that employer.

[codecov-commenter]

Codecov Report

Attention: Patch coverage is 84.81013% with 24 lines in your changes missing coverage. Please review.

Project coverage is 76.38%. Comparing base (2b41748) to head (1048b32).

Files	Patch %	Lines
idaes/core/util/model_diagnostics.py	84.81%	15 Missing and 9 partials ⚠️

Additional details and impacted files

@@            Coverage Diff             @@
##             main    #1454      +/-   ##
==========================================
+ Coverage   76.37%   76.38%   +0.01%     
==========================================
  Files         393      393              
  Lines       65085    65242     +157     
  Branches    14426    14466      +40     
==========================================
+ Hits        49708    49838     +130     
- Misses      12815    12834      +19     
- Partials     2562     2570       +8     

☔ View full report in Codecov by Sentry.
📢 Have feedback on the report? Share it here.

[dallan-keylogic]

@Robbybp , is there any possibility of combining a tool like this with block triangularization? My intuition is that a constraint in which catastrophic cancellation occurs (taking the difference between two large numbers to get a small number) might show up differently in the block structure than one in which it does not (adding a small number to a large number to get another large number), even though individual constraint might be identical.

[Robbybp]

@dallan-keylogic My only initial thought is that a "catastrophic cancellation" has different implications depending on where it appears appears in the block triangular decomposition. If it is in a diagonal block, it could be a problem (as we're relying on that term for nonsingularity). Otherwise, it is fine, as the entry could be zero and the matrix's singularity would not change.

[dallan-keylogic]

Okay, a very basic example to see the issues at play. We're not concerned when we add a small number to a big number to get another big number:

This system has the matrix representation

However, the situation is different when we take the difference of two big numbers to get a small number:

This system has the matrix representation

Looking at these two systems, there doesn't seem to be any structural difference between them. The variable b is chosen to be a pivot in both cases. However, when we define a scaled variable

we get a different picture.

The first system then becomes

with a matrix representation

whose system matrix is well-conditioned with determinant -1.

The second system, on the other hand, becomes

with a matrix representation

whose system matrix is ill-conditioned with determinant -epsilon. It can be shown that the condition number blows up in the limit as epsilon goes to zero.

[Robbybp]

@dallan-keylogic I would say that scaling b by 1/eps is not the right choice here. A block triangularization on the scaled system could definitely identify that. We could potentially identify that b, as a column-singleton, should not be scaled up.

[dallan-keylogic]

@Robbybp Whether or not we should scale by 1/eps requires a judgement call about what precision we need b. If having b to the same absolute precision as a and c is sufficient, we should not scale it. If we need b to the same relative precision as a and c, we must scale it. We should scale it when, for example, we're then going to divide it by a small number as part of a finite difference approximation to a derivative (note, @andrewlee94 , this sort of action is taken in the ControlVolume1D).

[jsiirola]

Overall looks pretty reasonable. Some questions and comments (that should probably be addressed)

[jsiirola]

I would be happier if the "collect" routine didn't do formatting (conversion to a string).

[andrewlee94]

This one would take some more work and require changing other collect methods as well. I think that might be better as a separate issue/PR.

[jsiirola]

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)():

[andrewlee94]

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

[jsiirola]

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?

[dallan-keylogic]

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.

[andrewlee94]

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.

[dallan-keylogic]

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.

[andrewlee94]

@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)?

[jsiirola]

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

[andrewlee94]

I had forgotten to make that correction - thank you.

[jsiirola]

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

[andrewlee94]

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

[jsiirola]

Almost all of your nonlinear handlers could be replaced by a single callback:

def _check_general_expr(self, node, child_data):
    mismatch, cancelling, const = self._perform_checks(node, child_data)
    val = node._apply_operation(list(map(self._get_value_for_sum_subexpression, child_data)))
    return [val], mismatch, cancelling, const

[andrewlee94]

Thank you - I did not know I could do that.

[dallan-keylogic]

Something that occurred to me is that it might be useful for the user to call this not on the entire model, but on a subset of constraints (perhaps identified by the SVD toolbox, perhaps identified some other way). That probably doesn't need to be present during the initial rollout, but is something to keep in mind.

[dallan-keylogic]

Do we have a test case for a standard linking constraint like this with the variables unfixed? We want to make sure that arc constraints aren't going to get lumped in here.

[andrewlee94]

As it stands, Arc constraints will get lumped in here, however from a strict sense they should fall in here (as they can be trivially eliminated). However, from a structural stand point we probably want to ignore these.

What do we want to allow? One option would be the special case of a=b where a and b are leaf nodes (i.e. not expressions themselves).

[andrewlee94]

@dallan-keylogic I believe I have fixed this - the cancellation checks are only run if at least one of the sides in an (in)equality is a sum expression. This means that constraints with the form a == b and a == f(b) will not be flagged as cancelling.

[andrewlee94]

Also, regarding you initial comment, this a tool to do that could be written fairly easily as it would just iterate over a ComponentSet containing the Constraint of interest and call this walker for each Constraint. This walker is not part of the DiagnosticsToolbox class, so it can easily be called from other places too (such as the SVDToolbox) if we can come up with a way to identify a set of Constraints` of interest.

[dallan-keylogic]

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.

[mrmundt]

I did not review the tests cases yet because I do have some suggestions for changes in the main file that may require changes in test cases.

[mrmundt]

2 is a magic number! Why 2?

[andrewlee94]

We are looking for any combination of terms which cancel, thus the minimum number of terms to consider is 2 (a single term cannot cancel with itself). I can add a comment.

[dallan-keylogic]

What this block of code does is go through combinations of terms in an expression first in groups of 2, then in groups of 3, all the way up to groups of len(values_list). So if we have the sum expression a+b+c+d, we first iterate through (a, b), (a, c), (a, d), (b, c), (b, d), (c, d), then iterate through (a, b, c,), (a, b, d), (a, c, d), (b, c, d), then (a,b,c,d).

However, the number of terms you're checking grows exponentially in expression size. In particular, if len(values_list) == m, you'll be checking 2 ** m - m -1 terms. I expect that this sort of check will take an extremely long time on any model with Expressions of any significant length, much less an extreme chonker like eNRTL.

[andrewlee94]

We can probably make this more efficient by:

1. Stripping any term with a value of 0
2. Breaking the loop at the first failure - we do not count how many cancellations there are, just if there is at least 1.

[andrewlee94]

@mrmundt I think I have addressed all of your comment. If you have time to do another review, it would be appreciated.

[andrewlee94]

@mrmundt Thank you - I fixed the left over print statements.