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.

Sign up for GitHub

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

Jump to bottom

Add parallel constraint/variable check to report_numerical_issues #1385

Merged
merged 28 commits into from
May 3, 2024
Merged

Add parallel constraint/variable check to report_numerical_issues #1385

Show file tree
Hide file tree
Changes from 27 commits
Commits
Show all changes
28 commits
Select commit Hold shift + click to select a range
ed3eb3a
update imports of native_types and pyomo_constant_types (which is dep…
Robbybp Mar 6, 2024
842650b
Adding next batch of diagnostics tests
andrewlee94 Mar 8, 2024
6d09b5c
Next batch of daignsotics tests
andrewlee94 Mar 10, 2024
1389783
Merge branch 'consttype-import' of https://github.com/Robbybp/idaes-p…
andrewlee94 Mar 11, 2024
01e010c
Work around for ASL issue
andrewlee94 Mar 14, 2024
cc9c240
Adding more diagnostics checks
andrewlee94 Mar 14, 2024
2189347
Last unit model diagnostics tests
andrewlee94 Mar 14, 2024
7755690
Merge branch 'main' into diagnostics_testing
andrewlee94 Mar 14, 2024
8a9a7ee
Fixing typo
andrewlee94 Mar 14, 2024
92b54f4
Merge branch 'diagnostics_testing' of https://github.com/andrewlee94/…
andrewlee94 Mar 14, 2024
1b91594
Improving fix for ASL issue
andrewlee94 Mar 15, 2024
7f4a976
Better implementation of fix
andrewlee94 Mar 15, 2024
5af6d99
Merge branch 'main' of https://github.com/IDAES/idaes-pse into diagno…
andrewlee94 Mar 15, 2024
448081a
Fixing pylint and Python 3.8 failures
andrewlee94 Mar 15, 2024
2d8afda
Removing old implementation of workaround
andrewlee94 Mar 15, 2024
a4a84dd
Fixing noisy test
andrewlee94 Mar 15, 2024
4cf5fb2
Moving registration of external functions
andrewlee94 Mar 18, 2024
55258cf
Reverting to version that works on Windows
andrewlee94 Mar 18, 2024
cd47ae2
Trying another way to get binary files
andrewlee94 Mar 18, 2024
dc7f8cc
add parallel variable/constraint checks to report_numerical_issues; c…
Robbybp Mar 19, 2024
6c5a6c4
update tests
Robbybp Mar 23, 2024
0a1894a
Merge branch 'main' of https://github.com/idaes/idaes-pse into diagno…
Robbybp Mar 23, 2024
06c239d
Merge branch 'main' of https://github.com/idaes/idaes-pse into diagno…
Robbybp Apr 10, 2024
ad00bf9
tighten parallel_component_tolerance to 1e-8
Robbybp Apr 10, 2024
69ba71f
adjust model to make parallel variable test less sensitive
Robbybp Apr 10, 2024
8960343
update test to reflect new default tolerance of 1e-8
Robbybp Apr 10, 2024
c34099b
consistent format for displaying parallel tolerance
Robbybp Apr 10, 2024
d088358
Merge branch 'main' into diagnostics-parallel
Robbybp Apr 25, 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
44 changes: 40 additions & 4 deletions idaes/core/util/model_diagnostics.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -286,7 +286,7 @@ def svd_sparse(jacobian, number_singular_values):
CONFIG.declare(
"parallel_component_tolerance",
ConfigValue(
default=1e-4,
default=1e-8,
domain=float,
description="Tolerance for identifying near-parallel Jacobian rows/columns",
),
Expand Down Expand Up @@ -1179,6 +1179,29 @@ def _collect_numerical_warnings(self, jac=None, nlp=None):
self.display_constraints_with_extreme_jacobians.__name__ + "()"
)

# Parallel variables and constraints
partol = self.config.parallel_component_tolerance
par_cons = check_parallel_jacobian(
self._model, tolerance=partol, direction="row", jac=jac, nlp=nlp
)
par_vars = check_parallel_jacobian(
self._model, tolerance=partol, direction="column", jac=jac, nlp=nlp
)
if par_cons:
p = "pair" if len(par_cons) == 1 else "pairs"
warnings.append(
f"WARNING: {len(par_cons)} {p} of constraints are parallel"
f" (to tolerance {partol:.1E})"
)
next_steps.append(self.display_near_parallel_constraints.__name__ + "()")
if par_vars:
p = "pair" if len(par_vars) == 1 else "pairs"
warnings.append(
f"WARNING: {len(par_vars)} {p} of variables are parallel"
f" (to tolerance {partol:.1E})"
)
next_steps.append(self.display_near_parallel_variables.__name__ + "()")

return warnings, next_steps

def _collect_numerical_cautions(self, jac=None, nlp=None):
Expand Down Expand Up @@ -1441,7 +1464,6 @@ def report_numerical_issues(self, stream=None):
lines_list=next_steps,
title="Suggested next steps:",
line_if_empty=f"If you still have issues converging your model consider:\n"
f"{TAB*2}display_near_parallel_constraints()\n{TAB*2}display_near_parallel_variables()"
f"\n{TAB*2}prepare_degeneracy_hunter()\n{TAB*2}prepare_svd_toolbox()",
footer="=",
)
Expand Down Expand Up @@ -3619,14 +3641,25 @@ def ipopt_solve_halt_on_error(model, options=None):
)


