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.
In version 0.23, the parameter optimization objects were moved into their own submodule called paramopt
[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.23.1'
[2]:
FLD = 'Methane'
ppo = teqp.paramopt.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.paramopt.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)= 2640.2565791642164
differential_evolution step 2: f(x)= 962.3756948786976
differential_evolution step 3: f(x)= 962.3756948786976
differential_evolution step 4: f(x)= 962.3756948786976
differential_evolution step 5: f(x)= 317.8225079750391
differential_evolution step 6: f(x)= 317.8225079750391
differential_evolution step 7: f(x)= 317.8225079750391
differential_evolution step 8: f(x)= 317.8225079750391
differential_evolution step 9: f(x)= 317.8225079750391
differential_evolution step 10: f(x)= 265.3026660890965
differential_evolution step 11: f(x)= 265.3026660890965
differential_evolution step 12: f(x)= 265.3026660890965
differential_evolution step 13: f(x)= 265.3026660890965
differential_evolution step 14: f(x)= 265.3026660890965
differential_evolution step 15: f(x)= 265.3026660890965
differential_evolution step 16: f(x)= 153.92786840948762
differential_evolution step 17: f(x)= 153.92786840948762
differential_evolution step 18: f(x)= 153.92786840948762
differential_evolution step 19: f(x)= 153.92786840948762
differential_evolution step 20: f(x)= 153.92786840948762
differential_evolution step 21: f(x)= 153.92786840948762
differential_evolution step 22: f(x)= 153.92786840948762
differential_evolution step 23: f(x)= 153.92786840948762
differential_evolution step 24: f(x)= 143.4728945770965
differential_evolution step 25: f(x)= 115.17668169040033
differential_evolution step 26: f(x)= 115.17668169040033
differential_evolution step 27: f(x)= 115.17668169040033
differential_evolution step 28: f(x)= 115.17668169040033
differential_evolution step 29: f(x)= 115.17668169040033
differential_evolution step 30: f(x)= 115.1166132737544
differential_evolution step 31: f(x)= 115.1166132737544
differential_evolution step 32: f(x)= 98.88614368389673
differential_evolution step 33: f(x)= 98.88614368389673
differential_evolution step 34: f(x)= 98.88614368389673
differential_evolution step 35: f(x)= 98.88614368389673
differential_evolution step 36: f(x)= 74.45782154462813
differential_evolution step 37: f(x)= 74.45782154462813
differential_evolution step 38: f(x)= 73.0782298154598
differential_evolution step 39: f(x)= 54.50814101608117
differential_evolution step 40: f(x)= 28.236972307787255
differential_evolution step 41: f(x)= 28.236972307787255
differential_evolution step 42: f(x)= 28.236972307787255
differential_evolution step 43: f(x)= 28.236972307787255
differential_evolution step 44: f(x)= 28.236972307787255
differential_evolution step 45: f(x)= 28.236972307787255
differential_evolution step 46: f(x)= 28.236972307787255
differential_evolution step 47: f(x)= 28.236972307787255
differential_evolution step 48: f(x)= 28.236972307787255
differential_evolution step 49: f(x)= 28.236972307787255
differential_evolution step 50: f(x)= 27.807109540500456
differential_evolution step 51: f(x)= 27.807109540500456
differential_evolution step 52: f(x)= 27.807109540500456
differential_evolution step 53: f(x)= 27.807109540500456
differential_evolution step 54: f(x)= 27.807109540500456
differential_evolution step 55: f(x)= 27.807109540500456
differential_evolution step 56: f(x)= 27.807109540500456
differential_evolution step 57: f(x)= 27.807109540500456
differential_evolution step 58: f(x)= 27.807109540500456
differential_evolution step 59: f(x)= 27.807109540500456
differential_evolution step 60: f(x)= 27.807109540500456
differential_evolution step 61: f(x)= 23.498644693715875
differential_evolution step 62: f(x)= 23.498644693715875
differential_evolution step 63: f(x)= 23.498644693715875
differential_evolution step 64: f(x)= 23.498644693715875
differential_evolution step 65: f(x)= 23.498644693715875
differential_evolution step 66: f(x)= 23.498644693715875
differential_evolution step 67: f(x)= 23.498644693715875
differential_evolution step 68: f(x)= 21.16019457623664
differential_evolution step 69: f(x)= 21.16019457623664
differential_evolution step 70: f(x)= 19.620489270537846
differential_evolution step 71: f(x)= 19.286816629629598
differential_evolution step 72: f(x)= 19.286816629629598
differential_evolution step 73: f(x)= 19.286816629629598
differential_evolution step 74: f(x)= 19.286816629629598
differential_evolution step 75: f(x)= 19.286816629629598
differential_evolution step 76: f(x)= 19.286816629629598
differential_evolution step 77: f(x)= 19.286816629629598
differential_evolution step 78: f(x)= 19.286816629629598
differential_evolution step 79: f(x)= 19.286816629629598
differential_evolution step 80: f(x)= 19.286816629629598
differential_evolution step 81: f(x)= 19.055993887333994
differential_evolution step 82: f(x)= 19.055993887333994
differential_evolution step 83: f(x)= 19.055993887333994
differential_evolution step 84: f(x)= 19.055993887333994
differential_evolution step 85: f(x)= 19.055993887333994
differential_evolution step 86: f(x)= 19.055993887333994
differential_evolution step 87: f(x)= 19.055993887333994
differential_evolution step 88: f(x)= 19.055993887333994
differential_evolution step 89: f(x)= 19.055993887333994
differential_evolution step 90: f(x)= 19.055993887333994
differential_evolution step 91: f(x)= 19.055993887333994
differential_evolution step 92: f(x)= 19.055993887333994
differential_evolution step 93: f(x)= 19.055993887333994
differential_evolution step 94: f(x)= 19.055993887333994
differential_evolution step 95: f(x)= 19.055993887333994
differential_evolution step 96: f(x)= 19.055993887333994
differential_evolution step 97: f(x)= 18.987681202407188
differential_evolution step 98: f(x)= 18.987681202407188
differential_evolution step 99: f(x)= 18.355792749571748
differential_evolution step 100: f(x)= 18.355792749571748
differential_evolution step 101: f(x)= 18.355792749571748
differential_evolution step 102: f(x)= 18.355792749571748
differential_evolution step 103: f(x)= 18.355792749571748
differential_evolution step 104: f(x)= 17.732999573616617
differential_evolution step 105: f(x)= 17.732999573616617
differential_evolution step 106: f(x)= 17.732999573616617
differential_evolution step 107: f(x)= 17.396950587217354
differential_evolution step 108: f(x)= 17.396950587217354
differential_evolution step 109: f(x)= 17.396950587217354
differential_evolution step 110: f(x)= 17.396950587217354
differential_evolution step 111: f(x)= 17.283336737521402
differential_evolution step 112: f(x)= 17.047863176037033
differential_evolution step 113: f(x)= 17.047863176037033
differential_evolution step 114: f(x)= 17.047863176037033
differential_evolution step 115: f(x)= 16.70189296244538
differential_evolution step 116: f(x)= 16.70189296244538
differential_evolution step 117: f(x)= 16.70189296244538
differential_evolution step 118: f(x)= 16.679215439754024
differential_evolution step 119: f(x)= 16.585527215317267
differential_evolution step 120: f(x)= 16.585527215317267
differential_evolution step 121: f(x)= 16.440880985703856
differential_evolution step 122: f(x)= 16.440880985703856
differential_evolution step 123: f(x)= 16.440880985703856
differential_evolution step 124: f(x)= 16.440356282043112
differential_evolution step 125: f(x)= 16.440356282043112
differential_evolution step 126: f(x)= 16.440356282043112
differential_evolution step 127: f(x)= 16.440356282043112
differential_evolution step 128: f(x)= 16.440356282043112
differential_evolution step 129: f(x)= 16.440356282043112
differential_evolution step 130: f(x)= 16.3860525355149
differential_evolution step 131: f(x)= 16.255845102047942
differential_evolution step 132: f(x)= 16.255845102047942
differential_evolution step 133: f(x)= 16.237798228248426
differential_evolution step 134: f(x)= 16.04149141698661
differential_evolution step 135: f(x)= 16.04149141698661
differential_evolution step 136: f(x)= 16.04149141698661
differential_evolution step 137: f(x)= 15.861265882923531
differential_evolution step 138: f(x)= 15.683342142054734
differential_evolution step 139: f(x)= 15.631350161630701
differential_evolution step 140: f(x)= 15.241571118774463
differential_evolution step 141: f(x)= 14.907975717145758
differential_evolution step 142: f(x)= 14.907975717145758
differential_evolution step 143: f(x)= 14.683725820963293
differential_evolution step 144: f(x)= 14.274009967668707
differential_evolution step 145: f(x)= 14.274009967668707
differential_evolution step 146: f(x)= 13.945513832128142
differential_evolution step 147: f(x)= 13.059689802523561
differential_evolution step 148: f(x)= 12.887628269047887
differential_evolution step 149: f(x)= 12.887628269047887
differential_evolution step 150: f(x)= 12.887628269047887
differential_evolution step 151: f(x)= 12.887628269047887
differential_evolution step 152: f(x)= 12.887628269047887
differential_evolution step 153: f(x)= 12.887628269047887
differential_evolution step 154: f(x)= 12.887628269047887
differential_evolution step 155: f(x)= 12.887628269047887
differential_evolution step 156: f(x)= 12.887628269047887
differential_evolution step 157: f(x)= 12.887628269047887
differential_evolution step 158: f(x)= 12.78991787613433
differential_evolution step 159: f(x)= 12.78991787613433
differential_evolution step 160: f(x)= 12.78991787613433
differential_evolution step 161: f(x)= 12.78991787613433
differential_evolution step 162: f(x)= 12.78991787613433
differential_evolution step 163: f(x)= 12.78991787613433
differential_evolution step 164: f(x)= 12.78991787613433
differential_evolution step 165: f(x)= 12.78991787613433
differential_evolution step 166: f(x)= 12.78991787613433
differential_evolution step 167: f(x)= 12.78991787613433
differential_evolution step 168: f(x)= 12.78991787613433
differential_evolution step 169: f(x)= 12.78991787613433
differential_evolution step 170: f(x)= 12.78991787613433
differential_evolution step 171: f(x)= 12.78991787613433
Polishing solution with 'L-BFGS-B'
message: Optimization terminated successfully.
success: True
fun: 12.78991787613433
x: [ 1.185e+00 3.519e-10 1.145e+02 1.205e+01
5.303e+00]
nit: 171
nfev: 7006
population: [[ 1.185e+00 3.519e-10 ... 1.205e+01 5.303e+00]
[ 1.187e+00 3.515e-10 ... 1.201e+01 5.277e+00]
...
[ 1.185e+00 3.518e-10 ... 1.206e+01 5.300e+00]
[ 1.184e+00 3.519e-10 ... 1.204e+01 5.305e+00]]
population_energies: [ 1.279e+01 1.291e+01 ... 1.279e+01 1.279e+01]
194.8172792201637 9203.132811337495
[2]:
[Text(0, 0.5, '$(w_{fit}/w_{\\rm pexp}-1)\\times 100$'),
Text(0.5, 0, '$T$ / K')]
