Design of SOFC model based on SOEC unit present in idaes, issue with DOF

[slargher]

Hi everyone! We are trying to design a SOFC model on IDAES using the SoecDesign present here in the repository under power generation units (idaes/models_extra/power_generation/unit_models/soec_design.py). Is anyone else working on something similar or did in the past? We are encountering infeasability errors in the model convergence and think we might have a scaling factor issue. Also the SOEC model when running has DOF=1, is that correct? From my understanding the DOF should be 0 for the unit models present in the SOFC but DOF = 1 or the whole flowsheet? It is possible that infeasability errors arise from just the scaling factors and not from errors in the Constraints? Thank you very much in advance for any help and I am open to discuss further about SOFCs!

[dallan-keylogic]

I don't have personal experience with the SoecDesign unit model. Although the current tests for that model do not assert that there are zero degrees of freedom, I was able to verify that there are, in fact, zero degrees of freedom when solving it. I would recommend starting from the flowsheet provided in that test, solving it as-is, then fixing/unfixing variables to suit your purpose.

Bad scaling can certainly result in infeasibility when the underlying model is feasible. If you have tiny mole fractions of certain components, IPOPT may think a bound is active when the mole fraction takes a small-but-nonzero value. Also, the enthalpy balance contains numbers many orders of magnitude greater than anything else in the problem, which makes it hard to get a residual small enough to meet IPOPT's tolerance without scaling. For a brief overview of scaling, check out Section 2 of this paper.

[slargher]

Thanks a lot for your reply!! We fixed the unfeasibility issue and now we are able to run and find an optimal solution on our sofc model based on the soec current tests. But we are encountering a weird behaviour of the solver, we get a different result (actual output values) every time we run the code and restart the kernel. Even changing the path of the .py file will cause a big difference between runs. We are using IPOPT, so the seeding is not the issue here. This does not happen with the models present on the IDAES repo so we are missing for sure some important characteristic of the solver behavioiur. Any idaes? :)

[dallan-keylogic]

Did you manage to resolve your degrees of freedom issue? I'm not sure how IPOPT acts with nonzero degrees of freedom and no objective function.

[slargher]

Hi! Thanks for your reply! Yes we did manage to get DOF = 0 for the whole model (to check we use print(degrees_of_freedom(model))), but when we print the DOF for each unit model those are not zero, and vary for each one. Is that a normal behaviour in you opinion? DOF = 0 should be for the whole model?
I don't want to be too forward or take advantage of your help, but I am not sure if I can explain myself clearly enough so if it is easier the code I am referring to is in my public repo (SOFC_uni.py for the model and SOFC_run_script.ipynb) Here we are still having some issues with the mixer as we can't manage to make it properly work as an ideal mixer conserving properties like the number of moles. We would like for the mixer to a sort of average of the input streams temperatures. Apart from the mixer, we think we are using standard constraints so we are not sure if we are missing some important equation.
Any help is appreciated, but of course, mindful of your time.

[dallan-keylogic]

when we print the DOF for each unit model those are not zero, and vary for each one. Is that a normal behaviour in you opinion?

This is normal. Those degrees of freedom correspond to the variables on the inlet ports of each unit model and are taken care of by arc constraints at the flowsheet level.

        self.i = Var(
            self.flowsheet().time, 
            initialize=0, 
            bounds=(0, 0.4),
            units=pyo.units.ampere / (pyo.units.m**2)/10000, 
            doc="Current density",
        )
...
        self.cell_area = Var(
            domain=pyo.NonNegativeReals,
            units= (pyo.units.m**2)/10000, 
            doc="Single-cell active area",
        )

Note that 1) Pyomo supports centimeters as a unit, 2) I don't know how it will react to having numbers in the units field, and 3) you've misplaced parentheses in the units for i. Be careful when using units that are not SI base units in equations, but it should be fine in this case because both variables are using cm**2.

        self.separator.split_fraction[:, "o2_rich_strm", "O2"].fix(1)
        self.separator.split_fraction[:, "o2_rich_strm", "N2"].fix(0)

        self.reactor_separator.split_fraction[:, "o2_strm", "O2"].fix(1)
        self.reactor_separator.split_fraction[:, "o2_strm", "H2"].fix(0)
        self.reactor_separator.split_fraction[:, "o2_strm", "H2O"].fix(0)

