Degrees of freedom drop to -5 after initialization when adding Heater unit in a simple system

[Hanli131]

Hi everyone,

I've built a simple RSOC system, as shown in Figure 1.

The system first preheats the fuel and air streams (from room temperature to 1073 K), and then feeds them into the RSOC for hydrogen production. A Heater unit is used to ensure that the inlet gas temperature reaches 1073 K.

After initializing the system (see lines 284–313 in the code), the degrees of freedom drop from 0 to -5.

I found that this issue is caused by the oxygenheater unit. When I remove the oxygenheater and directly connect the oxygenfeed to the stack_oxygen_inlet, the system looks like Figure 2, and initialization completes without issues—the degrees of freedom remain at 0 and everything runs correctly.

I'm quite confused and would really appreciate your insights:
Why does the degrees of freedom drop after initialization when the heater is included, and how can this be resolved properly?

Additional information for reference:

1. The RSOC model and configuration are based on the following example:
   https://github.com/IDAES/idaes-pse/blob/main/idaes/models_extra/power_generation/unit_models/soc_submodels/tests/test_herring_replication.py
   I’ve verified that this RSOC model runs correctly under the given flowsheet conditions.

2. Despite the negative degrees of freedom after initialization, the model can still be solved and produces the same results as the version without the oxygenheater.

3. The code is included below.
   SOC_system.zip
   To reproduce the version shown in Figure 2 (which runs successfully), please:
   Enable the commented-out lines (213, 229),
   Disable the oxygenheater-related lines (193, 216–217, 228, 232).

Thank you in advance for your help!

import numpy as np
import pyomo.environ as pyo
from idaes.core import FlowsheetBlock
from idaes.models.properties import iapws95
from idaes.core.util.model_statistics import degrees_of_freedom
from idaes.core.solvers import get_solver

from idaes.models.properties.modular_properties.base.generic_property import (
    GenericParameterBlock,
)
from idaes.models_extra.power_generation.properties.natural_gas_PR import (
    get_prop,
    EosType,
)
from idaes.models_extra.power_generation.unit_models.soc_submodels import (
    SolidOxideModuleSimple,
)
from idaes.models_extra.power_generation.unit_models.soc_submodels.common import (
    _comp_enthalpy_expr,
    _comp_entropy_expr,
)
from idaes.models.unit_models.heat_exchanger import HeatExchangerFlowPattern
from idaes.models.unit_models import Feed
from idaes.models.unit_models.heater import Heater
from idaes.models.unit_models.heat_exchanger import HeatExchanger, delta_temperature_amtd_callback
import idaes.logger as idaeslog


def cccm_to_mps(V):
    # Takes cccm, returns moles per second
    V = V * 1e-6 / 60
    P = 101300
    R = 8.31
    T = 273.15
    return V * P / (R * T)

def mps_to_cccm(n):
    P = 101300
    R = 8.31
    T = 273.15
    V = n * R * T / P
    return V * 1e6 * 60


N_cell = 10
P = 101325
_constF = 96485

zfaces = np.linspace(0, 1, 11).tolist()

xfaces_electrode = [0.0, 1.0]
xfaces_electrolyte = [0.0, 1.0]

fuel_comps = ["H2", "H2O"]
fuel_stoich_dict = {
    "H2": -0.5,
    "H2O": 0.5,
    "Vac": 0.5,
    "O^2-": -0.5,
    "e^-": 1.0,
}
oxygen_comps = ["O2", "N2"]
oxygen_stoich_dict = {"O2": -0.25, "Vac": -0.5, "O^2-": 0.5, "e^-": -1.0}

