Degree of freedom not equal zero for PDE problems using pyomo.dae

[lxhowl]

Summary

I use pyomo.dae to solve PDE problems but found the DOF is not equal to 0 and also changes using different discretization transformations.

Steps to reproduce the issue

# example.py
import pyomo.environ as pyo
import pyomo.dae as dae

# parameters
tf = 1.00
Ys = 10.0
alpha = 5.0
beta = 5.0

m = pyo.ConcreteModel()

m.t = dae.ContinuousSet(bounds=(0, tf))
m.x = dae.ContinuousSet(bounds=(0, 1))

m.y = pyo.Var(m.t, m.x)
m.s = pyo.Var(m.t, m.x)

m.dydt = dae.DerivativeVar(m.y, wrt=m.t)
m.dydx = dae.DerivativeVar(m.y, wrt=m.x)
m.d2ydx2 = dae.DerivativeVar(m.y, wrt=(m.x, m.x))

@m.Constraint(m.t, m.x)
def pde(m, t, x):
    return beta * m.dydt[t, x]  == alpha *  m.d2ydx2[t, x]

@m.Constraint(m.t)
def bc1(m, t):
    return m.y[t, 0] == Ys

@m.Constraint(m.t)
def bc2(m, t):
    return m.dydx[t, 1] == 0

@m.Constraint(m.x)
def ic(m, x):
    if x == 0:
        return pyo.Constraint.Skip
    return m.y[0, x] == 0.0

# transform and solve
pyo.TransformationFactory('dae.finite_difference').apply_to(m, wrt=m.x, nfe=20)
pyo.TransformationFactory('dae.finite_difference').apply_to(m, wrt=m.t, nfe=30)
pyo.SolverFactory('ipopt').solve(m).write()
...

Error Message

Problem:
Lower bound: -inf
Upper bound: inf
Number of objectives: 1
Number of constraints: 2572
Number of variables: 2573
Sense: unknown

Information on your system

Pyomo version: 6.6.2
Python version: 3.10.12
Operating system: Linux 5.15.120+
How Pyomo was installed (PyPI, conda, source): idaes_pse
Solver (if applicable): Ipopt 3.13.2

Additional information

Colab notebook for the minimal problem
Same problem when switching discretizers in pyomo PDE example

[Robbybp]

Hi @lxhowl, user-support type questions like this are probably better suited for the Google group or Stack Overflow, but we do have some tools available for debugging degree of freedom issues. I can add the following lines to the end of your script:

from ideas.core.util import DiagnosticsToolbox
dt = DiagnosticsToolbox(m)
dt.report_structural_issues()

This gives me:

====================================================================================
Model Statistics

        Activated Blocks: 1 (Deactivated: 0)
        Free Variables in Activated Constraints: 2573 (External: 0)
            Free Variables with only lower bounds: 0
            Free Variables with only upper bounds: 0
            Free Variables with upper and lower bounds: 0
        Fixed Variables in Activated Constraints: 0 (External: 0)
        Activated Equality Constraints: 2572 (Deactivated: 0)
        Activated Inequality Constraints: 0 (Deactivated: 0)
        Activated Objectives: 0 (Deactivated: 0)

------------------------------------------------------------------------------------
2 WARNINGS

    WARNING: 1 Degree of Freedom
    WARNING: Structural singularity found
        Under-Constrained Set: 4 variables, 2 constraints
        Over-Constrained Set: 3 variables, 4 constraints

------------------------------------------------------------------------------------
1 Cautions

    Caution: 682 unused variables (0 fixed)

------------------------------------------------------------------------------------
Suggested next steps:

    display_underconstrained_set()
    display_overconstrained_set()

====================================================================================

This indicates that you have two degree-of-freedom issues in the model you've written: A set of 2 constraints that is attempting to solve for 4 variables, and a set of 4 constraints that is attempting to solve for 3 variables, for a net of 1 degree of freedom. We can see what variables/constraints are in these under and over-constrained sets by calling the suggested methods:

dt.display_underconstrained_set()
dt.display_overconstrained_set()

which gives me:

====================================================================================
Dulmage-Mendelsohn Under-Constrained Set

    Independent Block 0:

        Variables:

            d2ydx2[0,0]
            dydt[0,0]

        Constraints:

            pde[0,0]

    Independent Block 1:

        Variables:

            d2ydx2[0,0.05]
            dydt[0,0.05]

        Constraints:

            pde[0,0.05]

====================================================================================
====================================================================================
Dulmage-Mendelsohn Over-Constrained Set

    Independent Block 0:

        Variables:

            dydx[0,1]
            y[0,0.95]
            y[0,1]

        Constraints:

            bc2[0]
            ic[0.95]
            ic[1]
            dydx_disc_eq[0,1]

====================================================================================

I did not look closely enough at the model to tell you exactly why these sets are under and over-constrained, but I would probably start by asking why the pde constraint is attempting to "solve for" two variables at the (0,0) index. For more information about what these "under and over-constrained sets" are, you can see the Wikipedia page for the Dulmage-Mendelsohn decomposition or this paper.

[adowling2]

@Robbybp Thanks! We'll look at this more closely. This could even become a good example for the Diagnostics toolbox in IDAES. How to debug an over or underconstrained PDE model.

Tagging: @andrewlee94 and @blnicho for their awareness.