Infeasibility diagnostic tool (#1409)

* add wrapper for pyomo.contrib.iis.mis.compute_infeasibility_explanation

* testing basic stream functionality

* integrate into model diagnostics toolbox

* spelling

* remove separate function

.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"
+iy = "iy"

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(

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