def check_parallel_jacobian(model, tolerance: float = 1e-4, direction: str = "row"):
def check_parallel_jacobian(
model,
tolerance: float = 1e-4,
direction: str = "row",
jac=None,
nlp=None,
):
"""
Check for near-parallel rows or columns in the Jacobian.

Near-parallel rows or columns indicate a potential degeneracy in the model,
as this means that the associated constraints or variables are (near)
duplicates of each other.

For efficiency, the ``jac`` and ``nlp`` arguments may be provided if they are
already available. If these are provided, the provided model is not used. If
either ``jac`` or ``nlp`` is not provided, a Jacobian and ``PyomoNLP`` are
computed using the model.

This method is based on work published in:

Klotz, E., Identification, Assessment, and Correction of Ill-Conditioning and
Expand All @@ -3637,6 +3670,8 @@ def check_parallel_jacobian(model, tolerance: float = 1e-4, direction: str = "ro
model: model to be analysed
tolerance: tolerance to use to determine if constraints/variables are parallel
direction: 'row' (default, constraints) or 'column' (variables)
jac: model Jacobian as a ``scipy.sparse.coo_matrix``, optional
nlp: ``PyomoNLP`` of model, optional

Returns:
list of 2-tuples containing parallel Pyomo components
Expand All @@ -3651,7 +3686,8 @@ def check_parallel_jacobian(model, tolerance: float = 1e-4, direction: str = "ro
"Must be 'row' or 'column'."
)

jac, nlp = get_jacobian(model, scaled=False)
if jac is None or nlp is None:
jac, nlp = get_jacobian(model, scaled=False)

# Get vectors that we will check, and the Pyomo components
# they correspond to.
Expand Down
29 changes: 17 additions & 12 deletions idaes/core/util/tests/test_model_diagnostics.py
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -975,10 +975,10 @@ def test_display_near_parallel_variables(self):
model.v3 = Var()
model.v4 = Var()

model.c1 = Constraint(expr=model.v1 == model.v2 - 0.99999 * model.v4)
model.c2 = Constraint(expr=model.v1 + 1.00001 * model.v4 == 1e-8 * model.v3)
model.c1 = Constraint(expr=1e-8 * model.v1 == 1e-8 * model.v2 - 1e-8 * model.v4)
model.c2 = Constraint(expr=1e-8 * model.v1 + 1e-8 * model.v4 == model.v3)
model.c3 = Constraint(
expr=1e8 * (model.v1 + model.v4) + 1e10 * model.v2 == 1e-6 * model.v3
expr=1e3 * (model.v1 + model.v4) + 1e3 * model.v2 == model.v3
Comment on lines -978 to +981
Copy link
Member Author

Choose a reason for hiding this comment

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

This test was pretty fragile. The coefficients on v3 are small relative to those on other variables, so its dot product is vacuously parallel with anything. This didn't show up previously because we were relying on the absolute tolerance of 1e-4 to give it an empty sparsity pattern. Maybe this raises the question of what to do when one of the vectors we are testing is zero. I've removed the coefficients on v3 so that this isn't an issue for now.

Also, the 1e-8 ratio between small and large coefficients for v1 and v2 lead to relative "distances-from-parallel" that are very close to our new tolerance of 1e-8, so I've adjusted some coefficients to give us a little more buffer.

Copy link
Contributor

Choose a reason for hiding this comment

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

This is really because the relative tolerance should be comparing to is diff <= tolerance * unorm *vnorm, not diff <= tolerance * max(unorm, vnorm). The absolute tolerance probably can get phased out entirely.

Copy link
Member Author

Choose a reason for hiding this comment

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

Maybe unorm * vnorm is the right denominator for the relative tolerance. I would still like to have an absolute tolerance, though, as I wouldn't want (1e-9, 2e-9) to be considered parallel to (1, 2). (Unless we want to consider a zero vector as "parallel to everything", which I think is debatable.)

Copy link
Contributor

Choose a reason for hiding this comment

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

I support filtering out rows/columns with norms that are effectively zero, because otherwise they'll just spam the log as being parallel to every other row. No need to use an absolute tolerance for the inner product, though, we can just choose to skip displaying the row/column if one of the norms is too small.

In the diagnostics toolbox, the tolerance for this filtering should be the CONFIG.jacobian_small_value_warning, the same tolerance as used in display_variables_with_extreme_jacobians and display_constraints_with_extreme_jacobians to display a warning message. That way it's displayed (once) to the user there.

)

dt = DiagnosticsToolbox(model=model)
Expand Down Expand Up @@ -1121,7 +1121,7 @@ def test_collect_numerical_warnings_jacobian(self):

warnings, next_steps = dt._collect_numerical_warnings()

assert len(warnings) == 3
assert len(warnings) == 4
assert (
"WARNING: 2 Variables with extreme Jacobian values (<1.0E-08 or >1.0E+08)"
in warnings
Expand All @@ -1132,7 +1132,7 @@ def test_collect_numerical_warnings_jacobian(self):
)
assert "WARNING: 1 Constraint with large residuals (>1.0E-05)" in warnings

assert len(next_steps) == 3
assert len(next_steps) == 4
assert "display_variables_with_extreme_jacobians()" in next_steps
assert "display_constraints_with_extreme_jacobians()" in next_steps
assert "display_constraints_with_large_residuals()" in next_steps
Expand Down Expand Up @@ -1338,8 +1338,7 @@ def test_report_numerical_issues_ok(self):
Suggested next steps:

If you still have issues converging your model consider:
display_near_parallel_constraints()
display_near_parallel_variables()

prepare_degeneracy_hunter()
prepare_svd_toolbox()

Expand Down Expand Up @@ -1369,9 +1368,11 @@ def test_report_numerical_issues_exactly_singular(self):
Jacobian Condition Number: Undefined (Exactly Singular)

------------------------------------------------------------------------------------
1 WARNINGS
3 WARNINGS

WARNING: 2 Constraints with large residuals (>1.0E-05)
WARNING: 1 pair of constraints are parallel (to tolerance 1.0E-08)
WARNING: 1 pair of variables are parallel (to tolerance 1.0E-08)

------------------------------------------------------------------------------------
0 Cautions
Expand All @@ -1382,6 +1383,8 @@ def test_report_numerical_issues_exactly_singular(self):
Suggested next steps:

display_constraints_with_large_residuals()
display_near_parallel_constraints()
display_near_parallel_variables()

====================================================================================
"""
Expand Down Expand Up @@ -1433,8 +1436,8 @@ def test_report_numerical_issues_jacobian(self):
model.v2 = Var(initialize=0)
model.v3 = Var(initialize=0)