Despite it being present in the version of the model in the repo, fixing Vars during unit model construction is discouraged because it causes issues if Pyomo adds additional points to the time set during, for example, discretization for dynamics. Since this is a steady-state model, it should be fine here, but I thought I should point it out for future reference.

        @self.Constraint(time)
        def dT_upper(b, t):
            return ( b.mixer_out.mixed_state[t].temperature
                - b.translator_h2_out.properties_out[t].temperature ) <= 10

        @self.Constraint(time)
        def dT_lower(b, t):
            return ( b.mixer_out.mixed_state[t].temperature
                - b.translator_h2_out.properties_out[t].temperature ) >= -10 

Generally we like to keep performance constraints like this at the flowsheet level rather than having them exist at the unit model level.

        for pp in [
            self.config.anode_side_prop_package,     # H2/H2O
            self.config.cathode_side_prop_package,   # O2/N2 (air)
            self.reaction_props,         # common reaction basis
        ]:
            pp.set_default_scaling("flow_mol", 1e1)
            pp.set_default_scaling("flow_mol_phase", 1e1)

You may need to scale "flow_mol_phase_comp" as well.

the mixer as we can't manage to make it properly work as an ideal mixer conserving properties like the number of moles.

From looking at the Jupyter notebook, IPOPT is still returning infeasible:

WARNING: Loading a SolverResults object with a warning status into
model.name="unknown";
    - termination condition: infeasible
    - message from solver: Ipopt 3.13.2\x3a Converged to a locally infeasible
      point. Problem may be infeasible.

That would explain why the mixer isn't conserving material---the equations aren't being solved to the prescribed tolerance. As for why you might be having infeasibility, the current density you want to use, 0.3939 A/cm^2, is close to the upper bound of 0.4. That's close enough that IPOPT may mistake it as an active constraint. You might decrease fuel utilization or loosen that bound. The delta T bounds might also be causing trouble, so try deactivating them.
Otherwise, you can investigate the DiagnosticsToolbox to determine why your model is infeasible.

There are more sophisticated ways of determining exactly which variable bounds or constraints IPOPT think are active by looking at the dual variables, but presently Pyomo makes it something of a pain to get them. Streamlining that process is on the list of future additions to the diagnostics toolbox, though.

[slargher]

Thank you very much for the very explanative and fast reply! We implemented your suggestions in SOFC_unit.py and were able to fix the mixer behaviour. Now our Jupyter Notebook (SOFC_run_script.ipynb) returns an EXIT: Optimal Solution Found. but we are not sure why this happens if Model Statistics Degrees of Freedom: 1 When we investigate the Diagnostic tool:
`Dulmage-Mendelsohn Under-Constrained Set

Independent Block 0:

    Variables:

        fs.sofc.translator_h2_out.properties_in[0.0].temperature
        fs.sofc.reactor_separator.water_strm_state[0.0].temperature
        fs.sofc.translator_h2_out.properties_out[0.0].temperature
        fs.sofc.reactor_separator.mixed_state[0.0].temperature
        fs.sofc.reactor.control_volume.properties_out[0.0].temperature
        fs.sofc.reactor_separator.o2_strm_state[0.0].temperature
        fs.sofc.reactor.control_volume.heat[0.0]
        fs.sofc.translator_o2_out.properties_in[0.0].temperature
        fs.sofc.heater.control_volume.heat[0.0]
        fs.sofc.translator_o2_out.properties_out[0.0].temperature
        fs.sofc.heater.control_volume.properties_out[0.0].temperature
        fs.sofc.mixer_out.o2_strm_final_state[0.0].temperature
        fs.sofc.mixer_out.o2_poor_strm_state[0.0].temperature
        fs.sofc.mixer_out.mixed_state[0.0].temperature

    Constraints:

        fs.sofc.separator_fuel_to_translator_expanded.temperature_equality[0.0]
        fs.sofc.translator_h2_out.temperature_eqn[0.0]
        fs.sofc.reactor_separator.temperature_equality_eqn[0.0,water_strm]
        fs.sofc.react_to_sep_expanded.temperature_equality[0.0]
        fs.sofc.reactor_separator.temperature_equality_eqn[0.0,o2_strm]
        fs.sofc.reactor.control_volume.enthalpy_balances[0.0]
        fs.sofc.separator_o2rich_to_translator_expanded.temperature_equality[0.0]
        fs.sofc.sofc_energy_balance[0.0]
        fs.sofc.translator_o2_out.temperature_eqn[0.0]
        fs.sofc.heater.control_volume.enthalpy_balances[0.0]
        fs.sofc.mix_o2_final_expanded.temperature_equality[0.0]
        fs.sofc.heater_to_mixer_out_expanded.temperature_equality[0.0]
        fs.sofc.mixer_out.enthalpy_mixing_equations[0.0]` 

