Add wrapper for pyomo.contrib.iis.compute_infeasibility_explanation to diagnositics toolbox

[bknueven]

The function compute_infeasibility_explanationcomputes sets of mutually infeasible constraints, which can be much more useful for diagnosing infeasibility than the current approach of checking constraint residuals after Ipopt converges to a point of local infeasibility. It was originally developed for WaterTAP (watertap-org/watertap#1289), but a Pyomo developer took the code and put it into pyomo.contrib.iis (Pyomo/pyomo#3172).

It seems like it would be a useful tool to add to the diagnostics toolbox so IDAES users more broadly are aware of is availability. @hunterbarber has successfully leveraged it to identify infeasibility explanations for WaterTAP models.

Below is a basic demo.

import logging

from pyomo.environ import *
from pyomo.contrib.iis import compute_infeasibility_explanation

from idaes.core.solvers import get_solver

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

# for pretty-printed output; should probably fix in Pyomo
l = logging.getLogger("mis_logger")
l.setLevel(logging.INFO)
h = logging.StreamHandler()
h.setLevel(logging.INFO)
h.setFormatter(logging.Formatter(fmt=""))
l.addHandler(h)

# run the tool, using the default tolerance of `1e-6` in IDAES:
compute_infeasibility_explanation(m, solver=get_solver(), logger=l, tolerance=1e-6)

output:

WARNING: Loading a SolverResults object with a warning status into
model.name="test_infeas";
    - termination condition: infeasible
    - message from solver: Ipopt 3.13.2\x3a Converged to a locally infeasible
      point. Problem may be infeasible.
...
WARNING: Loading a SolverResults object with a warning status into
model.name="test_infeas";
    - termination condition: infeasible
    - message from solver: Ipopt 3.13.2\x3a Converged to a locally infeasible
      point. Problem may be infeasible.
Model test_infeas may be infeasible. A feasible solution was found with only the following variable bounds relaxed:
	ub of var x[1] by 0.0006513573895143176
	lb of var x[2] by 0.9993488499882835
Another feasible solution was found with only the following variable bounds relaxed:
	lb of var x[1] by 0.7071068260436827
	ub of var x[2] by 0.4142137153683113
	ub of var y by 0.707106906946579
Another feasible solution was found with only the following inequality constraints, equality constraints, and/or variable bounds relaxed:
	constraint: c by 1.0000001198315143
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:

[hunterbarber]

Recently when simulating a model, I was ramping a model variable in an attempt to approach the known physical limits (maximum possible performance) of the system. The model became infeasible at a region before the anticipated physical limits in practice.

Knowing this infeasibility was likely caused by something in the model structure, I initialized the model and first optimally solved it at the boundary of the feasible region, then took a small step into the infeasible region and applied the diagnostic tools available.

Based on the approach the model took to attempt to optimally solve the infeasible problem, the tools available (such as display_variables_near_bounds, display_constraints_with_large_residuals, and display_extreme_jacobian_entries) were not pointing to the underlying problem, and instead pointing to very distantly related equations that I could not identify the root cause. Here, compute_infeasibility_explanation clearly pointed to the problem (in this case an indexed variable was too strictly bounded and a subset of the indices were hitting a bound) that current available methods in the toolbox were not indicating.

[andrewlee94]

I think this would be a valuable tool to include; the main issues I see right now are determining how best to integrate it and finding someone with the time to do so.

Given my understanding of the tool, this is probably something that you would use if and when you get an infeasible solution returned by the solver (i.e. this is not a tool you use when you get an optimal solution but suspect you could get better performance). As such, it probably does not fit directly into the current report_structural_issues or report_numerical_issues method, but stands more as a tool to use if you encounter a specific solver failure (in this case infeasible). The question then would be more about how to inform the user that maybe they should try this method, and adding thin wrapper around this tool so it can be called directly from the toolbox.

The one other thing I see is that compute_infeasibility_explanation uses a logger for all output, which is a bit different from the rest of the diagnostics tools which use a user defined stream. This is probably a minor issue however, although it would mean a little more work.

[bknueven]

I think this would be a valuable tool to include; the main issues I see right now are determining how best to integrate it and finding someone with the time to do so.

Given my understanding of the tool, this is probably something that you would use if and when you get an infeasible solution returned by the solver (i.e. this is not a tool you use when you get an optimal solution but suspect you could get better performance). As such, it probably does not fit directly into the current report_structural_issues or report_numerical_issues method, but stands more as a tool to use if you encounter a specific solver failure (in this case infeasible). The question then would be more about how to inform the user that maybe they should try this method, and adding thin wrapper around this tool so it can be called directly from the toolbox.

At first thought, advanced IDAES users already know to check display_variables_near_bounds and display_constraints_with_large_residuals if they suspect their model might be infeasible. Perhaps we could suggest the compute_infeasibility_explanation as a follow up if neither of these methods readily identify the cause of the potential infeasibility? Admittedly it would be nice to just parse the Ipopt output and look for "converged to a point of local infeasiblity".

The one other thing I see is that compute_infeasibility_explanation uses a logger for all output, which is a bit different from the rest of the diagnostics tools which use a user defined stream. This is probably a minor issue however, although it would mean a little more work.

Maybe we should just change this in Pyomo? Given this was just merged into Pyomo/main, I don't foresee an issue with slightly modifying the API.

[andrewlee94]

1. We've tried to keep the diagnostics separate from the solver, in part because each solver is different and we want to keep the tools separate. One option however is to add compute_infeasibility_explanation as a suggested step if we see constraints with large residuals (although this can blur with poor scaling, but that just means the user needs to fix scaling first).
2. If you think changing the code in Pyomo is an option, it would definitely make life easier on the IDAES end. This is probably less critical however, as it would mostly just affect formatting of the output. I.e. we would probably just call compute_infeasibility_explanation either way and display the output to the user in whatever form it comes in.

[Robbybp]

I support adding compute_infeasibility_explanation as a suggestion if we see large residuals.

[bknueven]

On the stream front, it seems like we could wrap the pyomo.contrib version as follows (https://stackoverflow.com/questions/1383254/logging-streamhandler-and-standard-streams):

import logging
from pyomo.contrib.iis import mis

def compute_infeasibility_explanation(model, solver, tee=False, tolerance=1e-6, stream=None):
    h = logging.StreamHandler(stream)
    h.setLevel(logging.INFO)
    h.setFormatter(logging.Formatter(fmt=""))

    l = logging.getLogger("compute_infeasibility_explanation")
    l.setLevel(logging.INFO)
    l.addHandler(h)
    mis.compute_infeasibility_explanation(model, solver, tee=tee, tolerance=tolerance, logger=l)

One option however is to add compute_infeasibility_explanation as a suggested step if we see constraints with large residuals (although this can blur with poor scaling, but that just means the user needs to fix scaling first).

I support adding compute_infeasibility_explanation as a suggestion if we see large residuals.

If we're converged on this I can add the above method to the toolbox, and place the suggestion in the appropriate spot.

[Robbybp]

Why do you need to set a formatter for the StreamHandler?

[bknueven]

Why do you need to set a formatter for the StreamHandler?

Good question -- it isn't needed.
Here is the revised implementation, which creates a temporary logger instead.

import logging
from pyomo.contrib.iis import mis

def compute_infeasibility_explanation(model, solver, tee=False, tolerance=1e-6, stream=None):
    h = logging.StreamHandler(stream)
    h.setLevel(logging.INFO)

    l = logging.Logger(name="compute_infeasibility_explanation")
    l.setLevel(logging.INFO)
    l.addHandler(h)
    
    mis.compute_infeasibility_explanation(model, solver, tee=tee, tolerance=tolerance, logger=l)

[andrewlee94]

I think this plan sounds good - we can finish up any details in the PR.

[Robbybp]

Sounds good to me too.