Fit Pure Fluid ParametersΒΆ
This example shows how to use the new fitting class of version 0.21 of teqp to fit the model parameters for the SAFT-VR-Mie model based on fitting pseudo-experimental data obtained from a reference equation of state.
[1]:
import teqp
import numpy as np
import matplotlib.pyplot as plt
import CoolProp.CoolProp as CP
import scipy.optimize
display(teqp.__version__)
nonpolar = {
"kind": "SAFT-VR-Mie",
"model": {
"coeffs": [
{
"name": "R32",
"BibTeXKey": "Bell",
"m": 1.2476268271391935,
"sigma_m": 3.6080717234117107e-10,
"epsilon_over_k": 172.53065054286867,
"lambda_r": 14.634722358167384,
"lambda_a": 6
}
]
}
}
template = {
'kind': 'genericSAFT',
'model': {
'nonpolar': nonpolar
}
}
# NOTE: '/' inside field names MUST be escaped as ~1; see https://datatracker.ietf.org/doc/html/rfc6901#section-3
pointers = [
'/model/nonpolar/model/coeffs/0/m',
'/model/nonpolar/model/coeffs/0/sigma_m',
'/model/nonpolar/model/coeffs/0/epsilon_over_k',
'/model/nonpolar/model/coeffs/0/lambda_r',
'/model/nonpolar/model/coeffs/0/lambda_a'
]
x0 = [1.5, 3e-10, 150, 19, 5.7]
bounds = [(1,5), (2e-10,5e-10), (100,400), (12,50), (5.1, 6.0)]
'0.21.0'
[2]:
FLD = 'Methane'
ppo = teqp.PureParameterOptimizer(template, pointers)
Ts = np.linspace(100, 170, 10)
# Generate some artificial pseudo-experimental data from
# the reference EOS as implemented in CoolProp
for T in Ts:
pt = teqp.SatRhoLPWPoint()
rhoL, rhoV = [CP.PropsSI('Dmolar','Q',Q,'T',T,FLD) for Q in [0,1]]
p = CP.PropsSI('P','Q',0,'T',T,FLD)
w = CP.PropsSI('speed_of_sound','Q',0,'T',T,FLD)
pt.T = T
# Measurands (here, pseudo-experimental values)
pt.p_exp = p
pt.rhoL_exp = rhoL
pt.w_exp = w
# Additional parameters
pt.rhoL_guess = rhoL
pt.rhoV_guess = rhoV
pt.R = 8.31446261815324
pt.M = CP.PropsSI('molemass',FLD)
AS = CP.AbstractState('HEOS',FLD)
AS.update(CP.DmolarT_INPUTS, 1e-10, T)
pt.Ao20 = AS.tau()**2*AS.d2alpha0_dTau2() # -cv0/R
# Weights (multiplied by 100 to put on a percentage basis)
pt.weight_p = 1*100
pt.weight_rho = 1*100
pt.weight_w = 0.25*100
ppo.add_one_contribution(pt)
def cost_function(x):
# return ppo.cost_function_threaded(x, 10) # This is an option if you have lots of threads
return ppo.cost_function(x)
r = scipy.optimize.differential_evolution(cost_function, bounds=bounds, disp=True, maxiter=10000, popsize=8)
print(r)
x = r.x
model = teqp.make_model(ppo.build_JSON(x))
Tc, rhoc = model.solve_pure_critical(400, 5000)
print(Tc, rhoc)
anc = teqp.build_ancillaries(model, Tc, rhoc, 0.5*Tc)
def _get_SOS(model, T, rho, z, *, R, M, Ao20):
""" Helper function to calculate speed of sound """
Ar0n = model.get_Ar02n(T, rho, z)
Ar01 = Ar0n[1]; Ar02 = Ar0n[2]
Ar11 = model.get_Ar11(T, rho, z)
Ar20 = model.get_Ar20(T, rho, z)
# M*w^2/(R*T) where w is the speed of sound
# from the definition w = sqrt(dp/drho|s)
Mw2RT = 1 + 2*Ar01 + Ar02 - (1 + Ar01 - Ar11)**2/(Ao20 + Ar20)
if Mw2RT < 0:
return 1e6
w = (Mw2RT*R*T/M)**0.5
return w
TcREF = CP.PropsSI('Tcrit', FLD)
Tsverify = np.linspace(100, TcREF*0.99999, 1000)
RHOL, RHOV, PPP, WWW = [],[],[],[]
z = np.array([1.0])
for T in Tsverify:
rhoL, rhoV = model.pure_VLE_T(T, anc.rhoL(T)*1.01, anc.rhoV(T)*0.9, 10)
pL = rhoL*model.get_R(z)*T*(1+model.get_Ar01(T, rhoL, z))
pV = rhoV*model.get_R(z)*T*(1+model.get_Ar01(T, rhoV, z))
RHOL.append(rhoL)
RHOV.append(rhoV)
PPP.append(pL)
AS = CP.AbstractState('HEOS',FLD)
AS.update(CP.DmolarT_INPUTS, 1e-10, T)
Ao20 = AS.d2alpha0_dTau2()*AS.tau()**2
WWW.append(_get_SOS(model, T, rhoL, z, R=8.314462618,M=CP.PropsSI('molemass',FLD),Ao20=Ao20))
# Plot the T-rho for VLE
line, = plt.plot(np.array(RHOL), Tsverify)
plt.plot(np.array(RHOV), Tsverify, color=line.get_color())
for Q in [0, 1]:
D = CP.PropsSI('Dmolar','T',Tsverify,'Q',Q,FLD)
plt.plot(D, Tsverify, lw=2, color='k')
plt.gca().set(xlabel=r'$\rho$ / mol/m$^3$', ylabel='$T$ / K')
# And a deviation plot, much closer to the critical point
# than the fitting region
fig, (ax1, ax2, ax3) = plt.subplots(3,1, figsize=(10,5), sharex=True)
ax1.plot(Tsverify, (np.array(PPP)/CP.PropsSI('P','T',Tsverify,'Q',0,FLD)-1)*100)
ax1.set(ylabel=r'$(p_{fit}/p_{\rm pexp}-1)\times 100$')
ax2.plot(Tsverify, (np.array(RHOL)/CP.PropsSI('Dmolar','T',Tsverify,'Q',0,FLD)-1)*100)
ax2.set(ylabel=r'$(\rho_{fit}/\rho_{\rm pexp}-1)\times 100$')
ax3.plot(Tsverify, (np.array(WWW)/CP.PropsSI('speed_of_sound','T',Tsverify,'Q',0,FLD)-1)*100)
ax3.set(ylabel=r'$(w_{fit}/w_{\rm pexp}-1)\times 100$', xlabel='$T$ / K')
differential_evolution step 1: f(x)= 5818.421956412952
differential_evolution step 2: f(x)= 499.1768555361148
differential_evolution step 3: f(x)= 487.6265883142837
differential_evolution step 4: f(x)= 457.7391656468499
differential_evolution step 5: f(x)= 457.7391656468499
differential_evolution step 6: f(x)= 457.7391656468499
differential_evolution step 7: f(x)= 418.0474196262808
differential_evolution step 8: f(x)= 418.0474196262808
differential_evolution step 9: f(x)= 418.0474196262808
differential_evolution step 10: f(x)= 338.1635147679501
differential_evolution step 11: f(x)= 338.1635147679501
differential_evolution step 12: f(x)= 338.1635147679501
differential_evolution step 13: f(x)= 338.1635147679501
differential_evolution step 14: f(x)= 338.1635147679501
differential_evolution step 15: f(x)= 338.1635147679501
differential_evolution step 16: f(x)= 330.10096482774765
differential_evolution step 17: f(x)= 330.10096482774765
differential_evolution step 18: f(x)= 330.10096482774765
differential_evolution step 19: f(x)= 330.10096482774765
differential_evolution step 20: f(x)= 330.10096482774765
differential_evolution step 21: f(x)= 275.1544998794892
differential_evolution step 22: f(x)= 275.1544998794892
differential_evolution step 23: f(x)= 241.81895513378925
differential_evolution step 24: f(x)= 241.81895513378925
differential_evolution step 25: f(x)= 241.81895513378925
differential_evolution step 26: f(x)= 241.81895513378925
differential_evolution step 27: f(x)= 241.81895513378925
differential_evolution step 28: f(x)= 241.81895513378925
differential_evolution step 29: f(x)= 241.81895513378925
differential_evolution step 30: f(x)= 241.81895513378925
differential_evolution step 31: f(x)= 241.81895513378925
differential_evolution step 32: f(x)= 241.81895513378925
differential_evolution step 33: f(x)= 241.81895513378925
differential_evolution step 34: f(x)= 241.81895513378925
differential_evolution step 35: f(x)= 237.2929366395986
differential_evolution step 36: f(x)= 237.2929366395986
differential_evolution step 37: f(x)= 237.2929366395986
differential_evolution step 38: f(x)= 237.2929366395986
differential_evolution step 39: f(x)= 237.2929366395986
differential_evolution step 40: f(x)= 230.68631809484037
differential_evolution step 41: f(x)= 230.20609180473951
differential_evolution step 42: f(x)= 209.69348566431282
differential_evolution step 43: f(x)= 209.69348566431282
differential_evolution step 44: f(x)= 209.69348566431282
differential_evolution step 45: f(x)= 182.23116564860715
differential_evolution step 46: f(x)= 182.23116564860715
differential_evolution step 47: f(x)= 146.77960931385533
differential_evolution step 48: f(x)= 146.77960931385533
differential_evolution step 49: f(x)= 144.65095762969202
differential_evolution step 50: f(x)= 122.62697814799093
differential_evolution step 51: f(x)= 63.07678081956036
differential_evolution step 52: f(x)= 63.07678081956036
differential_evolution step 53: f(x)= 63.07678081956036
differential_evolution step 54: f(x)= 63.07678081956036
differential_evolution step 55: f(x)= 63.07678081956036
differential_evolution step 56: f(x)= 63.07678081956036
differential_evolution step 57: f(x)= 46.53073192085213
differential_evolution step 58: f(x)= 46.53073192085213
differential_evolution step 59: f(x)= 46.53073192085213
differential_evolution step 60: f(x)= 46.53073192085213
differential_evolution step 61: f(x)= 46.53073192085213
differential_evolution step 62: f(x)= 46.53073192085213
differential_evolution step 63: f(x)= 39.1040859298758
differential_evolution step 64: f(x)= 39.1040859298758
differential_evolution step 65: f(x)= 39.1040859298758
differential_evolution step 66: f(x)= 39.1040859298758
differential_evolution step 67: f(x)= 39.1040859298758
differential_evolution step 68: f(x)= 39.1040859298758
differential_evolution step 69: f(x)= 39.1040859298758
differential_evolution step 70: f(x)= 39.1040859298758
differential_evolution step 71: f(x)= 35.91213962638507
differential_evolution step 72: f(x)= 35.91213962638507
differential_evolution step 73: f(x)= 35.91213962638507
differential_evolution step 74: f(x)= 31.991928955829998
differential_evolution step 75: f(x)= 25.69596469618272
differential_evolution step 76: f(x)= 25.69596469618272
differential_evolution step 77: f(x)= 25.69596469618272
differential_evolution step 78: f(x)= 25.69596469618272
differential_evolution step 79: f(x)= 25.69596469618272
differential_evolution step 80: f(x)= 25.69596469618272
differential_evolution step 81: f(x)= 25.69596469618272
differential_evolution step 82: f(x)= 25.69596469618272
differential_evolution step 83: f(x)= 25.69596469618272
differential_evolution step 84: f(x)= 25.69596469618272
differential_evolution step 85: f(x)= 25.40977267008008
differential_evolution step 86: f(x)= 25.08654968764155
differential_evolution step 87: f(x)= 21.732540180244946
differential_evolution step 88: f(x)= 21.732540180244946
differential_evolution step 89: f(x)= 21.732540180244946
differential_evolution step 90: f(x)= 21.732540180244946
differential_evolution step 91: f(x)= 21.732540180244946
differential_evolution step 92: f(x)= 21.732540180244946
differential_evolution step 93: f(x)= 21.732540180244946
differential_evolution step 94: f(x)= 21.732540180244946
differential_evolution step 95: f(x)= 21.5950406845984
differential_evolution step 96: f(x)= 21.565973142204154
differential_evolution step 97: f(x)= 21.565973142204154
differential_evolution step 98: f(x)= 21.565973142204154
differential_evolution step 99: f(x)= 21.530760903539687
differential_evolution step 100: f(x)= 21.530760903539687
differential_evolution step 101: f(x)= 21.316727944016527
differential_evolution step 102: f(x)= 21.25499340116124
differential_evolution step 103: f(x)= 20.559872546488993
differential_evolution step 104: f(x)= 20.020765011191784
differential_evolution step 105: f(x)= 20.020765011191784
differential_evolution step 106: f(x)= 20.020765011191784
differential_evolution step 107: f(x)= 20.020765011191784
differential_evolution step 108: f(x)= 19.549324351102406
differential_evolution step 109: f(x)= 19.532246269269777
differential_evolution step 110: f(x)= 19.532246269269777
differential_evolution step 111: f(x)= 19.01797540712929
differential_evolution step 112: f(x)= 19.01797540712929
differential_evolution step 113: f(x)= 18.895398559399503
differential_evolution step 114: f(x)= 18.776806394241568
differential_evolution step 115: f(x)= 18.749622014848157
differential_evolution step 116: f(x)= 18.749622014848157
differential_evolution step 117: f(x)= 18.56057810916112
differential_evolution step 118: f(x)= 18.47984249770556
differential_evolution step 119: f(x)= 18.47984249770556
differential_evolution step 120: f(x)= 18.47984249770556
differential_evolution step 121: f(x)= 18.346929054564406
differential_evolution step 122: f(x)= 18.060391807258394
differential_evolution step 123: f(x)= 18.060391807258394
differential_evolution step 124: f(x)= 18.060391807258394
differential_evolution step 125: f(x)= 17.744114151073173
differential_evolution step 126: f(x)= 17.68251103024557
differential_evolution step 127: f(x)= 17.68251103024557
differential_evolution step 128: f(x)= 17.68251103024557
differential_evolution step 129: f(x)= 17.52261388130418
differential_evolution step 130: f(x)= 17.2817439273557
differential_evolution step 131: f(x)= 16.54892785495758
differential_evolution step 132: f(x)= 16.54892785495758
differential_evolution step 133: f(x)= 16.31867885756894
differential_evolution step 134: f(x)= 15.911061208426595
differential_evolution step 135: f(x)= 15.614991653902301
differential_evolution step 136: f(x)= 15.614991653902301
differential_evolution step 137: f(x)= 14.789691475222233
differential_evolution step 138: f(x)= 14.58866564278831
differential_evolution step 139: f(x)= 14.58866564278831
differential_evolution step 140: f(x)= 14.58866564278831
differential_evolution step 141: f(x)= 14.349150893719104
differential_evolution step 142: f(x)= 13.423324045866506
differential_evolution step 143: f(x)= 13.21096960040352
differential_evolution step 144: f(x)= 13.040709831794738
differential_evolution step 145: f(x)= 13.040709831794738
differential_evolution step 146: f(x)= 13.040709831794738
differential_evolution step 147: f(x)= 13.040709831794738
differential_evolution step 148: f(x)= 13.040709831794738
differential_evolution step 149: f(x)= 13.040709831794738
differential_evolution step 150: f(x)= 13.040709831794738
differential_evolution step 151: f(x)= 13.040709831794738
differential_evolution step 152: f(x)= 13.040709831794738
differential_evolution step 153: f(x)= 13.040709831794738
differential_evolution step 154: f(x)= 13.040709831794738
differential_evolution step 155: f(x)= 12.8860672106379
differential_evolution step 156: f(x)= 12.8860672106379
differential_evolution step 157: f(x)= 12.8860672106379
differential_evolution step 158: f(x)= 12.8860672106379
differential_evolution step 159: f(x)= 12.8860672106379
differential_evolution step 160: f(x)= 12.8860672106379
differential_evolution step 161: f(x)= 12.884275329003621
differential_evolution step 162: f(x)= 12.884275329003621
differential_evolution step 163: f(x)= 12.884275329003621
differential_evolution step 164: f(x)= 12.884275329003621
differential_evolution step 165: f(x)= 12.862938445406227
differential_evolution step 166: f(x)= 12.862938445406227
differential_evolution step 167: f(x)= 12.862938445406227
differential_evolution step 168: f(x)= 12.843987764028336
differential_evolution step 169: f(x)= 12.843987764028336
differential_evolution step 170: f(x)= 12.843987764028336
differential_evolution step 171: f(x)= 12.843987764028336
differential_evolution step 172: f(x)= 12.843987764028336
differential_evolution step 173: f(x)= 12.843987764028336
Polishing solution with 'L-BFGS-B'
message: Optimization terminated successfully.
success: True
fun: 12.843987764028336
x: [ 1.159e+00 3.550e-10 1.153e+02 1.284e+01
5.170e+00]
nit: 173
nfev: 7086
population: [[ 1.159e+00 3.550e-10 ... 1.284e+01 5.170e+00]
[ 1.151e+00 3.559e-10 ... 1.291e+01 5.169e+00]
...
[ 1.174e+00 3.531e-10 ... 1.260e+01 5.160e+00]
[ 1.147e+00 3.563e-10 ... 1.298e+01 5.156e+00]]
population_energies: [ 1.284e+01 1.299e+01 ... 1.335e+01 1.298e+01]
194.75971762562753 9153.1514317028
[2]:
[Text(0, 0.5, '$(w_{fit}/w_{\\rm pexp}-1)\\times 100$'),
Text(0.5, 0, '$T$ / K')]
