1D Reactor Models in IDAES

[denbigh]

Hi,

I'm interested in developing custom, 1D reactor models (e.g. steady state catalytic packed bed reactor with heat- and mass-balances, 1D absorber models, 1D+time unsteady state pressure swing adsorption models, etc.) I would normally use specialized DE solvers in scipy or DifferentialEquations.jl to solve these systems (via Method of Lines, or specialized IVP or BVP solvers, depending on the problem) but our project is considering using the IDAES framework, due to the integration with pyomo. In the examples I have seen, the custom models still follow quite a 'rigid' form; they don't allow you to define an arbitrary python function which utilizes complex packages (e.g. scipy) and provides a vector output. Is this impression correct? Is the IDAES framework flexible enough to handle 1D ODE problems solved with arbitrary solvers? Even a simple example with a 1D plug flow reactor would be very useful. Thanks!

[kanishka-ghosh]

The plug flow reactor unit model (link) in IDAES is a 1D model that can be used for steady-state and dynamic reactor modeling with ODEs. The ODEs are solved as IVPs using solver IPOPT.

Following are some tests/examples that show how to use the PFR model:
https://github.com/IDAES/idaes-pse/blob/main/idaes/models/unit_models/tests/test_pfr.py
https://github.com/IDAES/examples-pse/blob/main/src/Examples/UnitModels/Reactors/plug_flow_reactor_testing.ipynb

[andrewlee94]

We also have tools to allow you to use integrators such as those in the PETSc solver libraries to solve ODEs: https://github.com/IDAES/examples-pse/tree/main/src/Examples/DAE

[denbigh]

Hi all,

Thanks for the replies! I imagine that an arbitrary black-box unit operation model (i.e. a python script that defines a function mapping inputs to outputs) is not possible, as that wouldn't play well with Pyomo? These options do seem quite flexible though.
Thanks!

[denbigh]

Hi all,

Yes, this was my understanding. I suspect that the surrogate model option may be a viable alternative if the ODEs do not conform to the models pre-built into IDAES.

Thanks

[andrewlee94]

Note that IDAES tried to keep the standards/requirements to a minimum - fundamentally the end result is a Pyomo model and the IDAES tools are just there to try to help you along the way. so you do not need to conform to the pre-built models in anyway.

However, I will add one final note (at the risk of sending you off somewhere else) - you should always match the tool you are using to the problem you are trying to solve. IDAES and Pyomo are fundamentally built with optimization in mind; if all you are trying to do is run simulations of the model, then a stand-alone integrator will be easier to use and more reliable. However, if you ever want to optimize a model (especially one with ODEs) - that is when you need tools like Pyomo and IDAES. Building on that; ODEs are inherently difficult to work with, so you should always ask if a surrogate will be sufficient for your needs.

[denbigh]

@andrewlee94 Just following up on this - I've been playing with some of the surrogate modeling tools in IDAES (e.g. ALAMO) and have been impressed with them. However, the examples here,

https://idaes-pse.readthedocs.io/en/1.11.0/user_guide/modeling_extensions/index.html

mainly focus on building surrogate models. I don't think there are any examples of how to use a surrogate model within IDAES. For example, if you use ALAMO(PY) to build a simple surrogate of a chemical reactor, how would you incorporate this as a unit operation in a larger simulation? Are there any examples of this work flow in the docs?

Thanks!

[andrewlee94]

Hi,

There are a few examples on using surrogates in a flowsheet here: https://idaes-examples.readthedocs.io/en/latest/docs/surrogates/alamo/alamo_flowsheet_optimization_doc.html

The SkeletonModel is also intended for this purpose, but I don;t think we have a good example that combines it with a surrogate model object. There is this example, but the "surrogate" used there an empirical correlation: https://idaes-examples.readthedocs.io/en/latest/docs/unit_models/operations/skeleton_unit_doc.html

[denbigh]

Thanks, I wasn't aware of these tutorial notebooks!

[denbigh]

