Use of REFPROP as a property package

[eric-hintz]

(TLDR near bottom)
Hello,

I have been reading through the IDAES-PSE documentation (which is quite well done, by the way, thank you for the effort put into that!) to determine if IDAES would be a good fit for my modeling project. It seems like a really great software with a lot of features that I've been trying to implement manually just using SciPy, so I am very glad to find that there is software that solves these problems for me.

However, the one concern I have is the property package. Various sections of my modeling conditions include difficult to model thermodynamic conditions for certain components and complicated mixtures in other sections. The only thermodynamics package that resolves these conditions properly is REFPROP (or more accurately for the complicated mixtures, a pre-release beta version of REFPROP sent to me by the developers at NIST). However, REFPROP is distributed as Fortran code with a DLL and wrappers for various languages including python, making it a complete black box for modeling purposes. From what I gather reading the documentation and other discussions posted here about similar things (especially #1265), this poses a challenge for IDAES (or more accurately, Pyomo/IPOPT).

If I was only calculating properties for pure fluids, I could simply use the Helmholtz property package and load in the coefficients manually. However, the mixtures I am working with are extremely non-ideal, meaning it is very important that I use the Helmholtz mixing rules used by REFPROP (and GERG and others like them), with the accurate coefficients used by REFPROP. Here are my initial ideas to solve this issue based on what I've read:

1. Write a custom property package that essentially wraps calls to the REFPROP DLL in a function that also returns the first (and maybe second) derivatives using some automatic differentiation package for the derivatives not explicitly available from REFPROP (perhaps JAX to use its functional paradigm to effectively "cache" derivatives). However, this is likely to be slow and I'm not totally sure it's even possible.
2. Use the Helmholtz EOS property package and manually load in the coefficient values from REFPROP for the components that I will use. This is probably the method that best integrates into the IDAES ecosystem, and wouldn't be too difficult to implement (theoretically at least). However, the concern I have here is the mixtures I will need to deal with. If possible, I am happy to load in the mixing coefficients for all of the binary pairs in my model, but I am not sure that the mixing rules used by REFPROP and similar Helmholtz EOS are defined anywhere in the IDAES ecosystem.
3. Write a custom property package that defines the pure fluid and mixture EOS from scratch, and then load in the coefficients. I am not sure this has any advantage over the Helmholtz property model except for handling mixtures, but if I could just add on the mixing rules to the Helmholtz property package, that would probably save some significant effort.

For a little context, the conditions/mixtures I am modeling that require this effort are cryogenic hydrogen and ammonia-containing mixtures, specifically ammonia, hydrogen, and nitrogen. The overall model is an ammonia synthesis process. Or rather, many ammonia synthesis processes requiring a large range of thermodynamic states.

TLDR: I need to use REFPROP, or at least the Helmholtz EOS with REFPROP's coefficients, for some complex thermodynamic calculations including highly non-ideal mixtures, which appears problematic within the IDAES ecosystem.

Are there better ways to resolve this problem that I am not aware of? Or is this issue enough of a non-starter that I will not be able to use IDAES with this restriction?

[andrewlee94]

@eric-hintz There are a few option for what you want to do, each with their own advantages and disadvantages. The biggest problem you face here is that the Helmholtz equations of state themselves are high-order polynomials which require specialized root finding routines to give the right answer, which is why IDAES packaged them as external function calls (it is impractical to implement these as true equation-oriented constraints).

1. It you have at least first partial derivatives, you can use Pyomo's Grey Box interface to wrap these function calls into an IDAES property package. The WaterTAP team (https://github.com/watertap-org/watertap) has had a lot of success using this approach to integrating a tool called Reaktoro (https://reaktoro.org/ - this might even do some of what you want) using this approach. Whilst it is slower (about x5 in the current implementation I am told), it is not too slow to use. @bknueven might be able to provide you with some more details of this approach. As a whole, IDAES and its related projects are starting to look to this approach more to make use of the large number of excellent tools that already exist for calculating thermodynamic properties.
2. Pyomo also has a Trust Region tool which can be used to wrap external functions which lack partial derivatives. This tool works by using a small number of function calls to develop a low-fidelity surrogate of the function which is then used in the model and refined as it converges to the optimal point.
3. Implementing a new property package is definitely possible, and the good news is that you should only need to do the mixing rules part; you should be able to use the existing pure component Helmholtz functions inside this property package to do the pure component part. I am not sure what the mixing rules look like so I cannot say how easy this would be to write and solve, but it is definitely possible in theory.

[eric-hintz]

Thanks for the detailed reply!

I've heard about Reaktoro before, but I'm not too familiar, so I'll look into that more. Other than that, REFPROP does supply a decent amount of partial derivatives itself, so maybe it wouldn't be too bad to fill in the gaps and use the grey box interface that way.

I have read about the trust region solver, and it seems relatively easy to implements, but from what I understand (which may not be much), it seems to be more biased towards local, non-global solutions than other solvers. Because of that, it might be a lot more work on the back end than other solutions. Then again, it might be possible to bound the problem space sufficiently that this isn't a problem, so I might look into this at some point.

That's great to hear that I could build off the existing Helmholtz functions. Would the way to do this be making a custom package or using the modular package approach?
As far as the mixing rules, they are of similar or higher complexity to the pure fluid Helmholtz EOS. The CoolProp documentation on the topic is the best documentation of these equations that I have found (in case you aren't too familiar with CoolProp/REFPROP, they follow the exact same foundations, and are even developed by the same people, but REFPROP, as the commercial version, has more input/output flexibility and a more complete coefficient database). This is probably the most "ideal" solution, but rereading the CoolProp documentation now, I am worried this might require far more time than I can really dedicate at this point.

For now, I think I'll start with trying to use the grey box interface with REFPROP. I'll check out the WaterTAP project to see their implementation.

As an aside, if I end up with a solution of high enough quality, I would be happy to refactor the code into IDAES standards and contribute it to the project.

[andrewlee94]

Something I forgot to mention in my previous comments: the grey-box approach is limited in that it will only work with cyipopt as the solver (which IDAES is hoping to distribute soon as part of our binary packages). It is theoretically possible to link other solvers, but that would require effort at the Pyomo level.

For implementing the mixing rules, it would be a custom property package to start with (the modular framework could be extended to cover this, but that is a lot more work).

Finally, if you do have some success with this we would love to hear about it and accept a contribution.

[eric-hintz]

Cyipopt shouldn't be a problem I think, I probably would have been using IPOPT instead anyways.

Good to know, thank you again!

[adowling2]

If you are interesting in modeling refrigerants in IDAES, take a look at these papers from my research group:

https://www.sciencedirect.com/science/article/pii/S0378381223001139

https://pubs.acs.org/doi/10.1021/acs.iecr.2c01928

In these papers, we fit a cubic EoS within IDAES.