Parameter estimation and its tutorial

[jamiecripwell]

Hi devs

So I am wanting to start using Clapeyron as a full suite - including parameterisation - to encourage my students to do the same outside of our in-house software. As it stands though, I am struggling to reproduce even simple literature SAFT parameter sets - even using the .ipynb tutorial file for methane. There are a couple of things that stand out as potential sources here:

• the reference data and units: to my understandingm everything in Clapeyron is SI based so (e.g.) vapour pressure and liquid density should be in $Pa$ and $\frac{mol}{m^3}$ correct? In which case, is the methane saturation data in the example just mock data? Or in weird units?
• unfamiliarity with the optimisation routines of Metaheuristics.jl: thus far I have simply tried increasing/removing the iterations in the method as described in the tutorial markdown, but this makes little difference and the final optimal parameters are bad. Could this be a local minimum problem? I can certainly read further into the different algorithms, but given my uncertainty over the data above, I don't want to go down the wrong rabbit hole (and a tutorial should give correct answers to simple problems right?). I did try replace it with my own data in SI units, but the parameters are still bad.
• I could also spend time checking the estimation.jl file, but I would just like to make sure I'm getting the basic functionality right first.

Any help would be much appreciated!
Jamie

[pw0908]

Hi Jamie,

Happy to hear you're interested in using Clapeyron! To answer each of your questions:

