diff --git a/.github/workflows/typos.toml b/.github/workflows/typos.toml index c0633eb4d3..b335b9a959 100644 --- a/.github/workflows/typos.toml +++ b/.github/workflows/typos.toml @@ -36,6 +36,9 @@ equil = "equil" astroid = "astroid" # delt == delta delt = "delt" +# Minimal Infeasible System +mis = "mis" +MIS = "MIS" # TEMPORARY - Remove after example update upadate = "upadate" @@ -44,4 +47,4 @@ RTO = "RTO" PN = "PN" hd = "hd" Tge = "Tge" -iy = "iy" \ No newline at end of file +iy = "iy" diff --git a/idaes/core/util/model_diagnostics.py b/idaes/core/util/model_diagnostics.py index 7fb08393c0..92ca93da3b 100644 --- a/idaes/core/util/model_diagnostics.py +++ b/idaes/core/util/model_diagnostics.py @@ -23,6 +23,7 @@ from math import log, isclose, inf, isfinite import json from typing import List +import logging import numpy as np from scipy.linalg import svd @@ -73,10 +74,12 @@ 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.contrib.iis import mis from pyomo.common.deprecation import deprecation_warning from pyomo.common.errors import PyomoException from pyomo.common.tempfiles import TempfileManager +from idaes.core.solvers.get_solver import get_solver from idaes.core.util.model_statistics import ( activated_blocks_set, deactivated_blocks_set, @@ -291,6 +294,14 @@ def svd_sparse(jacobian, number_singular_values): description="Tolerance for identifying near-parallel Jacobian rows/columns", ), ) +CONFIG.declare( + "absolute_feasibility_tolerance", + ConfigValue( + default=1e-6, + domain=float, + description="Feasibility tolerance for identifying infeasible constraints and bounds", + ), +) SVDCONFIG = ConfigDict() @@ -725,6 +736,50 @@ def display_constraints_with_large_residuals(self, stream=None): footer="=", ) + def compute_infeasibility_explanation(self, stream=None, solver=None, tee=False): + """ + This function attempts to determine why a given model is infeasible. It deploys + two main algorithms: + + 1. Relaxes the constraints of the problem, and reports to the user + some sets of constraints and variable bounds, which when relaxed, creates a + feasible model. + 2. Uses the information collected from (1) to attempt to compute a Minimal + Infeasible System (MIS), which is a set of constraints and variable bounds + which appear to be in conflict with each other. It is minimal in the sense + that removing any single constraint or variable bound would result in a + feasible subsystem. + + Args: + stream: I/O object to write report to (default = stdout) + solver: A pyomo solver object or a string for SolverFactory + (default = get_solver()) + tee: Display intermediate solves conducted (False) + + Returns: + None + + """ + if solver is None: + solver = get_solver() + if stream is None: + stream = sys.stdout + + h = logging.StreamHandler(stream) + h.setLevel(logging.INFO) + + l = logging.Logger(name=__name__ + ".compute_infeasibility_explanation") + l.setLevel(logging.INFO) + l.addHandler(h) + + mis.compute_infeasibility_explanation( + self._model, + solver, + tee=tee, + tolerance=self.config.absolute_feasibility_tolerance, + logger=l, + ) + def get_dulmage_mendelsohn_partition(self): """ Performs a Dulmage-Mendelsohn partitioning on the model and returns @@ -1123,6 +1178,7 @@ def _collect_numerical_warnings(self, jac=None, nlp=None): next_steps.append( self.display_constraints_with_large_residuals.__name__ + "()" ) + next_steps.append(self.compute_infeasibility_explanation.__name__ + "()") # Variables outside bounds violated_bounds = _vars_violating_bounds( diff --git a/idaes/core/util/tests/test_model_diagnostics.py b/idaes/core/util/tests/test_model_diagnostics.py index 3b1312521a..d502c4ebf8 100644 --- a/idaes/core/util/tests/test_model_diagnostics.py +++ b/idaes/core/util/tests/test_model_diagnostics.py @@ -1088,8 +1088,9 @@ def test_collect_numerical_warnings(self, model): assert "WARNING: 1 Constraint with large residuals (>1.0E-05)" in warnings assert "WARNING: 1 Variable at or outside bounds (tol=0.0E+00)" in warnings - assert len(next_steps) == 2 + assert len(next_steps) == 3 assert "display_constraints_with_large_residuals()" in next_steps + assert "compute_infeasibility_explanation()" in next_steps assert "display_variables_at_or_outside_bounds()" in next_steps @pytest.mark.component @@ -1137,10 +1138,11 @@ def test_collect_numerical_warnings_jacobian(self): ) assert "WARNING: 1 Constraint with large residuals (>1.0E-05)" in warnings - assert len(next_steps) == 4 + assert len(next_steps) == 5 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 + assert "compute_infeasibility_explanation()" in next_steps @pytest.mark.component def test_collect_numerical_cautions(self, model): @@ -1392,6 +1394,7 @@ def test_report_numerical_issues_exactly_singular(self): Suggested next steps: display_constraints_with_large_residuals() + compute_infeasibility_explanation() display_near_parallel_constraints() display_near_parallel_variables() @@ -1431,6 +1434,7 @@ def test_report_numerical_issues(self, model): Suggested next steps: display_constraints_with_large_residuals() + compute_infeasibility_explanation() display_variables_at_or_outside_bounds() ==================================================================================== @@ -1479,6 +1483,7 @@ def test_report_numerical_issues_jacobian(self): Suggested next steps: display_constraints_with_large_residuals() + compute_infeasibility_explanation() display_variables_with_extreme_jacobians() display_constraints_with_extreme_jacobians() display_near_parallel_variables() @@ -4010,3 +4015,38 @@ def test_columns(self, afiro): (afiro.X15, pytest.approx(-0.042451666, rel=1e-5)), (afiro.X37, pytest.approx(-0.036752232, rel=1e-5)), ] + + +class TestComputeInfeasibilityExplanation: + + @pytest.fixture(scope="class") + def model(self): + # create an infeasible model for demonstration + m = ConcreteModel() + + m.name = "test_infeas" + m.x = Var([1, 2], bounds=(0, 1)) + m.y = Var(bounds=(0, 1)) + + m.c = Constraint(expr=m.x[1] * m.x[2] == -1) + m.d = Constraint(expr=m.x[1] + m.y >= 1) + + return m + + @pytest.mark.component + @pytest.mark.solver + def test_output(self, model): + dt = DiagnosticsToolbox(model) + + stream = StringIO() + + dt.compute_infeasibility_explanation(stream=stream) + + expected = """Computed Minimal Intractable System (MIS)! +Constraints / bounds in MIS: + lb of var x[2] + lb of var x[1] + constraint: c +Constraints / bounds in guards for stability: +""" + assert expected in stream.getvalue()