Unavoidable InfeasibleConstraintException in NL writer for variable fixed outside of bounds

[Robbybp]

Summary

If a variable is fixed to a value outside of its bounds, the NL writer throws an InfeasibleConstraintException in cache_fixed_var. This may be reasonable default behavior, but I would argue there are cases where you would want to ignore this.

My particular use case is an initialization method where I solve a sequence of small subproblems, sometimes with solvers that do not support bounds (e.g. calculate_variable_from_constraint or pyo.SolverFactory("scipy.fsolve")). The solutions may violate bounds, and I may fix some of the resulting variable values before a subsequent solve. If this subsequent solve uses the NL writer and any of the fixed variables violate their bounds, I get the above error.

The workaround is to temporarily remove bounds on fixed variables before each sub-problem solve. This is not difficult, but I thought this might be a common enough issue to warrant an option in the NL writer.

Steps to reproduce the issue

A basic example of my workflow that gives me the error is:

import pyomo.environ as pyo                                                                                                                                                                                                                                                     
from pyomo.util.calc_var_value import calculate_variable_from_constraint                                                                                                                                                                                                        
                                                                                                                                                                                                                                                                                
m = pyo.ConcreteModel()                                                                                                                                                                                                                                                         
                                                                                                                                                                                                                                                                                
m.x = pyo.Var([1,2], initialize=1.0, bounds=(0, None))                                                                                                                                                                                                                          
m.u = pyo.Var([1,2], initialize=1.0, bounds=(0, None))                                                                                                                                                                                                                          
                                                                                                                                                                                                                                                                                
m.eq = pyo.Constraint(pyo.PositiveIntegers)                                                                                                                                                                                                                                     
m.eq[1] = m.u[1] == 5                                                                                                                                                                                                                                                           
m.eq[2] = m.u[1] + 2 * m.u[2] == -1 + m.x[2]                                                                                                                                                                                                                                    
m.eq[3] = m.x[1] * m.x[2] == m.u[1] + m.u[2] + 7                                                                                                                                                                                                                                
                                                                                                                                                                                                                                                                                
m.obj = pyo.Objective(expr=m.x[1]**2 + m.x[2]**2)                                                                                                                                                                                                                               
                                                                                                                                                                                                                                                                                
calculate_variable_from_constraint(m.u[1], m.eq[1])                                                                                                                                                                                                                             
m.u[1].fix()                                                                                                                                                                                                                                                                    
m.eq[1].deactivate()                                                                                                                                                                                                                                                            
calculate_variable_from_constraint(m.u[2], m.eq[2])                                                                                                                                                                                                                             
m.u[2].fix()                                                                                                                                                                                                                                                                    
m.eq[2].deactivate()                                                                                                                                                                                                                                                            
pyo.SolverFactory("ipopt").solve(m)                                                                                                                                                                                                                                             

Error Message

This gives me the error:

WARNING (W1002): Setting Var 'u[2]' to a numeric value `-2.5` outside the
bounds (0, None).
  ...
  File "/path/to/pyomo/pyomo/repn/plugins/nl_writer.py", line 2877, in cache_fixed_var
    raise InfeasibleConstraintException(
pyomo.common.errors.InfeasibleConstraintException: model contains a trivially infeasible variable 'u[2]' (fixed value -2.5 outside bounds [0, None]).

Information on your system

Pyomo version: Recent main
Python version: 3.9.6

[michaelbynum]

I think an option in the NL writer is a great idea.