This implementation can be rewritten using linear algebra matrix notation, this can simplify analytical derivations, improve readability (although for the trained mind), maintainability, reduce human error in the coding (less typing therefore less typos), and improve execution speed (if we make use of libraries that can automatically convert matrix notation to optimal implementation of the underlying loops - e.g. NumPy). For details on the derivations of the matrix notation for the NRTL model, as well as for UNIQUAC, UNIFAC and COSMO see the lectures notes of Abreu, C. R. A. in our library here.
- - - - - - - - - - - - - - Renon & Prausnitz - - - - - - - - - - - - - - | - - - - - - - - - - - - - - Abreu - - - - - - - - - - - - - - |
---|---|
$ \frac{g^E}{RT}=\sum_{i=1}^n \left[ x_i\frac{\sum_{j=1}^n \tau_{j,i} G_{j,i} x_{j}}{\sum_{k=1}^n G_{k,i} x_k} \right] $ | $ \frac{g^E}{RT}=\underline{x}^T\underline{\underline{E}} \underline{x}$ |
Where
- - - - - - - - - - - - - - Renon & Prausnitz - - - - - - - - - - - - - - | - - - - - - - - - - - - - - Abreu - - - - - - - - - - - - - - |
---|---|
$\tau_{i,j}= \frac{g_{i,j}-g_{j,j}}{RT}=\frac{A_{i,j}}{T}$ | $\underline{\underline{\tau}}=-T^{-1}\underline{\underline{A}}$ |
$G_{i,j}=\mathrm{exp}(-\alpha_{i,j} \tau_{i,j})$ | $\underline{\underline{G}}= \mathrm{exp}(\underline{\underline{\alpha}} \circ \underline{\underline{\tau}})$ |
--- | $\underline{\underline{\Lambda}}=(\underline{\underline{\tau}} \circ \underline{\underline{G}})$ |
--- | $\underline{\underline{E}}=\underline{\underline{\Lambda}}\mathscr{D}^{-1}(\underline{\underline{G}}^T\underline{x})$ |
Therefore
- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - Renon & Prausnitz - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - | - - - - - - - - - - - Abreu - - - - - - - - - - - |
---|---|
$ln(\gamma_i)= \frac{\sum_{j=1}^n\left[\tau_{j,i} G_{j,i} x_{j}\right]}{\sum_{k=1}^n\left[G_{k,i}x_{k}\right]} + \sum_{j=1}^n\left[ \left(\frac{\ G_{i,j} x_{j}}{\sum_{k=1}^n\left[G_{k,j}x_{k}\right]}\right) \left(\tau_{i,j}-\frac{\sum_{j=1}^n\left[\tau_{i,j} G_{i,j} x_{i}\right]}{\sum_{k=1}^n\left[G_{k,j}x_{k}\right]} \right) \right] $ | $ ln(\underline{\gamma})=\left(\underline{\underline{E}}^S-\underline{\underline{L}}\mathscr{D}\underline{x}\underline{\underline{E}}^T \right)\underline{x}$ |
Where
Replace summations of matrices in one of its dimensions
$C_{i,j} = \sum_k{\left[A_{i,k} \times B_{k,j}\right]}$
$d_{i} = \sum_k{\left[A_{i,k} \times b_{k}\right]}$
$e_{j} = \sum_k{\left[a_{k} \times B_{k,j}\right]}$
$\underline{\underline{C}}=\underline{\underline{A}}\underline{\underline{B}}$
$\underline{d}=\underline{\underline{A}}\underline{b}$
$\underline{e}=(\underline{a^T}\underline{\underline{B}})^T$
C = A @ B
d = A @ b
e = (a.T @ B).T
Condensate operations that apply analogously in all elements of one or more matrices
$C_{i,j} = A_{i,j} \times B_{i,j}$
$\underline{\underline{C}}=\underline{\underline{A}} \circ \underline{\underline{B}}$
C = A * B
Diagonalization is used to represent element-wise multiplication between: two single line matrixes, two single columns, one single line matrix and each line of a (nl,nc) matrix, or one single column matrix and each column of a (nl,nc) matrix. This is useful in analytical differentiation presented in the lecture notes of Abreu C. R. A..
On the other hand, this notation may be dropped later, at implementation time, for a numerical solution, depending on the programming enviroment, see below the usage and correspondence of notations:
$C_{i,j} = A_{i,j} \times b_{j}$
$D_{i,j} = b_{i} \times A_{i,j}$
$e_{i} = a_{i} \times b_{i}$
$f_{j} = b_{j} \times a_{j} $
$\underline{\underline{C}}=\underline{\underline{A}} \left( \mathscr D \underline{b} \right )$
$\underline{\underline{D}}= \left( \mathscr D \underline{b} \right ) \underline{\underline{A}}$
$\underline{e}= \left( \mathscr D \underline{a} \right ) \underline{b}$
$\underline{f}=\left( \underline{b}^{T} \left( \mathscr D \underline{a} \right ) \right)^T$
C = A * b.T
D = b * A
e = a * b
f = (b.T * a.T).T
where we chose to use:
- upper case letters to represent full (nl,nc) matrixes
- lower case letters to represent single column matrixes)
while from Python/NumPy syntax:
M.T
is the numpy native syntax for transposition of a matrixM
A \* B
is the numpy native syntax for the element-wise multiplication between matrixesA
andB
Note that in matrix algebra in python, single columns and single line matrixes should be represented by 2d arrays having shape (nl,1) and (1,nc), respectively. This is not quite the same thing as 1d arrays. Understand the differences seeing usage examples here.
where:
- nl stands for number of lines in a single column (nl,1) matrix or in a full (nl,nc) matrix
- nc stands for number of columns in a single line (1,nc) matrix or in a full (nl,nc) matrix
Renon et al. (1969) fitted 1 $\alpha$ valid for all 3 binary interactions $\{(1,2),(1,3),(2,3)\}$ parameter and 6 $A_{i,j}$ parameters, two for each binary interaction filling a non symmetric $A$ matrix.
# Ethyl acetate (1) + water (2) + ethanol (3)
alpha12 = 0.4
alpha23 = 0.3
alpha13 = 0.3
# 6 binary Aij parameters
Dg12 = 1335 * 4.184 #J/K
Dg21 = 2510 * 4.184 #J/K
Dg23 = 976 * 4.184 #J/K
Dg32 = 88 * 4.184 #J/K
Dg13 = 301 * 4.184 #J/K
Dg31 = 322 * 4.184 #J/K
we will assemble the parameters in a matrix structure so that we can access each parameter by its index, as in
A[0,0]
and A[0,1]
rather than as A11
and A12
, so we can loop trough all of them using an iterator, see below:
import numpy as np
from scipy.constants import R
print(R)
8.3144598
#assemble matrix with regressed parameters Dg_i,j, according to the model all diagonal terms are zero
Dg = np.array([[0, Dg12, Dg13],
[Dg21, 0, Dg23],
[Dg31, Dg32, 0]])
A = Dg/R
#assemble symmetric matrix alpha
alpha = np.array([[0, alpha12, alpha13],
[alpha12, 0, alpha23],
[alpha13, alpha23, 0]])
We previously had:
def Gamma(T,x,alpha,A):
tau=np.zeros([3,3])
for j in range(3):
for i in range(3):
tau[j,i]=A[j,i]/T
G=np.zeros([3,3])
for j in range(3):
for i in range(3):
G[j,i]=np.exp((-alpha[j,i]*tau[j,i]))
Gamma=np.zeros([3])
for i in range(3):
Sj1=0
Sj2=0
Sj3=0
for j in range(3):
Sj1 += tau[j,i]*G[j,i]*x[j]
Sj2 += G[j,i]*x[j]
Sk1=0
Sk2=0
Sk3=0
for k in range(3):
Sk1 += G[k,j]*x[k]
Sk2 += x[k]*tau[k,j]*G[k,j]
Sk3 += G[k,j]*x[k]
Sj3 += ((x[j]*G[i,j])/(Sk1))*(tau[i,j]-(Sk2)/(Sk3))
Gamma[i]=np.exp(Sj1/Sj2 + Sj3)
return Gamma
#test it to see if results match
#trial temperature and composition:
T = 293.15 #K
x=np.array([.1,.3,.6]) #normalized
ans=Gamma(T,x,alpha,A)
print(ans) #ttest using those trial input
[ 1.4967744 1.28850578 1.01628367]
def Gamma_linalg(T,c_x,q_alpha, q_A): # here we chose to use the starting letters s, c, l, and Q to identify scalar variables, single column matrixes, single line matrixes and square matrixes to the reader
# e_T should be an scalar value for temperature
# c_x should be a single-column matrix(2d array) representing composition
# q_alpha should be two matrix(2d array) representing component dependent parameters inferred from experimental data
# q_tau should be two matrix(2d array) representing component dependent parameters inferred from experimental data
q_tau = q_A/T #element-wise division by scalar
q_at = q_alpha*q_tau #M2d * N2d yields element-wise multiplication
q_minusat = -q_at #element wise signal change
q_G = np.exp(q_minusat) #element wise exponentiation
q_Lambda = (q_tau*q_G) #element-wise multiplication
q_GT = q_G.T #M.T yields the transpose matrix of M;
c_den = q_GT @ c_x #M @ N yields the matrix multiplication between M and N
c_invden = 1/c_den #scalar 1 is broadcast for element-wise division
l_invden = c_invden.T #transposition of a single column matrix yields a single line matrix
q_E = q_Lambda * l_invden #element wise multiplication between (nl,nc) matrix M with (1,nc) matrix l broadcasts the element-wise multiplication of each (1,nc) submatrix of M with the unique (1,nc) matrix l
q_L = q_G * l_invden #broadcasting of element-wise multiplication
l_x = c_x.T #transposition of a single column matrix yields a single line matrix
q_Lx = q_L * l_x #broadcasting of element-wise multiplication
q_ET = q_E.T #transposition of square matrix
q_LxET = q_Lx @ q_ET #matrix multiplication
q_ES = q_E+q_ET #element-wise sum
q_ESminusLxET = q_ES-q_LxET #element-wise subtraction
q_ESminusLxETx = q_ESminusLxET @ c_x #matrix multiplication
gamma = np.exp(q_ESminusLxETx) #element-wise exponentiation
return gamma
#a test case for the function
#where x was the composition represented in a 1d array
#and now line and x_as_column is a single lne and a single column matrix, respectively to represent composition
#We build it using the array function to wrap the 1d array in another 1d aray, hence a 2d array
x_as_line = np.array([x])
#We transpose x_as_line to creata x_as_column, which is the shape expected by the linalgGamma function
x_as_column = np.array([x]).T #we wrap x with an extra braket so it is now a 2d array, a matrix, as we did not add any extra lines it is a single-line matrix, we tranpose to generate a single-column matrix (1d arrays cannot be transposed, there is no second dimension)
#print the output to see if errors occur and if values are coherent(between zero and infinity, tending to 1 for ideal solutions)
print(Gamma_linalg(T,x_as_column,alpha,A)) #test using those trial input
[[ 1.4967744 ] [ 1.28850578] [ 1.01628367]]
def Gamma_linalg_tiny(T,c_x,q_alpha, q_A):
#note that we used many lines for didatics
#we can do it in few lines:
#note that some expression occur more than once below
#so it may be useful define it as a intermediary recurrent term here
#and calculate it once to use it then several times
q_tau = q_A/T
q_G = np.exp(-(q_alpha*q_tau))
l_D = ((1/((q_G.T) @ c_x)).T)
q_E = (q_tau*q_G) * l_D
gamma = np.exp(((q_E+(q_E.T))-(((q_G * l_D) * (c_x.T)) @ (q_E.T))) @ c_x)
return gamma
#test it to see that the results are the same
print(Gamma_linalg_tiny(T,x_as_column,alpha,A)) #test using those trial input
[[ 1.4967744 ] [ 1.28850578] [ 1.01628367]]
Ipython provides a profiling tool, with %timeit you can evaluate the time for execution of a line of program (with all called dependencies). In the following cells, we use it to evaluate our function in version 1 with explicit for loops and in version 2a and 2b with linear algebra matrix operations
%timeit Gamma(T,x,alpha,A) #ttest using those trial input #My result was 90 micro seconds per loop
10000 loops, best of 3: 102 µs per loop
%timeit Gamma_linalg(T,x_as_column,alpha,A) #ttest using those trial input #My result was 25 micro seconds per loop
The slowest run took 198.08 times longer than the fastest. This could mean that an intermediate result is being cached. 10000 loops, best of 3: 28.1 µs per loop
%timeit Gamma_linalg_tiny(T,x_as_column,alpha,A) #ttest using those trial input #My result was 25 micro seconds per loop
The slowest run took 174.12 times longer than the fastest. This could mean that an intermediate result is being cached. 10000 loops, best of 3: 27.2 µs per loop
here more honest trials with random input
%%timeit
#approximately time the random number generation to subtract later
# ~21 micro seconds per loop here
N=3
x=np.random.rand(N,1)
x=x/sum(x)
alpha=np.random.rand(N,N)
A=np.random.rand(N,N)
T=(np.random.rand(1,1)+.5)*273
The slowest run took 431.19 times longer than the fastest. This could mean that an intermediate result is being cached. 10000 loops, best of 3: 21.7 µs per loop
%%timeit
# ~440 micro seconds per loop here (420 subtracting the random number generation)
N=3
x=np.random.rand(N,1)
x=x/sum(x)
alpha=np.random.rand(N,N)
A=np.random.rand(N,N)
T=(np.random.rand(1,1)+.5)*273
_=Gamma(
T,
x,
alpha,
A)
The slowest run took 28.12 times longer than the fastest. This could mean that an intermediate result is being cached. 1000 loops, best of 3: 426 µs per loop
%%timeit
# ~56 micro seconds per loop here (36 subtracting the random number generation)
N=3
x=np.random.rand(N,1)
x=x/sum(x)
alpha=np.random.rand(N,N)
A=np.random.rand(N,N)
T=(np.random.rand(1,1)+.5)*273
_=Gamma_linalg(
T,
x,
alpha,
A)
The slowest run took 189.96 times longer than the fastest. This could mean that an intermediate result is being cached. 10000 loops, best of 3: 55.7 µs per loop
%%timeit
# ~52 micro seconds per loop here (32 subtracting the random number generation)
N=3
x=np.random.rand(N,1)
x=x/sum(x)
alpha=np.random.rand(N,N)
A=np.random.rand(N,N)
T=(np.random.rand(1,1)+.5)*273
_=Gamma_linalg_tiny(
T,
x,
alpha,
A)
The slowest run took 195.10 times longer than the fastest. This could mean that an intermediate result is being cached. 10000 loops, best of 3: 56.2 µs per loop
Also, with minor effort, we can use Numba to further accelerate the code.
See below how it is able to accelerate the functions Gamma
and Gamma_linalg_tiny
and compare.
#These two lines is all that it takes to accelerate this function
from numba import jit
@jit
#now repeat the function with a different bname so we can compare
def Gamma_numba(T,x,alpha,A):
tau=np.zeros([3,3])
for j in range(3):
for i in range(3):
tau[j,i]=A[j,i]/T
G=np.zeros([3,3])
for j in range(3):
for i in range(3):
G[j,i]=np.exp((-alpha[j,i]*tau[j,i]))
Gamma=np.zeros([3])
for i in range(3):
Sj1=0
Sj2=0
Sj3=0
for j in range(3):
Sj1 += tau[j,i]*G[j,i]*x[j]
Sj2 += G[j,i]*x[j]
Sk1=0
Sk2=0
Sk3=0
for k in range(3):
Sk1 += G[k,j]*x[k]
Sk2 += x[k]*tau[k,j]*G[k,j]
Sk3 += G[k,j]*x[k]
Sj3 += ((x[j]*G[i,j])/(Sk1))*(tau[i,j]-(Sk2)/(Sk3))
Gamma[i]=np.exp(Sj1/Sj2 + Sj3)
return Gamma
from numba import jit
@jit
def lngammaNRTL(T,c_x,q_alpha, q_A):
q_tau = q_A/T
q_G = np.exp(-(q_alpha*q_tau))
l_D = ((1/((q_G.T)@
c_x)).T)
q_E = (q_tau*q_G)*l_D
return (((q_E+(q_E.T))-(((q_G*l_D)*(c_x.T))@
(q_E.T)))@
c_x)
#These two lines are all that it takes to accelerate this function
from numba import jit
@jit
#now repeat the function with a different name so we can compare them
def Gamma_linalg_tiny_numba(T,c_x,q_alpha, q_A):
q_tau = q_A/T
q_G = np.exp(-(q_alpha*q_tau))
l_D = ((1/((q_G.T) @ c_x)).T)
q_E = (q_tau*q_G) * l_D
gamma = np.exp(((q_E+(q_E.T))-(((q_G * l_D) * (c_x.T)) @ (q_E.T))) @ c_x)
return gamma
%%timeit
# ~370 micro seconds per loop here (350 subtracting the random number generation, versus 420 thats not much acceleration)
N=3
x=np.random.rand(N,1)
x=x/sum(x)
alpha=np.random.rand(N,N)
A=np.random.rand(N,N)
T=(np.random.rand(1,1)+.5)*273
_ = Gamma_numba(
T,
x,
alpha,
A)
The slowest run took 34.63 times longer than the fastest. This could mean that an intermediate result is being cached. 1000 loops, best of 3: 388 µs per loop
%%timeit
# ~34 micro seconds per loop here (14 subtracting the random number generation, versus 32 thats approximately half)
N=3
x=np.random.rand(N,1)
x=x/sum(x)
alpha=np.random.rand(N,N)
A=np.random.rand(N,N)
T=(np.random.rand(1,1)+.5)*273
_ = Gamma_linalg_tiny_numba(
T,
x,
alpha,
A)
The slowest run took 50.88 times longer than the fastest. This could mean that an intermediate result is being cached. 10000 loops, best of 3: 36.5 µs per loop
- gamma1 vs (x1,x2 | x3=1-x1-x2)
- gamma2 vs (x1,x2 | x3=1-x1-x2)
- gamma3 vs (x1,x2 | x3=1-x1-x2)
Liq-Liq Equilibria flash algorithm
2d x1-x2-x3 triangle plots of phase equilibria envelope at given T and P
Optimizing code with cython and numba, comparison with fortran and c
Problem inversion: Regression of alpha and tau parameters from experimental data.