our problem seems to lie in the temperature variables, woud it make sense to .fix a possible outlet temperature of the outlet streams and then .unfix them? We saw it is something that was done in the SoecDesign.py:
self.oxygen_side_inlet.fix() self.hydrogen_side_inlet.fix() self.oxygen_side_outlet.unfix() self.hydrogen_side_outlet.unfix()

and like this in its Jup Notebook:

m.fs.soec.hydrogen_side_outlet_temperature.fix(1023) m.fs.soec.oxygen_side_outlet_temperature.fix(1023)

Our aim would be to create a model that gives us the values of the outlet temperatures, therefore we are having a hard time understanding the logic of doing so. Would an upper boundary for the outlets' temperature work in the same way?
We initially thought that these value would derive form the internal energy balance, but at the moment we are trying to understand why our heat_reaction = 0 and heat_stoichiometric_reactor != 0, is there a way to inspect these equations?

Reactor - Inlet                       

STREAM PROPERTIES:
Flow [mol/s]        : 1.803970e-01
Temperature [K]     : 945.99
Pressure [Pa]       : 101325.00
Heat [W]            : -11309.60
Heat of reaction[W] : -0.00


           
 Reactor - Outlet                      

STREAM PROPERTIES:
Flow [mol/s]        : 1.542719e-01
Temperature [K]     : 1230.27
Pressure [Pa]       : 101325.00
Heat [W]            : -11309.60
Heat of reaction[W] : -0.00

We appreciate any help, but we also know we are asking for a lot of input and feedback. As this is our first time working with IDAES we are still very much in the learning process. Thanks a lot for the suggestions so far, they have been incredibly helpful.

[dallan-keylogic]

but we are not sure why this happens if Model Statistics Degrees of Freedom: 1

IPOPT can solve underconstrained problems to find a feasible point. However, if you're missing an equation, the solution may not be physically meaningful.

our problem seems to lie in the temperature variables, woud it make sense to .fix a possible outlet temperature of the outlet streams and then .unfix them? We saw it is something that was done in the SoecDesign.py: ... Our aim would be to create a model that gives us the values of the outlet temperatures, therefore we are having a hard time understanding the logic of doing so.

The SoecDesign was never designed to be predictive, but just something quick to handle material and energy balances until actually predictive models, the SolidOxideCell and SolidOxideModuleSimple, were created. I don't fully understand what's going on inside the SoecDesign in terms of energy balances at the moment.

In terms of your energy balances, you might try adding a constraint equating the fuel side and oxygen side temperatures. Hopefully that works, check the Dulmage-Mendelsohn partition to make sure your model isn't structurally singular (which is a necessary but not sufficient condition for your model to give physically meaningful results).

[slargher]

Thanks a lot for your reply! We kept working on our model and implement your suggestions, so this Q&A has been super helpful. We are trying to study now how the stoichiometric reactor works for two main reasons:

1. our Sofc energy balance (when double checked with another script) still doesn't close and is not 0, we think the issue lies in how we deal with the heat variables:

    self.reactor = um.StoichiometricReactor(
        property_package             = self.reaction_props,
        reaction_package             = self.rxn_props,
        has_pressure_change          = False,
        has_heat_transfer            = True,
        has_heat_of_reaction         = False,
    )
   

