Infeasibility diagnostic tool

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,6 +74,7 @@
 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
@@ -291,6 +293,14 @@
         description="Tolerance for identifying near-parallel Jacobian rows/columns",
     ),
 )
+CONFIG.declare(
+    "absolute_feasibility_tolerance",
+    ConfigValue(
+        default=1e-6,
+        domain=float,
+        description="Feasiblity tolerance for idenifying infeasible constraint and bounds",
+    ),
+)
 
 
 SVDCONFIG = ConfigDict()
@@ -725,6 +735,36 @@
             footer="=",
         )
 
+    def compute_infeasibility_explanation(self, solver, stream=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:
+            solver: A pyomo solver object or a string for SolverFactory
+            stream: I/O object to write report to (default = stdout)
+            tee:  Display intermediate solves conducted (False)
+
+        Returns:
+            None
+
+        """
+        compute_infeasibility_explanation(
+            self._model,
+            solver,
+            tolerance=self.config.absolute_feasibility_tolerance,
+            stream=stream,
+        )
+
     def get_dulmage_mendelsohn_partition(self):
         """
         Performs a Dulmage-Mendelsohn partitioning on the model and returns
@@ -1123,6 +1163,9 @@
             next_steps.append(
                 self.display_constraints_with_large_residuals.__name__ + "()"
             )
+            next_steps.append(
+                self.compute_infeasibility_explanation.__name__ + "(solver=)"
+            )
 
         # Variables outside bounds
         violated_bounds = _vars_violating_bounds(
@@ -3887,6 +3930,51 @@
     return ill_cond
 
 
+def compute_infeasibility_explanation(
+    model, solver, tee=False, tolerance=1e-6, stream=None
+):
+    """
+    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.
+
+    This function is a wrapper for the same capability in pyomo.contrib.iis.mis.
+
+    Args:
+        model: A pyomo block
+        solver: A pyomo solver object or a string for SolverFactory
+        tee:  Display intermediate solves conducted (False)
+        tolerance: The feasibility tolerance to use when declaring a
+            constraint feasible (1e-08)
+        stream: I/O object to write report to (default = stdout)
+
+    Returns:
+        None
+
+    """
+    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(
+        model, solver, tee=tee, tolerance=tolerance, logger=l
+    )
+
+
 # -------------------------------------------------------------------------------------------
 # Private module functions
 def _var_in_block(var, block):

idaes/core/util/tests/test_model_diagnostics.py

@@ -81,6 +81,7 @@
     _collect_model_statistics,
     check_parallel_jacobian,
     compute_ill_conditioning_certificate,
+    compute_infeasibility_explanation,
 )
 from idaes.core.util.parameter_sweep import (
     SequentialSweepRunner,
@@ -1088,8 +1089,9 @@
         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(solver=)" in next_steps
         assert "display_variables_at_or_outside_bounds()" in next_steps
 
     @pytest.mark.component
@@ -1137,10 +1139,11 @@
         )
         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(solver=)" in next_steps
 
     @pytest.mark.component
     def test_collect_numerical_cautions(self, model):
@@ -1392,6 +1395,7 @@
 Suggested next steps:
 
     display_constraints_with_large_residuals()
+    compute_infeasibility_explanation(solver=)
     display_near_parallel_constraints()
     display_near_parallel_variables()
 
@@ -1431,6 +1435,7 @@
 Suggested next steps:
 
     display_constraints_with_large_residuals()
+    compute_infeasibility_explanation(solver=)
     display_variables_at_or_outside_bounds()
 
 ====================================================================================
@@ -1479,6 +1484,7 @@
 Suggested next steps:
 
     display_constraints_with_large_residuals()
+    compute_infeasibility_explanation(solver=)
     display_variables_with_extreme_jacobians()
     display_constraints_with_extreme_jacobians()
     display_near_parallel_variables()
@@ -4010,3 +4016,37 @@
             (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):
+        stream = StringIO()
+
+        solver = get_solver()
+        compute_infeasibility_explanation(model, solver, 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()