model.c1 = Constraint(expr=model.v1 == model.v2)
model.c2 = Constraint(expr=model.v1 == 1e-8 * model.v3)
model.c1 = Constraint(expr=1e-2 * model.v1 == model.v2)
model.c2 = Constraint(expr=1e-2 * model.v1 == 1e-8 * model.v3)
model.c3 = Constraint(expr=1e8 * model.v1 + 1e10 * model.v2 == 1e-6 * model.v3)
Comment on lines -1436 to 1441
Copy link
Member Author

Choose a reason for hiding this comment

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

Similarly, I updated coefficients on v1 to give some buffer around the relative tolerance. The small coefficients on v3 are also causing it to appear as "parallel to" to the other variables, which again raises the question of how to handle zero vectors.


dt = DiagnosticsToolbox(model=model)
Expand All @@ -1445,14 +1448,15 @@ def test_report_numerical_issues_jacobian(self):
expected = """====================================================================================
Model Statistics

Jacobian Condition Number: 1.407E+18
Jacobian Condition Number: 1.118E+18

------------------------------------------------------------------------------------
3 WARNINGS
4 WARNINGS

WARNING: 1 Constraint with large residuals (>1.0E-05)
WARNING: 2 Variables with extreme Jacobian values (<1.0E-08 or >1.0E+08)
WARNING: 1 Constraint with extreme Jacobian values (<1.0E-08 or >1.0E+08)
WARNING: 3 pairs of variables are parallel (to tolerance 1.0E-08)

------------------------------------------------------------------------------------
4 Cautions
Expand All @@ -1468,6 +1472,7 @@ def test_report_numerical_issues_jacobian(self):
display_constraints_with_large_residuals()
display_variables_with_extreme_jacobians()
display_constraints_with_extreme_jacobians()
display_near_parallel_variables()

====================================================================================
"""
Expand Down
Loading