cell_config = {
    "has_holdup": True,
    "control_volume_zfaces": zfaces,
    "control_volume_xfaces_fuel_electrode": xfaces_electrode,
    "control_volume_xfaces_oxygen_electrode": xfaces_electrode,
    "control_volume_xfaces_electrolyte": xfaces_electrolyte,
    "fuel_component_list": fuel_comps,
    "fuel_triple_phase_boundary_stoich_dict": fuel_stoich_dict,
    "inert_fuel_species_triple_phase_boundary": [],
    "oxygen_component_list": oxygen_comps,
    "oxygen_triple_phase_boundary_stoich_dict": oxygen_stoich_dict,
    "inert_oxygen_species_triple_phase_boundary": ["N2"],
    "flow_pattern": HeatExchangerFlowPattern.cocurrent,
    "include_temperature_x_thermo": True,
    "include_contact_resistance": True,
}


def fix_cell_parameters(cell):
    cell.fuel_channel.length_x.fix(1.09e-3)
    cell.oxygen_channel.length_x.fix(1.09e-3)

    cell.length_y.fix(0.08)

    cell.length_z.fix(0.08)
    cell.fuel_channel.heat_transfer_coefficient.fix(100)

    cell.contact_flow_mesh_fuel_electrode.log_preexponential_factor.fix(
        pyo.log(0.55e-4)
    )
    cell.contact_flow_mesh_oxygen_electrode.log_preexponential_factor.fix(
        pyo.log(0.55e-4)
    )
    cell.contact_interconnect_fuel_flow_mesh.log_preexponential_factor.fix(
        pyo.log(0.15e-4)
    )
    cell.contact_interconnect_oxygen_flow_mesh.log_preexponential_factor.fix(
        pyo.log(0.15e-4)
    )

    cell.contact_flow_mesh_fuel_electrode.thermal_exponent_dividend.fix(0)
    cell.contact_flow_mesh_fuel_electrode.contact_fraction.fix(1)
    cell.contact_flow_mesh_oxygen_electrode.thermal_exponent_dividend.fix(0)
    cell.contact_flow_mesh_oxygen_electrode.contact_fraction.fix(1)

    cell.contact_interconnect_fuel_flow_mesh.thermal_exponent_dividend.fix(0)
    cell.contact_interconnect_fuel_flow_mesh.contact_fraction.fix(1)
    cell.contact_interconnect_oxygen_flow_mesh.thermal_exponent_dividend.fix(0)
    cell.contact_interconnect_oxygen_flow_mesh.contact_fraction.fix(1)

    cell.oxygen_channel.heat_transfer_coefficient.fix(100)

    cell.fuel_electrode.length_x.fix(23e-6)
    cell.fuel_electrode.porosity.fix(0.37)
    cell.fuel_electrode.tortuosity.fix(3)
    cell.fuel_electrode.solid_heat_capacity.fix(450)
    cell.fuel_electrode.solid_density.fix(3210.0)
    cell.fuel_electrode.solid_thermal_conductivity.fix(2.16)
    cell.fuel_electrode.resistivity_log_preexponential_factor.fix(pyo.log(8.856e-6))
    cell.fuel_electrode.resistivity_thermal_exponent_dividend.fix(0)

    cell.oxygen_electrode.length_x.fix(21e-6)
    cell.oxygen_electrode.porosity.fix(0.37)
    cell.oxygen_electrode.tortuosity.fix(3.0)
    cell.oxygen_electrode.solid_heat_capacity.fix(430)
    cell.oxygen_electrode.solid_density.fix(3030)
    cell.oxygen_electrode.solid_thermal_conductivity.fix(2.16)
    
        
    cell.oxygen_electrode.resistivity_log_preexponential_factor.fix(pyo.log(1.425e-4))
    cell.oxygen_electrode.resistivity_thermal_exponent_dividend.fix(0)

    cell.electrolyte.length_x.fix(1.4e-4)
    cell.electrolyte.heat_capacity.fix(470)
    cell.electrolyte.density.fix(5160)
    cell.electrolyte.thermal_conductivity.fix(2.16)
    cell.electrolyte.resistivity_log_preexponential_factor.fix(pyo.log(1.07e-4))
    cell.electrolyte.resistivity_thermal_exponent_dividend.fix(7237)
    
    # Values from Kazempoor and Braun
    cell.fuel_triple_phase_boundary.exchange_current_log_preexponential_factor.fix(
        pyo.log(5.5e8)
    )
    cell.fuel_triple_phase_boundary.exchange_current_activation_energy.fix(100e3)
    cell.fuel_triple_phase_boundary.activation_potential_alpha1.fix(1)
    cell.fuel_triple_phase_boundary.activation_potential_alpha2.fix(1)
    cell.fuel_triple_phase_boundary.exchange_current_exponent_comp["H2"].fix(0.5)
    cell.fuel_triple_phase_boundary.exchange_current_exponent_comp["H2O"].fix(0.1)

    cell.oxygen_triple_phase_boundary.exchange_current_log_preexponential_factor.fix(
        pyo.log(7e8)
    )
    cell.oxygen_triple_phase_boundary.exchange_current_activation_energy.fix(110e3)
    cell.oxygen_triple_phase_boundary.activation_potential_alpha1.fix(1)
    cell.oxygen_triple_phase_boundary.activation_potential_alpha2.fix(1)
    cell.oxygen_triple_phase_boundary.exchange_current_exponent_comp["O2"].fix(0.5)


