Degrees of freedom issue with linear presolve in new solver writers

[andrewlee94]

Summary

I was running a model with the new solver writers ("ipopt_v2") and found a case where the model would solve without the linear presolve but raised an exception about too few degrees of freedom if the presolve was turned on.

Steps to reproduce the issue

from pyomo.environ import ConcreteModel, Constraint, Expression, log10, Param, Set, SolverFactory, value, Var

from idaes.core.util.math import smooth_min


def build_model():
    m = ConcreteModel()

    # Smoothing parameter for complementarities
    m.eps = Param(default=1e-3)

    m.liquid_products = Set(initialize=["K", "HCO3", "CO3", "Cl", "H2O", "OH-"])
    m.solid_products = Set(initialize=["KHCO3", "K2CO3.1.5H2O", "KCl", "K2CO3.2KHCO3.1.5H2O"])
    m.feed_species = Set(initialize=["H2O", "KHCO3", "K2CO3", "Cl"])

    m.feed_flows = Var(m.feed_species, initialize=0)
    m.concentration = Var(m.liquid_products, initialize=0, bounds=(0, None))
    m.solid_flows = Var(m.solid_products, initialize=0, bounds=(0, None))

    m.pH = Var(initialize=10, bounds=(0, 15))

    # Material balances
    m.k_balance = Constraint(
        expr=m.feed_flows["KHCO3"]
        + 2*m.feed_flows["K2CO3"]
        == m.solid_flows["KHCO3"]
        + 2*m.solid_flows["K2CO3.1.5H2O"]
        + m.solid_flows["KCl"]
        + 4*m.solid_flows["K2CO3.2KHCO3.1.5H2O"]
        + m.concentration["K"]
    )
    m.cl_balance = Constraint(
        expr=m.feed_flows["Cl"]
        == m.solid_flows["KCl"]
        + m.concentration["Cl"]
    )
    m.c_balance = Constraint(
        expr=m.feed_flows["KHCO3"]
        + m.feed_flows["K2CO3"]
        == m.solid_flows["KHCO3"]
        + m.solid_flows["K2CO3.1.5H2O"]
        + 3*m.solid_flows["K2CO3.2KHCO3.1.5H2O"]
        + m.concentration["HCO3"]
        + m.concentration["CO3"]
    )
    m.o_balance = Constraint(
        expr=3*m.feed_flows["KHCO3"]
        + 3*m.feed_flows["K2CO3"]
        + m.feed_flows["H2O"]
        == 3*m.solid_flows["KHCO3"]
        + 4.5*m.solid_flows["K2CO3.1.5H2O"]
        + 10.5*m.solid_flows["K2CO3.2KHCO3.1.5H2O"]
        + 3*m.concentration["HCO3"]
        + 3*m.concentration["CO3"]
        + m.concentration["H2O"]
        + m.concentration["OH-"]
    )

    # pH Constraint
    m.pH_constraint = Constraint(expr=14+log10(m.concentration["OH-"]) == m.pH)

    # Carbonate equilibrium
    m.carbonate_equil = Constraint(expr=1e5*m.concentration["HCO3"]*m.concentration["OH-"]==1e5*3.07e-5*m.concentration["CO3"])

    # Reaction Quotients
    def rxn_rule(m, j):
        if j == "KHCO3":
            return 13.29680294 - m.concentration["K"]*m.concentration["HCO3"]
        if j == "K2CO3.2KHCO3.1.5H2O":
            return 1e-5*234181.2733 - 1e-5*m.concentration["K"]**4*m.concentration["CO3"]*m.concentration["HCO3"]**2
        if j == "K2CO3.1.5H2O":
            return 1e-3*1432.481319 - 1e-3*m.concentration["K"]**2*m.concentration["CO3"]
        if j == "KCl":
            return 26.46242932 - m.concentration["K"]*m.concentration["Cl"]

    m.Q = Expression(m.solid_products, rule=rxn_rule)

    # Complmentarity constraints
    def comp_rule(m, j):
        if j == "KHCO3":
            return 0 == smooth_min(m.Q["KHCO3"], m.solid_flows["KHCO3"], m.eps)
        if j == "K2CO3.2KHCO3.1.5H2O":
            return 0 == smooth_min(m.Q["K2CO3.2KHCO3.1.5H2O"], m.solid_flows["K2CO3.2KHCO3.1.5H2O"], m.eps)
        if j == "K2CO3.1.5H2O":
            return 0 == smooth_min(m.Q["K2CO3.1.5H2O"], m.solid_flows["K2CO3.1.5H2O"], m.eps)
        if j == "KCl":
            return 0 == smooth_min(m.Q["KCl"], m.solid_flows["KCl"], m.eps)
        return Constraint.Skip

    m.complementarity = Constraint(m.solid_products, rule=comp_rule)

    return m


if __name__ == "__main__":
    m = build_model()

    # Set feed conditions
    m.feed_flows["H2O"].fix(55.5081)
    m.feed_flows["KHCO3"].fix(5)
    m.feed_flows["K2CO3"].fix(7.5)
    m.feed_flows["Cl"].fix(0)
    m.pH.fix(10.3576)

    # Solve model
    # If "linear_presolve": False, this will run (if with poor numerics)
    solver = SolverFactory("ipopt_v2", writer_config={"linear_presolve": True})
    solver.solve(m, tee=True)