@adowling2 This surrogate model idea looks interesting, though I haven't used them much myself. If I have an arbitrary black box model in python that solves quite quickly (say < 1s), and takes a few scalars as input and spits a few scalars out as outputs, would it be reasonable to train a surrogate model that can mimic this model & is usable in IDAES? Honestly, I'd be happy with a spline fit if that were usable in IDAES :) There seem to be a lot of surrogate modeling tools in IDAES (ALAMO, PySMO, Keras, etc.) Is there one in particular that would work for this application?

[tk3016]

Hi all,
I have a follow-up question based on the previous discussion on this thread.
I managed to model a 1D reactor using a PFR. However, I would like to add custom energy and momentum balance equations.
Is it possible to do that for the PFR unit in IDAES, or will I have to create a custom unit model similar to the PFR?

Best wishes,
Tanuj

[andrewlee94]

When constructing the PFR, you can set energy_balance_type: EnergyBalanceType.none and momentum_balance_type: MomentumBalanceType.none (importing EnergyBalanceType and MomentumBalanceType from idaes.core). This will tell the PFR (and the underlying 1D control volume not to add any energy or momentum balance equations so that you can then go and add your own.

Note that adding your own energy and momentum balances may require you to write your own initialization routine as well (although there is a good chance the existing one will work).

Alternatively, it is is possible to turn off the default energy and momentum balances if you know their names (generally unit.control_volume.enthalpy_balances.deactivate() and unit.control_volume.pressure_balance.deactivate()); this lets you switch between the default terms and any custom constraints you want which might be helpful for initialization and/or testing.

[tk3016]

Dear Andrew,

As always, many thanks for your prompt response.
After looking through the documentation and other discussion threads in more detail, I figured that writing custom momentum/energy balances was not a wise thing to do. Instead, what I did was to define the deltaP term for momentum balance by defining additional constraints after setting the argument has_pressure_change = True and it worked well. I tried to do the same for energy balance, where I wanted to model a cooling jacket around the reactor. I defined the appropriate component enthalpy coefficients and heat of reaction terms while setting up the property package and then defined a constraint to specify heat_duty for each element in the PFR to be equal to UA(T_c - T(z)), here T_c is the temperature of the coolant fluid and T(z) is the temperature of the reactor at length z. On adding this constraint, I am unable to initialise my reactor but without this my reactor works fine under isothermal/adiabatic conditions. I wanted some further clarity if this is the right approach to model this or am I doing something wrong?

Looking forward to your response.

Best wishes,
Tanuj

[andrewlee94]

Hi Tanuj,

That is the way I would do it as well - the heat_duty, work (if required) and deltaP terms are intended to be used like that. Regarding initialization, it depends a bit on what the cause of the failure is. If you are seeing the solver report infeasible solutions, then check the bounds on your variables (especially T) to make sure your heat transfer is not causing them to be violated.

A few things you can try:

1. First of all, run your model throguh the diagnostics tools to look for any modeling issues (e.g. structural singularities): https://idaes-pse.readthedocs.io/en/stable/explanations/model_diagnostics/index.html. This will help make sure your model is actually solvable (i.e. separate a modeling error from just initialization). You can also try to use the toolbox on the model after the failed initialization to see what it can tell you about why it failed.
2. Next, try using this https://idaes-pse.readthedocs.io/en/stable/reference_guides/initialization/bt_initializer.html to initialize the model (I assume you are currently doing m.fs.pfr.initialize()). This is an automatic initializer that works for a lot of simpler models so it is a good place to start. If your Tc is a constant value, or you have co-current flow this will hopefully work. However, it will struggle with counter-flow systems.
3. If the above does not work, the you can try a multi-step initialization approach, where you first solve with a fixed heat_duty (using some initial guess) with your heat transfer constraint deactivated. The BTInitializer or default PFR initialization routine should solve this (assuming your heat duty does not try to push the temperature outside its bounds). Once this solves, then unfix the heat duty and activate your constraint and re-solve (just using solver.solve(), do not call the initialization again).

[tk3016]

Hi Andrew,

Thanks.
I tried the third suggestion and it worked! :) I did not know one could initialise like this - good to know!
I was sure that the model would simulate because I had developed a stand-alone pyomo model for my reactor previously and it worked without any issues.
Many thanks again for your help.

Best wishes,
Tanuj