Fitting interaction parameters of binary systems
Sometimes Equations of State needs an extra help to better predict behaviors.
Let’s see how we can optimize binary interaction parameters (BIPs) of a Cubic Equation, using equilibrium data.
Data structure
For all the default approaches of parameter optimization the data points must be in the format of a pandas DataFrame, with the following columns:
kind: Which kind of point is defined. The options are:bubbleTbubblePdewTdewPPTliquid-liquidcritical
T: Temperature (in Kelvin).P: Pressure (in bar).x1: Mole fraction of component 1 (lightest component) in heavy phase.y1: Mole fraction of component 1 (lightest component) in light phase.
Optional column:
weight: A column that assign a weight factor to that data point when calculating the error. If not provided, all weights are assumed as 1.
Depending on the specified data point kind the experimental error will be calculated differently. Moreover, different information will be needed. “Not used” information could be set as NaN (missing data) in the dataframe
|
Description |
Not Used |
|---|---|---|
|
A saturation temperature calculation is performed using The error is calculated between the calculated saturation temperature and |
|
|
A saturation pressure calculation is performed using The error is calculated between the calculated saturation pressure and |
|
|
A saturation temperature calculation is performed using The error is calculated between the calculated saturation temperature and |
|
|
A saturation pressure calculation is performed using The error is calculated between the calculated saturation pressure and |
|
|
A flash calculation is performed at The error is calculated between the calculated light phase composition and and between the calculated heavy phase composition and You must provide |
— |
|
Same as |
— |
|
The critical line is calculated. Since is a critical point both The error is calculated between the calculated pressure, temperature and composition and |
|
As shown, depending on the kind of each data point, some information is mandatory and other is not used. This behaviour allows the user to adjust the same data points to different kinds and enhance the fitting results.
In the following cell we import the data points from a .csv file
[ ]:
import pandas as pd
import numpy as np
df = pd.read_csv('./data/CO2_C6.csv')
df
This data set consist of multiple bubble pressure data points at two different temperatures. Both liquid (x1) and vapor (y1) light component mole compositions are provided.
As explained before, depending on the kind of data point, some columns are not used. In this case, since we are going to use bubbleP kind, the y1 column is not used, even if it is still provided in the data set.
Since the weight column is not provided, all data points will be equally weighted when calculating the error.
Default model prediction
Before starting to fit BIPs, let’s see how the model predicts without BIPs.
First, let’s define the model
[ ]:
import yaeos
Tc = [304.1, 504.0] # Critical temperature in K
Pc = [73.75, 30.12] # Critical pressure in bar
w = [0.4, 0.299] # Acentric factor
# Peng-Robinson (1976) model without interaction parameters
model = yaeos.PengRobinson76(Tc, Pc, w)
Calculation of Pxy diagrams
Now that the model is defined, lets calculate Pxy diagrams for each temperature.
First, we find the unique values of temperature to iterate over them later.
[ ]:
Ts = df["T"].unique()
Ts
[ ]:
import matplotlib.pyplot as plt
gpec = yaeos.GPEC(model)
for T in Ts:
msk = df["T"] == T
pxy = gpec.plot_pxy(T, color="black")
plt.scatter(df[msk]["x1"], df[msk]["P"], label="Bubble Points")
plt.scatter(df[msk]["y1"], df[msk]["P"], label="Dew Points")
plt.xlim(0, 1)
plt.legend()
plt.show()
Here we can see that the PengRobinson EoS present a negative deviation with respect to the experimental points. Let’s se how we can improve these predictions.
Definition of the optimization problem
Now that we already have our data points defined. We need to define the optimization problem.
The yaeos.fitting package includes the BinaryFitter object, this object encapsulates all the logic that is needed to fit binary systems.
[ ]:
import yaeos
from yaeos.fitting import BinaryFitter
help(BinaryFitter)
Model setter function
Looking at the object’s documentation we can see that we need to define a function model_setter.
This is a user-defined function that sets and returns the Equation of State in function of the fitting parameters and extra required arguments.
The model_setter signature is as follows:
def model_setter(parameters_to_fit, *args) -> yaeos.core.ArModel:
"""
Function to set the Equation of State model.
Parameters
----------
parameters_to_fit : array_like
Array that contains the parameters to fit.
*args : tuple
Extra arguments required by the model.
Returns
-------
yaeos.core.ArModel
The Equation of State model with the specified parameters.
"""
# Code that sets an Equation of State model.
The extra arguments of the model_setter are passed to the BinaryFitter object when it is initialized. Let’s see an example to fit the \(k_{ij}\) values of the quadratic mixing rule (QMR) for the Peng-Robinson (1976) EoS.
[ ]:
def fit_kij(x, Tcs, Pcs, ws):
"""Fit kij function of QMR.
Parameters
----------
x : array_like
List of interaction parameters (kij) to fit.
Tcs : array_like
Critical temperatures of the components in K.
Pcs : array_like
Critical pressures of the components in bar.
ws : array_like
Acentric factors of the components.
"""
# Create the kij matrix from the argument
kij = x[0]
kij_matrix = [
[0, kij],
[kij, 0]
]
# lij matrix is null, we are only fitting kij.
lij = [
[0, 0],
[0, 0]
]
# Setting of the mixing rule
mixrule = yaeos.QMR(kij=kij_matrix, lij=lij)
# Create the Peng-Robinson (1976) model with the mixing rule
model = yaeos.PengRobinson76(Tcs, Pcs, ws, mixrule=mixrule)
# return the model
return model
With the model_setter function defined we can instantiate a BinaryFitter object. Let’s see how we can do that.
[ ]:
# The parameters of the substances
Tc = [304.1, 504.0] # Critical temperature in K
Pc = [73.75, 30.12] # Critical pressure in bar
w = [0.4, 0.299] # Acentric factor
# Instantiate the BinaryFitter
problem = BinaryFitter(
model_setter=fit_kij, # model_setter function
model_setter_args=(Tc, Pc, w), # Extra arguments for the model_setter
data=df, # Data to fit
verbose=True # Print information during the fitting process
)
# Initial guess for the fitting parameters (kij)
x0 = [0.0]
# Call the fit method with the initial guess and fit the model
problem.fit(x0, bounds=None)
We can check the arguments of the fit method.
The first argument x0 is the initial guess for the fitting parameters.
The second argument, bounds, could be used to set the optimization bounds of the interaction parameters. Is defined as a tuple of two arrays, the first one contains the lower bounds and the second one contains the upper bounds for each parameter. To not specify bounds, set it to None.
The third argument, method, is used to select the optimization method. The default method is "Nelder-Mead". Check the scipy.optimize.minimize documentation for more information about the available methods.
[ ]:
help(problem.fit)
Optimization result
Now that we have fitted the parameter, let’s check the solution. For this, we can access the solution attribute of problem.
As shown we can see that the optimization terminated successfully with a kij value of around 0.113
[ ]:
problem.solution
Obtain the fitted model
Now that the problem has been optimized, let’s redefine our model with the solution.
For this, we use the get_model method inside the BinaryFitter object. This just uses the function that we have provided originally (fit_kij) to return the model with the specified parameters.
[ ]:
# Obtain the fitted model specifying the solution of the fitting and the parameters of the substances
model = problem.get_model(problem.solution.x, Tc, Pc, w)
Make predictions
Let’s repeat the calculation of Pxy diagrams that we have done earlier
[ ]:
import matplotlib.pyplot as plt
gpec = yaeos.GPEC(model)
for T in Ts:
msk = df["T"] == T
gpec.plot_pxy(T)
plt.scatter(df[msk]["x1"], df[msk]["P"], label="Bubble Points")
plt.scatter(df[msk]["y1"], df[msk]["P"], label="Dew Points")
plt.xlim(0, 1)
plt.legend()
plt.show()
Another example
Let’s check another example of fitting BIPs. First, we need to read the data.
[ ]:
import pandas as pd
df = pd.read_csv("data/C1_benzene.csv")
# Temperature of the experimental data
T = 373
df
Here we have a data set of bubble pressures and dew pressures all at 373 K. Moreover, the critical point is also provided at 373 K.
In this mixture of different data kind, the not used information is set as NaN (missing data) in the dataframe.
First, we are going to instantiate a Peng-Robinson EoS model without BIPs and see how it predicts the data.
[ ]:
# Substance properties
Tc = [190.564, 562.05] # Critical temperatures in K
Pc = [45.99, 48.95] # Critical pressures in bar
w = [0.0115478, 0.2103] # Acentric factors
Zc = [0.286, 0.268] # Critical compressibility factors
model = yaeos.RKPR(Tc, Pc, w, Zc)
Now let’s plot the experimental data and the model predictions.
[ ]:
import matplotlib.pyplot as plt
gpec = yaeos.GPEC(model)
pxy = gpec.plot_pxy(T)
msk_bub = df["kind"] == "bubbleP"
msk_dew = df["kind"] == "dewP"
msk_critical = df["kind"] == "critical"
plt.scatter(df[msk_bub]["x1"], df[msk_bub]["P"], label="Bubble Points")
plt.scatter(df[msk_dew]["y1"], df[msk_dew]["P"], label="Dew Points")
plt.scatter(df[msk_critical]["x1"], df[msk_critical]["P"], label="Critical Point")
[ ]:
def fit_kij(x, Tcs, Pcs, ws, Zc):
"""Fit kij function of QMR.
Parameters
----------
x : array_like
List of interaction parameters (kij) to fit.
Tcs : array_like
Critical temperatures of the components in K.
Pcs : array_like
Critical pressures of the components in bar.
ws : array_like
Acentric factors of the components.
"""
# Create the kij matrix from the argument
kij = x[0]
kij_matrix = [
[0, kij],
[kij, 0]
]
# lij matrix is null, we are only fitting kij.
lij = [
[0, 0],
[0, 0]
]
# Setting of the mixing rule
mixrule = yaeos.QMR(kij=kij_matrix, lij=lij)
# Create the Peng-Robinson (1976) model with the mixing rule
model = yaeos.RKPR(Tcs, Pcs, ws, Zc, mixrule=mixrule)
# return the model
return model
problem = BinaryFitter(
model_setter=fit_kij, # model_setter function
model_setter_args=(Tc, Pc, w, Zc), # Extra arguments for the model_setter
data=df, # Data to fit
verbose=True # Print information during the fitting process
)
[ ]:
x0 = np.array([0.0])
problem.fit(x0=x0, bounds=None)
[ ]:
model = problem.get_model(problem.solution.x, Tc, Pc, w, Zc)
[ ]:
import matplotlib.pyplot as plt
gpec = yaeos.GPEC(model)
pxy = gpec.plot_pxy(T)
msk_bub = df["kind"] == "bubbleP"
msk_dew = df["kind"] == "dewP"
msk_critical = df["kind"] == "critical"
plt.scatter(df[msk_bub]["x1"], df[msk_bub]["P"], label="Bubble Points")
plt.scatter(df[msk_dew]["y1"], df[msk_dew]["P"], label="Dew Points")
plt.scatter(df[msk_critical]["x1"], df[msk_critical]["P"], label="Critical Point")
plt.xlim(0, 1)
plt.legend()
plt.show()
[ ]:
def fit_nrtl_hv(x, Tc, Pc, ws, Zc):
g12, g21 = x[:2]
alpha = x[2]
gji = [[0, g12], [g21, 0]]
alpha = [[0, alpha], [alpha, 0]]
use_kij = [[False, False], [False, False]]
kij = [[0, 0], [0, 0]]
mixrule = yaeos.HVNRTL(alpha=alpha, gji=gji, use_kij=use_kij, kij=kij)
model = yaeos.RKPR(Tc, Pc, ws, Zc, mixrule=mixrule)
return model
[ ]:
Tc = [190.564, 562.05] # Critical temperatures in K
Pc = [45.99, 48.95] # Critical pressures in bar
w = [0.0115478, 0.2103] # Acentric factors
Zc = [0.286, 0.268] # Critical compressibility factors
problem = BinaryFitter(
model_setter=fit_nrtl_hv, # model_setter function
model_setter_args=(Tc, Pc, w, Zc), # Extra arguments for the model_setter
data=df, # Data to fit
verbose=True # Print information during the fitting process
)
[ ]:
x0 = [31.81663237, 4.10864439, -1.3615617]
x0 = np.zeros(3)
problem.fit(x0=x0, bounds=None, method="Powell")
[ ]:
problem.solution
[ ]:
model = problem.get_model(problem.solution.x, Tc, Pc, w, Zc)
[ ]:
import matplotlib.pyplot as plt
gpec = yaeos.GPEC(model)
pxy = gpec.plot_pxy(T)
msk_bub = df["kind"] == "bubbleP"
msk_dew = df["kind"] == "dewP"
msk_critical = df["kind"] == "critical"
plt.scatter(df[msk_bub]["x1"], df[msk_bub]["P"], label="Bubble Points")
plt.scatter(df[msk_dew]["y1"], df[msk_dew]["P"], label="Dew Points")
plt.scatter(df[msk_critical]["x1"], df[msk_critical]["P"], label="Critical Point")
plt.xlim(0, 1)
plt.legend()
plt.show()