Finding out reason of Pyomo model infeasibility

Question

I got a pyomo concrete model with lots of variables and constraints.

Somehow, one of the variable inside my model violates one constraint, which makes my model infeasible:

WARNING: Loading a SolverResults object with a warning status into model=xxxx;
    message from solver=Model was proven to be infeasible.

Is there a way to ask the solver, the reason of the infeasibility?

So for example, lets assume I got a variable called x, and if I define following 2 constraints, model will be ofc. infeasible.

const1:
    x >= 10

const2:
    x <= 5

And what I want to achieve that pointing out the constraints and variable which causes this infeasibility, so that I can fix it. Otherwise with my big model it is kinda hard to get what causing this infeasibility.

IN: write_some_comment
OUT: variable "x" cannot fulfill "const1" and "const2" at the same time.

Answer 1

Many solvers (including IPOPT) will hand you back the value of the variables at solver termination, even if the problem was found infeasible. At that point, you do have some options.

There is contributed code in pyomo.util.infeasible that might help you out. https://github.com/Pyomo/pyomo/blob/master/pyomo/util/infeasible.py

Usage:

from pyomo.util.infeasible import log_infeasible_constraints
...
SolverFactory('your_solver').solve(model)
...
log_infeasible_constraints(model)

Answer 2

I would not trust any numbers that the solver loads into the model after reporting "infeasible." I don't think any solvers come w/ guarantees on the validity of those numbers. Further, unless a package can divine the modeler's intent, it isn't clear how it would list the infeasible constraints. Consider 2 constraints:

C1:  x <= 5
C2:  x >= 10

X ∈ Reals, or Integers, ...

Which is the invalid constraint? Well, it depends! Point being, it seems an impossible task to unwind the mystery based on values the solver tries.

A possible alternate strategy: Load the model with what you believe to be a valid solution, and test the slack on the constraints. This "loaded solution" could even be a null case where everything is zero'ed out (if that makes sense in the context of the model). It could also be a set of known feasible solutions tried via unit test code.

If you can construct what you believe to be a valid solution (forget about optimal, just something valid), you can (1) load those values, (2) iterate through the constraints in the model, (3) evaluate the constraint and look for negative slack, and (4) report the culprits with values and expressions

An example:

import pyomo.environ as pe

test_null_case = True

m = pe.ConcreteModel('sour constraints')

# SETS
m.T = pe.Set(initialize=['foo', 'bar'])

# VARS
m.X = pe.Var(m.T)
m.Y = pe.Var()

# OBJ
m.obj = pe.Objective(expr = sum(m.X[t] for t in m.T) + m.Y)

# Constraints
m.C1 = pe.Constraint(expr=sum(m.X[t] for t in m.T) <= 5)
m.C2 = pe.Constraint(expr=sum(m.X[t] for t in m.T) >= 10)
m.C3 = pe.Constraint(expr=m.Y >= 7)
m.C4 = pe.Constraint(expr=m.Y <= sum(m.X[t] for t in m.T))

if test_null_case:
    # set values of all variables to a "known good" solution...
    m.X.set_values({'foo':1, 'bar':3})  # index:value
    m.Y.set_value(2)  # scalar
    for c in m.component_objects(ctype=pe.Constraint):
        if c.slack() < 0:  # constraint is not met
            print(f'Constraint {c.name} is not satisfied')
            c.display()  # show the evaluation of c
            c.pprint()   # show the construction of c
            print()
else:
    pass
    # instantiate solver & solve, etc...

Reports:

Constraint C2 is not satisfied
C2 : Size=1
    Key  : Lower : Body : Upper
    None :  10.0 :    4 :  None
C2 : Size=1, Index=None, Active=True
    Key  : Lower : Body            : Upper : Active
    None :  10.0 : X[foo] + X[bar] :  +Inf :   True

Constraint C3 is not satisfied
C3 : Size=1
    Key  : Lower : Body : Upper
    None :   7.0 :    2 :  None
C3 : Size=1, Index=None, Active=True
    Key  : Lower : Body : Upper : Active
    None :   7.0 :    Y :  +Inf :   True

Answer 3

Solution of Qi Chen wasn't working for me and it wasn't printing anything to my screen. I debugged and I saw that you have to use a logger as well with at least INFO level.

import logging
from pyomo.util.infeasible import log_infeasible_constraints

def log_pyomo_infeasible_constraints(model_instance):          
    # Create a logger object with DEBUG level
    logging_logger = logging.getLogger()
    logging_logger.setLevel(logging.DEBUG)
    # Create a console handler
    ch = logging.StreamHandler()
    ch.setLevel(logging.DEBUG)
    # add the handler to the logger
    logging_logger.addHandler(ch)
    # Log the infeasible constraints of pyomo object
    print("Displaying Infeasible Constraints")
    log_infeasible_constraints(model_instance, log_expression=True,
                         log_variables=True, logger=logging_logger)