if __name__ == "__main__":
    m = pyo.ConcreteModel()
    m.fs = FlowsheetBlock(dynamic=False, time_set=[0], time_units=pyo.units.s)

    m.fs.propertiesIapws95 = iapws95.Iapws95ParameterBlock()
    m.fs.prop_Iapws95 = iapws95.Iapws95StateBlock(parameters=m.fs.propertiesIapws95)
    
    m.fs.oxygen_params = GenericParameterBlock(
        **get_prop(oxygen_comps, {"Vap"}, eos=EosType.IDEAL), doc="Air-side parameters"
    )
    m.fs.fuel_params = GenericParameterBlock(
        **get_prop(fuel_comps, {"Vap"}, eos=EosType.IDEAL), doc="Fuel-side parameters"
    )
 
    for pp in [m.fs.fuel_params, m.fs.oxygen_params]:
        pp.set_default_scaling("flow_mol", 1e3)
        pp.set_default_scaling("flow_mol_phase", 1e3)
        pp.set_default_scaling("mole_frac_comp", 10)
        pp.set_default_scaling("mole_frac_phase_comp", 10)
        pp.set_default_scaling("enth_mol_phase", 1e-3)
        pp.set_default_scaling("pressure", 1e-5)
        pp.set_default_scaling("temperature", 1e-2)
        
    m.fs.fuelfeed = Feed(property_package=m.fs.fuel_params)
    m.fs.fuelheater = Heater(property_package=m.fs.fuel_params)
    m.fs.oxygenfeed = Feed(property_package=m.fs.oxygen_params)
    m.fs.oxygenheater = Heater(property_package=m.fs.oxygen_params)
    
    m.fs.fuelheat_exchanger = HeatExchanger(
            delta_temperature_callback=delta_temperature_amtd_callback,
            hot_side_name="shell",
            cold_side_name="tube",
            shell={"property_package": m.fs.fuel_params},
            tube={"property_package": m.fs.fuel_params})

    m.fs.stack = SolidOxideModuleSimple(
        solid_oxide_cell_config=cell_config,
        fuel_property_package=m.fs.fuel_params,
        oxygen_property_package=m.fs.oxygen_params,
    )
    
    from pyomo.network import Arc, SequentialDecomposition

    m.fs.s00 = Arc(source=m.fs.fuelfeed.outlet, destination=m.fs.fuelheat_exchanger.cold_side_inlet)
    m.fs.s01 = Arc(source=m.fs.fuelheat_exchanger.cold_side_outlet, destination=m.fs.fuelheater.inlet)
    m.fs.s02 = Arc(source=m.fs.fuelheater.outlet, destination=m.fs.stack.fuel_inlet)
    #m.fs.s03 = Arc(source=m.fs.oxygenfeed.outlet, destination=m.fs.stack.oxygen_inlet)
    m.fs.s04 = Arc(source=m.fs.stack.fuel_outlet, destination=m.fs.fuelheat_exchanger.hot_side_inlet)
    
    m.fs.s03 = Arc(source=m.fs.oxygenfeed.outlet, destination=m.fs.oxygenheater.inlet)
    m.fs.s05 = Arc(source=m.fs.oxygenheater.outlet, destination=m.fs.stack.oxygen_inlet)

    pyo.TransformationFactory("network.expand_arcs").apply_to(m)
    
    fix_cell_parameters(m.fs.stack.solid_oxide_cell)
    m.fs.stack.number_cells.fix(10)
    
    #set the initial conditions for oxygen side
    m.fs.oxygenfeed.flow_mol[0].fix(cccm_to_mps(3500))
    m.fs.oxygenfeed.mole_frac_comp[0, "O2"].fix(0.21)
    m.fs.oxygenfeed.mole_frac_comp[0, "N2"].fix(0.79)
    m.fs.oxygenfeed.temperature[0].fix(293)
    #m.fs.oxygenfeed.temperature[0].fix(1073)
    m.fs.oxygenfeed.pressure.fix(P)
    
    m.fs.oxygenheater.outlet.temperature.fix(1073)
    
    T_in = 800 + 273.15
    N_H2 = cccm_to_mps(410)
    N_H2O = cccm_to_mps(410)
    y_H2 = 0.5
    y_H2O = 0.5

    #caculate the cell's Nernst potential
    dHrxn = (
        _comp_enthalpy_expr(T_in, "H2")
        + 0.5 * _comp_enthalpy_expr(T_in, "O2")
        - _comp_enthalpy_expr(T_in, "H2O")
    )
    dSrxn = (
        _comp_entropy_expr(T_in, "H2")
        + 0.5 * _comp_entropy_expr(T_in, "O2")
        - _comp_entropy_expr(T_in, "H2O")
    )
    dGrxn = dHrxn - T_in * dSrxn

    V_nernst = pyo.value(
        1
        / (2 * _constF)
        * (
            dGrxn
            - 8.31
            * T_in
            * pyo.log(y_H2O / (y_H2 * 0.21**0.5) * (P / 101300) ** -0.5)
        )
    )
    
       
    m.fs.fuelfeed.temperature[0].fix(293)
    m.fs.fuelfeed.pressure.fix(P)
    m.fs.fuelfeed.flow_mol[0].fix(N_H2+N_H2O)
    m.fs.fuelfeed.mole_frac_comp[0, "H2"].fix(y_H2)
    m.fs.fuelfeed.mole_frac_comp[0, "H2O"].fix(y_H2O)
    
    m.fs.fuelheat_exchanger.overall_heat_transfer_coefficient[0].fix(500) # W/m2/K
    m.fs.fuelheat_exchanger.shell_outlet.temperature.fix(300)  
    
    
    m.fs.fuelheater.outlet.temperature.fix(T_in)
       
    
    m.fs.stack.potential_cell[0].fix(V_nernst+0.2)
    m.fs.stack.solid_oxide_cell.recursive_scaling()
    
    DOF = degrees_of_freedom(m)
    print(f"DOF before initialization is ({DOF})")
        
    seq = SequentialDecomposition()
    seq.options.select_tear_method = "heuristic"
    seq.options.tear_method = "Wegstein"
    seq.options.iterLim = 5
      
    G = seq.create_graph(m)
    heuristic_tear_set = seq.tear_set_arcs(G, method="heuristic")
    order = seq.calculation_order(G)
   
    for o in heuristic_tear_set:
        print(o.name) 
   
    for o in order:
        print(o[0].name)
        
    tear_guesses = {
            "mole_frac_comp": { 
                    (0, "H2"): 0.5,
                   (0, "H2O"): 0.5},
            "temperature": {0: 1000},
            "pressure": {0: 85910},
            "flow_mol":{0:0.0006}}

    # Pass the tear_guess to the SD tool
    seq.set_guesses_for(m.fs.fuelheater.inlet, tear_guesses)

    def function(unit):
            unit.initialize(outlvl=idaeslog.INFO)
 
    seq.run(m, function)
    
    DOF = degrees_of_freedom(m)
    print(f"DOF after initialization is ({DOF})")

    #assert degrees_of_freedom(m) == 0
    
    solver = get_solver("ipopt")
    status = solver.solve(
        m,
        tee=False,
        options={
            "tol": 1e-6,
            "max_iter": 100,
            "nlp_scaling_method": "user-scaling",
        },
    )
    
    #assert (
        #status.solver.termination_condition == pyo.TerminationCondition.optimal
    #)
    #assert status.solver.status == pyo.SolverStatus.ok
    
    #The following codes are used for show RSOC's calculation results    
    cell = m.fs.stack.solid_oxide_cell
    fuel_out = cell.fuel_outlet
    oxygen_out = cell.oxygen_outlet
    tracked_vars = {
        "voltage": cell.potential[0],
        "average_current_density": sum(cell.current_density[0, j] for j in cell.iznodes)
        / len(cell.iznodes),
        "average_temperature_z": sum(cell.temperature_z[0, j] for j in cell.iznodes)
        / len(cell.iznodes),
        "fuel_out_flow_mol": fuel_out.flow_mol[0],
        "fuel_out_temperature": fuel_out.temperature[0],
        "fuel_out_pressure": fuel_out.pressure[0],
        "fuel_out_mol_frac_comp[H2]": fuel_out.mole_frac_comp[0, "H2"],
        "fuel_out_mol_frac_comp[H2O]": fuel_out.mole_frac_comp[0, "H2O"],
        "oxygen_out_flow_mol": oxygen_out.flow_mol[0],
        "oxygen_out_temperature": oxygen_out.temperature[0],
        "oxygen_out_pressure": oxygen_out.pressure[0],
        "oxygen_out_mol_frac_comp[O2]": oxygen_out.mole_frac_comp[0, "O2"],
        "oxygen_out_mol_frac_comp[N2]": oxygen_out.mole_frac_comp[0, "N2"],
    }

    out={}
        
    for key in tracked_vars.keys():
        out[key] = pyo.value(tracked_vars[key])

