Reactor initialization error with GenericReactionParameterBlock

[am000-simreka]

Dear Ideas community,
I am trying to reproduce the hydrodealkylation conversion reactor from IDAES-PSE Tutorial 3: using the generic property and reaction packages using an equation of state. However, I am facing issues in the initialization step.
Our goal is to test idaes reactor models with equations of state so in the future we can use different chemicals.
To reproduce the problem, you can run the attached had_reactors.py file. I am also attaching the had_config file that contains the configuration file for the property and reaction package.
Can you help us to troubleshoot? The diagnostic tool box points to the enthalpy balance.
hda_files.zip

[andrewlee94]

Hi @am000-simreka, the answer in this case is actually very simple (and a very common occurance), but it is not what you think it is. As you noted you have a number of potential evaluation errors, but those are not actually causing problems for you. Unfortunately, many thermodynamic correlations have potential singularities (log of mole fractions, inverse relative temperatures, etc.) but as long as you are not operating near those points, you can use bounds to prevent them from being displayed.

However, the first thing you should do in the case of any solve (failed or not) is look at the solver logs, and in this case the termination message from the last solver step:

"Initialization Step 2 infeasible - Converged to a locally infeasible point. Problem may be infeasible..".

This tells you that is issue is that the problem is quite possibly infeasible (or very poorly scaled) meaning you need to look to the numerical diagnostics. To do this however, you will need to wrap the initialization in a try/except to get past the failure error (also for reference, I strongly recommend moving to the new Initializer approach for initialization as it helps avoid some of these issues due to better error handling). I replaced the last part of your script with the following:

Diagnostic1 = DiagnosticsToolbox(m)
Diagnostic1.report_structural_issues()
# Diagnostic1.display_potential_evaluation_errors()

try:
    m.fs.R101.initialize(outlvl=idaeslog.DEBUG)
except:
    pass

Diagnostic1.report_numerical_issues()
Diagnostic1.display_variables_at_or_outside_bounds()
Diagnostic1.compute_infeasibility_explanation()
#
# solver = SolverFactory('ipopt')
# status = solver.solve(m, tee=True) # tee outputs solver information.
# m.fs.R101.report()

If you look at the results of the numerical diagnostics, there is an obvious place to start - a variable at its bounds. If you then display this variable you will see that it is the outlet temperature and that it is sitting at the upper bound you have set. If you then look at the output of the infeasibility explainer, it tells you:

"Model unknown may be infeasible. A feasible solution was found with only the following variable bounds relaxed: ub of var fs.R101.control_volume.properties_out[0.0].temperature by 8.000313597318051"

I.e., the problem would be feasible if you relaxed the upper bound by a little over 8K, which I confirmed by changing the bound to 610 at which point the model solved. However, before doing this you should stop and ask why you set that bound in the first place, as it is often due to the range over which one or more correlations were fitted and thus relaxing the bound means extrapolating with those correlations. You need to decide if relaxing the bound is acceptable, or whether it means you should be looking for a different set of correlations to use.