• It's been a few years since I made these notebooks (I'm in the process of updating a lot of them). Nevertheless, I believe this data came from GERG-2008 (I don't even think the data is for methane). Essentially, it is mock data, but should still be close to experimental values. The units are indeed in SI (Pa and mol/m3).
• I went ahead and re-ran the notebooks and, for the data provided, the fit is excellent:
  
  
  With parameters:

• epsilon: 113.85112463206652
• sigma: 3.2873052264106515
• segment: 1.068112572714932

I'm not sure what you mean by the fitting is 'bad'. Do you mean in comparison to the base PC-SAFT parameters for methane? In which case, that's to be expected since the data used is not the same and, in this case, we're allowing the segment number to vary. If you're comparing to experimental data, then the same logic can be applied, with the added caveat that the fitting didn't use experimental data close to the critical point. And this is before one considers the degeneracies present when doing fitting of PC-SAFT parameters.

If you want more information regarding Metaheuristics.jl, they have an excellent set of docs: https://jmejia8.github.io/Metaheuristics.jl/stable/. In general, they provide global optimization algorithms using heuristic methods.

• It seems like your understanding of our estimation framework is fine, as long as we clear up the comments above. Are you interested in the particular details regarding how the parameters themselves are treated?

Best regards,

Pierre

[jamiecripwell]

Hi Pierre

Thanks so much for replying so quickly - one of the reasons I want my students to port to Clapeyron!

Your reply helped me find the two issues:

• I obviously got the same parameters, and stuck them into our in-house software for comparison, which produced the figure below (my meaning of bad):
  

• The data you provide (and therefore fit to, and compare the prediction against above) is not representative of methane at all - I include a comparison of that data compared to the (agreeing) data from NIST and DIPPR below. So I think an update of that data file would help the confusion.
  

• In trying to reproduce parameters, I was using the density as mass density in SI, which caused the other issue - this is standard practice for us as we use the DIPPR database, but I don't know if that is true for other groups. Perhaps a note in the tutorial notebook just reminding of this point?

Taken together (molar density and updated data), I can reproduce literature parameters so that is sorted, thanks!

While on the subject of parameterisation - would it require significant tinkering to the estimation.jl file to get the fitting metrics as part of the regression output (as is "usually" the case with the Metaheuristics package)?

Thanks!
Jamie

[pw0908]

Hi Jaime,

I now see what you mean. Yes, the data is definitely not that for methane (honestly, it's been so long since I wrote that code, I don't remember what molecule it was supposed to correspond to). I've now changed the data to correspond to methane and I am able to reproduce the PC-SAFT parameters.

As for the units, we've documented most of the units in our docs. You'll see that there is a mass_density function that obtains the density in mass units. You might also want to check out Unitful.jl which allows for easy unit conversions (we've also added an extension to make Clapeyron compatible with it).

The way the parameterisation is presented in the docs, we use the Clapeyron-Metaheuristics extension which allows users to just specify optimize(objective, estimator, method). This is more of a convenience function and, if you wanted to just use Metaheuristics directly and obtain its output, you can do so:

nparams = length(estimator.toestimate.upper)
upper = [estimator.toestimate.upper[i][1] for i ∈ 1:nparams]
lower = [estimator.toestimate.lower[i][1] for i ∈ 1:nparams]
bounds  = [lower upper]'
if verbose
    results = Metaheuristics.optimize(f,bounds,method;logger=logger)
else
    results = Metaheuristics.optimize(f,bounds,method;logger=(status) -> nothing 
end
  
params = minimizer(results)
model = return_model(estimator,estimator.model,params)

This is pretty much was the extension does automatically. See here.

Best regards,

Pierre

[jamiecripwell]

Thanks for your feedback Pierre!

I am getting the hang of the different aspects, but have two further questions I'd appreciated your help on with the parameterisation aspect:

1. There is an inconsistency with the "SI" approach for the association volume between the SAFT variants, specifically for SAFT-VR Mie: here, values in the database have units of cubic meters (e-29 takes it to angstrom), but the other SAFT variants seem to be in cubic angstrom (no huge exponential terms), although this is not clear from the documentation (only parameter with no explicit units). More generally, the approach is inconsistent with the presentation of SAFT-VR Mie's use of the range of association ($r_c^{ab}$) and the volume is thus temperature dependent - see equation 84 in Lafitte's paper. Importantly for prediction, the respective models do work with these different unit systems, although I have yet to find the point of divergence in the code - I am still tracing the association term's working more generally. This makes using SAFT-VR Mie and its variants a little challenging as a pseudo-bondvol parameter has to be used (rather than published $r_c^{ab}$) and indeed parameterisation using this approach just becomes inconsistent with the literature. This is something that I can always amend and contribute to the project when I am more thoroughly familiar, but I just wanted to bring it to your attention.
2. More for dev consideration - the parameterisation approach caters well for fitting SAFT parameters to either pure component data or mixture data, as illustrated in the example notebook, but the inherent EOSModel structure seems to make it impossible (probably just difficult) to consider both. This is because the object is inherently defined by its number of components and the pure component algorithms will not run on an object with two components - one can easily verify this by combining the bubble_point and saturation_p_rhol functions in the example notebook for a SAFT model with two components. I attempted to incorporate the following within the saturation_p_rhol function to "rebuild" the single-component object to be used in the saturation_pressure function. It successfully does that, but the very poor results of sending this through the optimisation suggest I am missing some important "carry-over" of the model parameters during the optimisation:

function saturation_p_rhol(model_multi::EoSModel,T) params_vector = fieldnames(typeof(model_multi.params))

user_locations_string = ";"
for i ∈ 1:length(params_vector)
        param_name = params_vector[i]
        if getfield(model_multi.params, params_vector[i]).values isa AbstractArray{Float64}
            if param_name == :sigma
                param_value = getfield(model_multi.params, param_name).values[1].*1e10
                user_locations_string *= "$param_name=$param_value,"
            else
                param_value = getfield(model_multi.params, param_name).values[1]
                user_locations_string *= "$param_name=$param_value,"
            end
        else
            param_value = getfield(model_multi.params, param_name).values.values[1]
            comp_name = model_multi.components[1]
            user_locations_string *= """$param_name=Dict((("$comp_name","e"),("$comp_name","H")) => $param_value),"""
        end
end
user_locations_string = user_locations_string[1:end-1]

 
model_type_name = nameof(typeof(model_multi));
comps = [model_multi.components[1]]

 expression = "$model_type_name($comps; userlocations=($user_locations_string))"

 model = eval(Meta.parse(expression))

sat = saturation_pressure(model,T)
return sat[1], 1/sat[2]`

Am I missing a simpler solution? Is there a way to instead feed a multicomponent EOSmodel object to the pure component algorithms to allow me to include both pure component properties and mixture data in the regression?

Thank you!

[jamiecripwell]

Ok, I seem to have figured out part 2 to some extent - simply leveraging the bubble_pressure algorithm and specifying the "mixture" as pure component 1.

[pw0908]

Hi Jaime,

Addressing each of your questions:

1. I think there's a bit of a mix-up here. There are inherently three versions of SAFT-VR Mie, each using what is referred to as a different "association kernel". In the original Lafitte 2013 paper, they used what is now called the hard-sphere kernel where the association length-scale parameter is $r_c^{ab}$ in Angstrom. However, this is the last time the Imperial group has used the hard-sphere kernel. Since then, they have implemented two more association kernels in Dufal's 2015 paper: the LJ and Mie kernels. These were developed to better line-up with the Mie reference system. These kernels both use the bonding volume, $K_{ij,ab}$, parameter (in Angstrom^3, hence why there is the 1e-30 multiplier). The 'default' association kernel now used in all SAFT-VR Mie publications since (and its successor, SAFT-gamma Mie) is the LJ kernel. To be consistent with these publications, we use the LJ kernel by default. That being said, there is one branch of SAFT-VR Mie users (eSAFT-VR Mie) which use the Mie kernel; this version is provided as SAFTVRMie15 but is only really used in eSAFTVRMie.

Generally, you'll find that the original developers of SAFT-VR Mie have made many changes through the years and haven't published all of them (e.g. their use of Gauss-Laguerre quadrature instead of Gauss-Legendre as they say in the 2013 paper). The only reason we are aware of these changes is because we've had many correspondences with them to make sure our implementation matches their's. You'll notice in our references for the SAFT-VR Mie model we cite both Lafitte's 2013 paper and Dufal's 2015 paper to highlight that our implementation follows from both.

2. Clearly this means we really need to develop more-rigorous examples... the code is currently implemented to make it easy to handle parameter estimation involving any number of components, even when the data used only involves pures or binaries. This will be a bit complicated to explain in words so I've made an example script based on the data you've sent, fitting the hexane-methanol binary mixture, as well as the saturation properties of hexane, simultaneously. For your own reference, you'll see in our example notebook for parameter estimation that, in the very last example (Group-contribution parameter estimation), it is possible to specify in the header of the CSV spreadsheet that the experimental data only involves a subset of the species in the model:

Clapeyron Estimator,,
[method = saturation_p_rhol, species = hexane],,
T,out_pre,out_rho

What this is does is tell Clapeyron, when evaluating saturation_p_rhol for this data, to reduce the model so that it only includes hexane. You can do this for any number of species (the separator is a blank space). I've personally used this to fit models involving over a dozen species at once, even though most of the data used is pure, binary or ternary.

I've attached the example here. You'll need to update your version of Clapeyron to the master branch as there is a typo within the parameter estimation framework.
example.zip

[jamiecripwell]

Thanks Pierre, that helps - I understand your point on the association kernels, and I can see how the divergent paths would make it difficult to cater to all approaches, I will just note for future for SAFT VR Mie.

Your example was helpful, but in many unexpected ways! The most notable being that you removed the pure component measurements from the binary mixture data set, and this made all the difference in the parameterisation. For some reason, if the pure component points are there, the optimisation seems to diverge (perhaps its the "0" composition?) - the logger helped show this as it stops after a single iteration with a very large OF value (this was the richer information I mentioned in my earlier post). So again, something I will note for my own future use, but I don't know if this is intentional? Some data sets contain the anchoring pure component points and others not, so this seemed to be poor luck in my first attempt.

Besides that, a couple of questions for my own clarity on the data in the csv files:

• Do the headers matter? In the saturation properties, pressure has the header out_pre while for the bubble algorithm it is out_p? And if it does matter, what are the headers for, say, speed of sound? Is there a "key" somewhere that has these mappings?
• Alternatively, do the column numbers in the csv file just need to match the output indexes for the OF? So for example, in the saturation_p_rhol function, pressure is the first output and density the second (return sat[1] 1/sat[2]) - and in the csv file, pressure is the first output column and density the second.
• You added the "species" keyword in the data csv files - are there other useful keywords for parameterisation? And are they documented somewhere?
• If you have multicomponent data for the regression (i.e. $x_2 \neq 1-x_1$), how do you set up the data file for these compositions? I suppose this is tied to my first question about whether the headers are important - if it just needs to match the OF function outputs, then I guess you would need to put $x_1$ as first output, first data column, $x_2$ as second output, second data column, etc. But confirmation on this would be great.

And thanks again for your responsiveness - I think it makes all the difference with adoption of open source packages.

[pw0908]

I'm glad this has helped clear some things up! You are not the only one who's contact us about inconsistencies in SAFT-VR Mie (but when there are 5 different version published in the literature, there's only so much we can do...).

Regarding the removing the pure-component information from the mixture data, as you may have guessed, the fact that the vapour composition in the gas phase will also be zero leads to a 0/0 situation which causes the objective function to diverge. One solution would be to modify the objective function form (which you can do in Clapeyron by changing the objective_form: https://clapeyronthermo.github.io/Clapeyron.jl/stable/api/estimation/#Clapeyron.Estimation). In addition, I would make the argument that including this pure-component data, especially we one already has a separate spreadsheet with the pure component data, biases the fitting somewhat.

For the remainder of your questions:

1. The headers of the columns do not matter much beyond including the out_ prefix. This helps Clapeyron distinguish between an input (no prefix) and an output (out_ prefix). You could have a column header be out_xyz and, Clapeyron will simply match the property based on the index only.
2. There aren't any other keywords other than species and method... we initially considered adding support for including property errors in case users wanted to provide uncertainty tracking (which you can technically do right now), but we never fully developed it. You'll find more information within the EstimationData struct: https://clapeyronthermo.github.io/Clapeyron.jl/stable/api/estimation/#Clapeyron.EstimationData
3. For multicomponent systems, you are correct, each composition would need to be a separate output. Now, if you wanted to do something fancy and define the error based on the vector distance between the experimental and predicted points in composition space, you'd need to define the objective function manually (i.e. not used the automatically output one from Estimation). This is doable and you might want to take a look at our objective function code to get an idea of the structure:

   Clapeyron.jl/src/estimation/estimation.jl

   Line 326 in 5ecca86

   function objective_function(estimation::Estimation,guesses)