[Hanli131]

Hi again,

I’ve done some further testing and discovered a new situation that might be helpful.

If I do not include the heat exchanger loop, the system can initialize and run successfully even with the oxygenheater present. This setup is shown in Figure 3.

In this case, I manually applied the stack outlet conditions to the hot side of the heat exchanger (see lines 288–292 in the code), instead of forming a full loop. With this adjustment, the degrees of freedom remain at 0 after initialization. The code is included below.
SOC_system_without_heat _loop.zip

Just wanted to share this additional observation in case it helps identify the root cause.

[dallan-keylogic]

I used the report_structural_issues and display_overconstrained_set methods of the IDAES diagnostics toolbox in order to determine the issue:

Dulmage-Mendelsohn Over-Constrained Set

    Independent Block 0:

        Variables:

            fs.oxygenheater.control_volume.properties_in[0.0].flow_mol
            fs.oxygenheater.control_volume.properties_in[0.0].flow_mol_phase[Vap]
            fs.oxygenheater.control_volume.properties_in[0.0].mole_frac_phase_comp[Vap,O2]
            fs.oxygenheater.control_volume.properties_out[0.0].flow_mol_phase[Vap]
            fs.oxygenheater.control_volume.properties_out[0.0].mole_frac_phase_comp[Vap,O2]
            fs.oxygenheater.control_volume.properties_in[0.0].mole_frac_comp[O2]
            fs.oxygenheater.control_volume.properties_out[0.0].flow_mol
            fs.oxygenheater.control_volume.properties_out[0.0].mole_frac_comp[O2]
            fs.oxygenheater.control_volume.properties_out[0.0].mole_frac_comp[N2]
            fs.oxygenheater.control_volume.properties_in[0.0].mole_frac_comp[N2]
            fs.oxygenheater.control_volume.properties_in[0.0].mole_frac_phase_comp[Vap,N2]
            fs.oxygenheater.control_volume.properties_out[0.0].mole_frac_phase_comp[Vap,N2]

        Constraints:

            fs.oxygenheater.control_volume.properties_in[0.0].total_flow_balance
            fs.oxygenheater.control_volume.material_balances[0.0,O2]
            fs.oxygenheater.control_volume.properties_in[0.0].component_flow_balances[O2]
            fs.oxygenheater.control_volume.properties_out[0.0].total_flow_balance
            fs.oxygenheater.control_volume.properties_out[0.0].component_flow_balances[O2]
            fs.s03_expanded.mole_frac_comp_equality[0.0,O2]
            fs.s05_expanded.flow_mol_equality[0.0]
            fs.oxygenheater.control_volume.properties_out[0.0].sum_mole_frac_out
            fs.s05_expanded.mole_frac_comp_equality[0.0,N2]
            fs.oxygenheater.control_volume.properties_in[0.0].component_flow_balances[N2]
            fs.oxygenheater.control_volume.material_balances[0.0,N2]
            fs.oxygenheater.control_volume.properties_out[0.0].component_flow_balances[N2]
            fs.s03_expanded.flow_mol_equality[0.0]
            fs.s03_expanded.mole_frac_comp_equality[0.0,N2]
            fs.s05_expanded.mole_frac_comp_equality[0.0,O2]

    Independent Block 1:

        Variables:

            fs.oxygenheater.control_volume.properties_in[0.0].pressure
            fs.oxygenheater.control_volume.properties_out[0.0].pressure

        Constraints:

            fs.oxygenheater.control_volume.pressure_balance[0.0]
            fs.s05_expanded.pressure_equality[0.0]
            fs.s03_expanded.pressure_equality[0.0]

    Independent Block 2:

        Variables:


        Constraints:

            fs.s05_expanded.temperature_equality[0.0]