Error Message

Ipopt 3.13.2:

******************************************************************************
This program contains Ipopt, a library for large-scale nonlinear optimization.
 Ipopt is released as open source code under the Eclipse Public License (EPL).
         For more information visit http://projects.coin-or.org/Ipopt

This version of Ipopt was compiled from source code available at
    https://github.com/IDAES/Ipopt as part of the Institute for the Design of
    Advanced Energy Systems Process Systems Engineering Framework (IDAES PSE
    Framework) Copyright (c) 2018-2019. See https://github.com/IDAES/idaes-pse.

This version of Ipopt was compiled using HSL, a collection of Fortran codes
    for large-scale scientific computation.  All technical papers, sales and
    publicity material resulting from use of the HSL codes within IPOPT must
    contain the following acknowledgement:
        HSL, a collection of Fortran codes for large-scale scientific
        computation. See http://www.hsl.rl.ac.uk.
******************************************************************************

This is Ipopt version 3.13.2, running with linear solver ma27.

Number of nonzeros in equality constraint Jacobian...:       31
Number of nonzeros in inequality constraint Jacobian.:        0
Number of nonzeros in Lagrangian Hessian.............:       18

Error in an AMPL evaluation. Run with "halt_on_ampl_error yes" to see details.
Error evaluating Jacobian of equality constraints at user provided starting point.
  No scaling factors for equality constraints computed!
Exception of type: TOO_FEW_DOF in file "IpIpoptApplication.cpp" at line 926:
 Exception message: status != TOO_FEW_DEGREES_OF_FREEDOM evaluated false: Too few degrees of freedom (rethrown)!

EXIT: Problem has too few degrees of freedom.
Traceback (most recent call last):
  File "C:\next_gen_platform\precipitation_test.py", line 107, in <module>
    solver.solve(m, tee=True)
    ~~~~~~~~~~~~^^^^^^^^^^^^^
  File "C:\Users\andrew.lee\AppData\Local\miniforge3\envs\workflow_manager\Lib\site-packages\pyomo\contrib\solver\base.py", line 598, in solve
    results: Results = super().solve(model)
                       ~~~~~~~~~~~~~^^^^^^^
  File "C:\Users\andrew.lee\AppData\Local\miniforge3\envs\workflow_manager\Lib\site-packages\pyomo\contrib\solver\ipopt.py", line 447, in solve
    raise RuntimeError(
    ...<2 lines>...
    )
RuntimeError: A feasible solution was not found, so no solution can be loaded.Please set opt.config.load_solutions=False to bypass this error.

Information on your system

Pyomo version: Pyomo 6.8.2
Python version: 3.13.1
Operating system: Windows 11
How Pyomo was installed: clone from source
Solver (if applicable): IDAES distribution of IPOPT

[jsiirola]

This is particularly nasty. If you turn on debugging (logging.getLogger('pyomo.repn').setLevel(logging.DEBUG)), you will see:

DEBUG: NL presolve: substituting concentration[Cl] := -1.0*solid_flows[KCl] + -0.0
DEBUG: NL presolve: bounds fixed solid_flows[KCl] := 0

Looking at m.cl_balance:

m.cl_balance = Constraint(
    expr=m.feed_flows["Cl"]
    == m.solid_flows["KCl"]
    + m.concentration["Cl"]
)

feed_flow["Cl"] is fixed to 0, so this constraint reduces to concentration[Cl] == -solid_flows[KCl]. The presolver does a bit of
bounds tightening here and realizes that since the concentration has a lower bound of 0, this implies an upper bound of 0 for solid_flows[KCl], which allows it to infer that the variable is actually fixed to zero.

Where we get weged is in the complimentariy constraint:

 KCl :   0.0 : 0.5*((26.46242932 - concentration[K]*concentration[Cl]) + solid_flows[KCl] - (((26.46242932 - concentration[K]*concentration[Cl]) - solid_flows[KCl])**2 + 1e-06)**0.5) :   0.0 :   True

After presolve, this can be simplified:

0 == 0.5*((26.46242932 - 0) + 0 - (((26.46242932 - 0) - 0)**2 + 1e-06)**0.5)
0 == 0.5*(26.46242932 - (26.46242932**2 + 1e-06)**0.5)
0 == -9.447356674741059e-09

The problem is the NL writer does not perform the simplification and recognize that after presolve that the constraint is constant. It is also confounded by the complementarity constraint using Expressions. So what is getting written out is

Q[KCl] := -(concentration[K] * 0) + 26.46242932
0 == 0.5 * (Q[KCl] + -(((-1 * 0) + Q[KCl]) ** 2 + 1e-6) ** 0.5 

This looks like a legitimate constraint over a legitimate defined variable ... until you realize that the defined variable is in fact constant, which makes the constraint trivial - and the constraint should be checked for feasibility and then deactivated and removes from the NL model.

All this is to say, I think I know what is going on -- but not how to fix it (without a whole lot of work re-processing the computed NL expressions, which we may have to do).