[Observation] Constraints on Feed Units May Not Take Effect During Initialization

[Hanli131]

I’d like to share some thoughts based on a simple example I created. The system diagram is shown below, and the objective is to calculate the outlet temperature.

In the model, I applied two constraints: one to fix the flow rate of H2feed and another to fix the outlet temperature of H2heater.

The mixer report before initialization is shown below:
As you can see, the specified values were not passed to the mixer before initialization — all parameters remain at their initial default values.

The report after initialization is as follows:

At this stage, the outlet temperature of H2heater has been successfully constrained. However, the flow rate of H2feed has not been constrained.

The report after solving the model is as follows:

Here, both constraints are correctly applied, and the flow rate of H2feed is as expected. This indicates that the constraint on H2feed itself is correctly defined but was not applied during initialization.

My guess is that because H2feed is located at the beginning of the system, it was not included in the initialization process. In this simple case, the abnormal flow value at inlet2 after initialization does not affect the final results. However, in more complex systems, such incorrect initialization values could lead to convergence failure or even model crashes (which I have already encountered in practice).

Suggested Solution:
Adding the following code during the initial parameter setup can assign a reasonable initial value to H2feed:

m.fs.H2feed.flow_mol.fix(10)
m.fs.H2feed.flow_mol.unfix()

Summary:
Constraints applied to feed units may not take effect during the initialization process, which can potentially cause issues in later calculations. Therefore, it is important to pay special attention to this when building models.

Share the codes here:

import pyomo.environ as pyo
from pyomo.environ import Objective, maximize
from pyomo.environ import Constraint, value
from idaes.core import FlowsheetBlock
from idaes.core.util.model_statistics import degrees_of_freedom
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.unit_models import Mixer, MomentumMixingType, Feed
from idaes.models.unit_models.heater import Heater
import idaes.logger as idaeslog

from pyomo.network import Arc, SequentialDecomposition

fuel_comps = ["H2", "H2O"]
oxygen_comps = ["O2", "N2"]

m = pyo.ConcreteModel()
m.fs = FlowsheetBlock(dynamic=False, time_set=[0], time_units=pyo.units.s)   

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"
)

#add units
m.fs.waterfeed = Feed(property_package=m.fs.fuel_params)

m.fs.waterfeed2 = Feed(property_package=m.fs.fuel_params)

m.fs.waterheater = Heater(property_package=m.fs.fuel_params)

m.fs.H2feed = Feed(property_package=m.fs.fuel_params)

m.fs.H2heater = Heater(property_package=m.fs.fuel_params)

m.fs.fuelmixer = Mixer(property_package=m.fs.fuel_params,
                     num_inlets=3,
                     momentum_mixing_type=MomentumMixingType.minimize)


#add arcs
m.fs.wf2wh = Arc(source=m.fs.waterfeed.outlet, destination=m.fs.waterheater.inlet)

m.fs.wh2fm = Arc(source=m.fs.waterheater.outlet, destination=m.fs.fuelmixer.inlet_1)

m.fs.hf2hh = Arc(source=m.fs.H2feed.outlet, destination=m.fs.H2heater.inlet)

m.fs.hh2fm = Arc(source=m.fs.H2heater.outlet, destination=m.fs.fuelmixer.inlet_2)

m.fs.wf22fm = Arc(source=m.fs.waterfeed2.outlet, destination=m.fs.fuelmixer.inlet_3)

pyo.TransformationFactory("network.expand_arcs").apply_to(m)


#set initial inputs
def flow_mol_rule(b):
    return m.fs.H2feed.flow_mol[0] == m.fs.waterfeed.flow_mol[0]/10

m.fs.H2feed.flow_mol_con = Constraint(rule=flow_mol_rule)

m.fs.waterfeed.flow_mol[0].fix(100)
m.fs.waterfeed.pressure[0].fix(101325)
m.fs.waterfeed.temperature[0].fix(400)
m.fs.waterfeed.mole_frac_comp[0, "H2"].fix(1e-6)    
m.fs.waterfeed.mole_frac_comp[0, "H2O"].fix(1-1e-6)

m.fs.waterfeed2.flow_mol[0].fix(90)
m.fs.waterfeed2.pressure[0].fix(101325)
m.fs.waterfeed2.temperature[0].fix(1000)
m.fs.waterfeed2.mole_frac_comp[0, "H2"].fix(1e-6)    
m.fs.waterfeed2.mole_frac_comp[0, "H2O"].fix(1-1e-6)

m.fs.H2feed.flow_mol.fix(10)
m.fs.H2feed.flow_mol.unfix()
m.fs.H2feed.pressure.fix(101325)
m.fs.H2feed.temperature[0].fix(400)
m.fs.H2feed.mole_frac_comp[0, "H2"].fix(1-1e-6)    
m.fs.H2feed.mole_frac_comp[0, "H2O"].fix(1e-6)

m.fs.waterheater.outlet.temperature.fix(550)
def H2_heater_rule(b):
    return m.fs.H2heater.outlet.temperature[0] == m.fs.waterheater.outlet.temperature[0] - 100

m.fs.H2heater.H2_heater_con = Constraint(rule=H2_heater_rule)


print(f"initial DOF is {degrees_of_freedom(m)}")
print("before the initialization:")
print(value(m.fs.H2feed.flow_mol[0]))
print(value(m.fs.waterfeed.flow_mol[0]))
m.fs.fuelmixer.report()


#initialization
seq = SequentialDecomposition()
seq.options.select_tear_method = "heuristic"
seq.options.tear_method = "Wegstein"
seq.options.iterLim = 500

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)
    
    
def function(unit):
        unit.initialize(outlvl=idaeslog.INFO)
 
seq.run(m, function)

print("after the initialization:")
print(value(m.fs.H2feed.flow_mol[0]))
print(value(m.fs.waterfeed.flow_mol[0]))
m.fs.fuelmixer.report()


#solve
solver = pyo.SolverFactory('ipopt')
results = solver.solve(
    m, 
    tee=False,
    options={
        "tol": 1e-8,
        "max_iter": 20000,
        "nlp_scaling_method": "user-scaling",
        "constr_viol_tol": 1e-6,
    })  

print("calculation results:")
print(value(m.fs.H2feed.flow_mol[0]))
print(value(m.fs.waterfeed.flow_mol[0]))
m.fs.fuelmixer.report()

[dallan-keylogic]

In general, model-specific initialization methods consider only the constraints created by the unit/property model and not additional ones added by the user. Adding additional constraints on the unit model before initialization can result in too few degrees of freedom. You should wait until after initialization to add any additional constraints that you have.