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.

score 8 · Accepted Answer · answered Jun 26 '18 at 13:52

8

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)

answered Jun 26 '18 at 13:52

Qi Chen

1,648
1
10
19

I could not have time for testing it. When I have, I accept ur answer mate. TY – oakca Jun 27 '18 at 15:32
for some reason I cannot import that package... I even installed pyomo.extras via conda. But still got no module name pyomo.util.infeasible, would you help :(? `ModuleNotFoundError: No module named 'pyomo.util.infeasible'` – oakca Aug 23 '18 at 11:57
1

https://stackoverflow.com/questions/51985436/pyomo-util-module-not-found – oakca Aug 23 '18 at 12:05
Has anyone got this working? Tried it but unable to get any kind of useful feedback along the lines of that requested. – John Curry Dec 10 '20 at 13:29
3

@JohnCurry bit of the same thing for me as well. Turns out if one constraint is infeasible, every other constraint having variables linked to that infeasible constraint will also turn out to be infeasible when using this, hence it is difficult to have just one particular constraint printed out. This is alot worse when using sets/indices as constraints are treated in blocks. I usually treat each constraint on a case by case basis deactivating them (constraint.deactivate()) serially, whilst ensuring the model still makes sense, until I narrow down the erring constraint. Hope this helps. – EJay Dec 11 '20 at 16:10
@EJay, OK thanks. I do the same thing. Would be nice to get some functionality on that as it can be a laborious trial and error process. But good to know there isnt any functionality that exists that I'm missing. – John Curry Dec 11 '20 at 17:05
The solvers cbc and gurobi both report my problem is infeasible, but `log_infeasible_constraints(model)` and `log_infeasible_bounds(model)` doesn't print anything to screen. – jmlarson Apr 21 '21 at 19:13
1

same problem of @jmlarson – user3254491 Feb 07 '22 at 14:22

score 2 · Answer 2 · answered Jan 29 '22 at 22:52

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

Finding out reason of Pyomo model infeasibility

2 Answers2

An example:

Reports:

Linked