We assume our reactor is adiabatic, it doesn't lose heat to the environment and the heat created from the esothermic reaction is entirely passed to the reaction products so generally we would put has_heat_of_reaction = True and has_heat_transfer = False. But the heat variable is only created when has_heat_transfer = True.
In your opinion, are we giving the arguments has_heat_transfer and has_heat_of_reaction the right meaning for our purpose?

2. Now that we have a H2-Sofc we would like to have NH3 as the fuel (new packages already created), this means that our stoichiometric reactor would need to process two different reactions: NH3-cracking which produces H2, and following H2-combustion as earlier, one after the other. In the documentation having multiple reactions seems possible, but we are having a hard time to find examples and figure out if sequential reactions are a possibility. Another option is to use one reactor after the other, in your experience, which one do you suggest is the "cleanest" way?

Any suggestion is appreciated, thank you for your time

[dallan-keylogic]

We assume our reactor is adiabatic, it doesn't lose heat to the environment and the heat created from the esothermic reaction is entirely passed to the reaction products so generally we would put has_heat_of_reaction = True and has_heat_transfer = False. But the heat variable is only created when has_heat_transfer = True.
In your opinion, are we giving the arguments has_heat_transfer and has_heat_of_reaction the right meaning for our purpose?

Are you assuming that the stoichiometric reactor object is adiabatic or that the entire SOC module is adiabatic? If the stoichiometric reactor is adiabatic, then you don't need a heat object and want has_heat_transfer=False. Passing has_heat_of_reaction creates an heat_of_reaction Expression on the ControlVolume0D that sums reactions[t].dh_rxn[r] for all the reactions defined in your reaction package. The value you need to use for dh_rxn, however, depends on what reference enthalpy you use for your components. If you include enthalpy of formation in your component enthalpy values, I believe dh_rxn should be zero and you can just pass has_heat_of_reaction=False.

In the documentation having multiple reactions seems possible, but we are having a hard time to find examples and figure out if sequential reactions are a possibility.

If you're using the modular property framework, you just add a second reaction with a second index. I don't think there's any specific example, but here I added a second (parallel) reaction for Methanol synthesis creating formaldehyde (with made up numbers for parameters).

    "rate_reactions": {
        "R1": {
            "stoichiometry": {
                ("Vap", "CO"): -1,
                ("Vap", "H2"): -2,
                ("Vap", "CH3OH"): 1,
            },
            "heat_of_reaction": constant_dh_rxn,
            "concentration_form": ConcentrationForm.moleFraction,
            "rate_constant": arrhenius,
            "rate_form": power_law_rate,
            "parameter_data": {
                "reaction_order": {("Vap", "CO"): 1, ("Vap", "H2"): 2},
                "arrhenius_const": (3.77287e19, pyunits.mol / pyunits.m**3 / pyunits.s),
                "dh_rxn_ref": (-90640, pyunits.J / pyunits.mol),
                "energy_activation": (109.2e3, pyunits.J / pyunits.mol),
            },
        }
        "R2": {
            "stoichiometry": {
                ("Vap", "CO"): -1,
                ("Vap", "H2"): -1,
                ("Vap", "CH2O"): 1,
            },
            "heat_of_reaction": constant_dh_rxn,
            "concentration_form": ConcentrationForm.moleFraction,
            "rate_constant": arrhenius,
            "rate_form": power_law_rate,
            "parameter_data": {
                "reaction_order": {("Vap", "CO"): 1, ("Vap", "H2"): 1},
                "arrhenius_const": (1.2345e6, pyunits.mol / pyunits.m**3 / pyunits.s),
                "dh_rxn_ref": (-98765, pyunits.J / pyunits.mol),
                "energy_activation": (314.1e5, pyunits.J / pyunits.mol),
            },
        }
    },

Reactions in series don't pose a problem writing equations, but you need to be careful that all flow rates remain positive when specifying extents of reaction.