Independent Block 2 was the easiest to look at, because it has a single constraint and zero variables. That constraint was the temperature equality constraint between the oxygen heater outlet, which was fixed as a design specification, and the solid oxide cell oxygen inlet. The oxygen inlet temperature is incorrectly fixed, as were all the other state variables on the port. Unfixing them should fix your immediate problem, but that doesn't explain why the issue occurred in the first place.

My first thought was that it might be a bug in the SOC initialization method, but we already have a test verifying that port variables are unfixed as expected. I didn't see any obvious reason in your code about why they might remain fix. Perhaps it's an issue in the Pyomo SequentialDecompositon?

[Hanli131]

Thank you for your reply!

I tried adding a heat exchanger loop to the oxygen side (so that the oxygen side now has a similar structure to the fuel side), and interestingly, the issue disappeared. The model can now maintain zero degrees of freedom both before and after initialization.

In the future, if I encounter similar cases, I think I will avoid using SequentialDecomposition and instead write custom initialization routines myself.

By the way, is there any tutorial or example in IDAES that explains how to write custom initialization procedures? I’d like to learn more about that. ：）

[dallan-keylogic]

Probably the best example is the initialization method for the SOC Flowsheet in Examples. It's not used in the public example but was used in a steady state optimization to generate the starting point and setpoints for dynamic optimization. (I hope to make an example notebook for that optimization public one day.)

[Hanli131]

Looking forward to your example notebook. 💐