2 years ago
#23678
Jaka Orehek
optimization problem scipy.optimize.minimize
I'm having a hard time to optimize my system... I'm working on an antisolvent continuous crystallization in a tubular system and my goal is to find antisolvent flow rates at each antisolvent injection locations throughout the system length to achieve desired median crystal size. So initial concentration, solvent flow rate, temperature, injection locations are constant, only flow rates vary.
I've written a code which is posted below - I know, it is a bit messy... The idea behind it is that initial guesses of flow rates are given initial_guesses = [25, 0, 45, 0, 0] which are optimized based on least-squares optimization procedure:
desired_size = 50 # um
objective = (desired_size - d50[-1])**2
However, I applied some restriction to the optimization proces:
- If the solute concentration at the end of the sector (sectors are tube section between antisolvent injection locations) is equal or higher than metastable concentration * 0.333 for given antisolvent composition, the antisolvent flow rate at the following injection point is 0 - to prevent further nucleation.
- If the added antisolvent dissolved the solution the the degree which causes the crystals to dissolve
c_raz[l,z] < c_s[l], than the initial guess is being reduced for 1 mL/mininitial_guesses[l] -= 1until the antisolvent flow rate at the concerned injection point is not small enough to prevent dissolving (this is stopped if the flow rate reach value <= 0).
However, it was found out that these restrictions are not working properly. Does anyone have better idea how should it be written to achieve better optimization?
The code:
from time import time
startTime = time()
#%% parameters
# length and diameter of the system
L_reaktor = (41.7+41.5+42.5+41.8+60)*0.01 # m
D_reaktor = 3e-3 # m
# solvent flow rate
fi_vtok_mL_min = 50 # mL/min
fi_v_top = fi_vtok_mL_min * 1e-6 / 60 # m^3/s
del_top_vtok = 0.55 # ethanol fraction in solvent, vol%
cistost_top = 0.96 # vol%
# antisolvent flow rates and injection points
lokacija = [0, 0.417, 0.832, 1.257, 1.675] # m
del_prot_pritok = 1 # water fraction in antisolvent, vol%
# solute concentration in the solvent at the inlet
y = 110 # g/L = kg/m^3 # 118.3
# solute mol mass
M_m = 122.123e-3 # kg/mol
# solute density
rho_krist = 1265.9 # kg/m^3
# antisolvent and solvent density at 20°C
rho_voda = 1000 # kg/m3
rho_etoh = 789 # kg/m3
# antisolvent and solvent viskozity at 20°C
mi_voda = 0.0010005
mi_etoh = 0.0011841
#%% liberies
from numpy import zeros, interp, linspace, log10, pi, concatenate
from matplotlib.pyplot import figure, plot, close, xlabel, ylabel, rc, xscale, legend, xlim
from scipy.optimize import minimize
close('all')
#%% other parameters
# crystal size range
L_min = 1 # um
L_max = 180 # um
# diskretization
Nx = 20
Nz_0 = 1e4
Nz = int(Nz_0)
L = zeros(len(lokacija))
for i in range(len(lokacija)-1):
L[i] = lokacija[i+1]-lokacija[i]
L[-1] = L_reaktor - lokacija[-1]
Nz = int(Nz_0)
dz = L/Nz_0
# log size distribution
min_L = log10(L_min) # log um
max_L = log10(L_max) # log um
velikosti = 10**(linspace(min_L,max_L,Nx)) * 1e-6 # m
Vi = 4 * pi * (velikosti/2)**3 / 3 # m^3
Ai = 4 * pi * (velikosti/2)**2 # m^2
dx = velikosti[1:] - velikosti[0:-1] # m
velikosti_um = velikosti * 1e6 # um
# system volume
A_reaktor = (pi * D_reaktor**2 / 4)
V_reaktor = L_reaktor * A_reaktor
# bilance part volume
V_bil = A_reaktor * dz
#%% kinetics
# nucleation
k_b1 = 1.1637e+11
Nb1 = 2.8994
# crystal growth
k_G = 1.46961e-05
Ng = 0.125455
# breakage
# C_b = 2.36235e-008 * 1.4 # breakage
# K4 = 4.144022492245547e-08 #C_b * 2**(3/2) * (3.14)**0.5/4 # breakage
#%% experiments
# metastable limit
V_del_izpad = [0,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8, 1]
c_izpad = [433.3,437.7550769,425.83833,401.5159019,368.36608,329.5792969,287.95813,245.9173019,205.48368,168.2962769,135.60625,108.2769019,86.78368,71.21417688,61.26813,56.25742188,55.10608,54.000006]
# solubility
V_del_nasic = [0,0.05,0.1,0.15,0.2,0.25,0.3,0.35,0.4,0.45,0.5,0.55,0.6,0.65,0.7,0.75,0.8, 1]
c_nasic = [383.3,387.7550769,375.83833,351.5159019,318.36608,279.5792969,237.95813,195.9173019,155.48368,118.2962769,85.60625,58.27690188,36.78368,21.21417688,11.26813,6.257421875,5.10608,2.9]
#%% REGRESsION ######################################################################################
# prealocation
Nfi = len(lokacija)
V_krist_na_s = zeros([Nfi,Nz+1]) # m^3/s
n_krist_na_s = zeros([Nfi,Nz+1]) # mol/s
c_raz = zeros([Nfi,Nz+1]) # mol/m^3
y_raz = zeros([Nfi,Nz+1])
n = zeros([Nfi,Nz+1,Nx])
n_kon = zeros([Nfi,Nx])
n_kom = zeros([Nfi,Nx+1])
n_kom_del = zeros([Nfi,Nx+1])
d10 = zeros(Nfi)
d50 = zeros(Nfi)
d90 = zeros(Nfi)
pretok_sum = zeros(Nfi)
fi_v_prot_sum = zeros(Nfi)
V_del_prot = zeros(Nfi)
V_del_top = zeros(Nfi)
m_mix = zeros(Nfi)
m_prot = zeros(Nfi)
m_top = zeros(Nfi)
rho_top = zeros(Nfi)
m_top_i = zeros(Nfi)
m_prot_i = zeros(Nfi)
m_krist_na_s = zeros([Nfi,Nz])
V_krist_na_s = zeros([Nfi,Nz]) # m^3/s
n_krist_na_s = zeros([Nfi,Nz]) # mol/s
fi_v_prot = zeros(Nfi)
c_s_0 = zeros(Nfi)
c_s = zeros(Nfi)
c_nuk_0 = zeros(Nfi)
c_nuk = zeros(Nfi)
rho_mix = zeros(Nfi)
mi_mix = zeros(Nfi)
Re = zeros(Nfi)
De = zeros(Nfi)
v = zeros(Nfi)
fi_v_prot_mlmin = zeros(Nfi)
# initial guesses
initial_guesses = [25, 0, 45, 0, 0]
# inlet solubility
c_s_inlet = y / M_m # mol/m^3
# inlet molarity
n_0 = c_s_inlet * fi_v_top # mol/s vstopa
# calculation
def ode_enacbe(initial_guesses):
for l in range(Nfi):
# initial conditions:
if l == 0:
n[0,0,:] = 0
else:
n[l,0,:] = n[l-1,-1,:]
#from mL/min to m3/s
fi_v_prot_mlmin[l] = initial_guesses[l]
fi_v_prot[l] = fi_v_prot_mlmin[l] * 1e-6 / 60 # m^3/s
if l > 0:
if c_raz[l-1,-2] >= c_nuk[l-1] * 0.333:
initial_guesses[l] = 0
fi_v_prot_mlmin[l] = 0
fi_v_prot[l] = 0
# total flow rate and velocity after each injection point
fi_v_prot_sum[l] = sum(fi_v_prot[0:l+1])
pretok_sum[l] = fi_v_top + fi_v_prot_sum[l]
v[l] = pretok_sum[l] / A_reaktor # m/s
# solvent/antisolvnt ratios
V_del_prot[l] = (fi_v_top * (1-del_top_vtok) + (fi_v_prot_sum[l] * (del_prot_pritok))) / pretok_sum[l] #1 - V_del_top # V_del_prot = sum(fi_v_prot[0:k+1] * del_prot_pritok) / pretok_sum # kontrola
V_del_top[l] = 1 - V_del_prot[l]
# solubility kg/m^3, mol/m^3
c_s_0[l] = interp(V_del_prot[l],V_del_nasic,c_nasic) # kg/m^3
c_s[l] = c_s_0[l] / M_m # mol/m^3
# metastable limit
c_nuk_0[l] = interp(V_del_prot[l],V_del_izpad,c_izpad) # kg/m^3
c_nuk[l] = c_nuk_0[l] / M_m # mol/m^3
# density and viscosity of solution
rho_mix[l] = V_del_prot[l] * rho_voda + V_del_top[l] * rho_etoh
mi_mix[l] = V_del_prot[l] * mi_voda + V_del_top[l] * mi_etoh
# Reynold' s and Dean's numbers
Re[l] = rho_mix[l] * v[l] * D_reaktor / mi_mix[l]
De[l] = Re[l] * 0.1**0.5
# composition before and after of injection points
if l == 0:
m_top_i[l] = rho_etoh * del_top_vtok * fi_v_top
m_prot_i[l] = rho_voda * (1-del_top_vtok) * fi_v_top
else:
m_top_i[l] = rho_etoh * V_del_top[l-1] * pretok_sum[l-1]
m_prot_i[l] = rho_voda * V_del_prot[l-1] * pretok_sum[l-1]
m_prot[l] = V_del_prot[l] * pretok_sum[l] * rho_voda
m_top[l] = V_del_top[l] * pretok_sum[l] * rho_etoh
# Euler's method
for z in range(0,Nz):
# mass balance
m_krist_na_s[l,z] = sum(n[l,z,:] * Vi) * fi_v_top * rho_krist # kg/s
# concentration balance
c_raz[l,z] = (n_0 - (m_krist_na_s[l,z] / M_m)) / pretok_sum[l] # mol/m^3
y_raz[l,z] = c_raz[l,z]*M_m
# supersaturation
S = c_raz[l,z] / c_s[l]
A = (m_prot[l]/m_top[l]) - (m_prot_i[l]/m_top_i[l]) # for nucleation
if c_raz[l,z] > c_s[l] and c_raz[l,z] < c_nuk[l]*0.333:
for i in range(Nx):
if i == 0:
n[l,z+1,i] = n[l,z,i] + 1/v[l] * dz[l] * (- (k_G * (S-1)**Ng) * n[l,z,i]/dx[i]) # št./m^3
elif i < Nx-1:
n[l,z+1,i] = n[l,z,i] + 1/v[l] * dz[l] * ((k_G * (S-1)**Ng) * n[l,z,i-1]/dx[i-1] - (k_G * (S-1)**Ng) * n[l,z,i]/dx[i]) # št./m^3
else:
n[l,z+1,i] = n[l,z,i] + 1/v[l] * dz[l] * (k_G * (S-1)**Ng) * (n[l,z,i-1]/dx[i-1]) # št./m^3
elif c_raz[l,z] >= c_nuk[l]*0.333:
B1 = k_b1 * (A)**Nb1 # nuclei/m^3 s
for i in range(Nx):
if i == 0:
n[l,z+1,i] = n[l,z,i] + 1/v[l] * dz[l] * (- (k_G * (S-1)**Ng) * n[l,z,i]/dx[i] + B1) # št./m^3
elif i < Nx-1:
n[l,z+1,i] = n[l,z,i] + 1/v[l] * dz[l] * ((k_G * (S-1)**Ng) * n[l,z,i-1]/dx[i-1] - (k_G * (S-1)**Ng) * n[l,z,i]/dx[i]) # št./m^3
else:
n[l,z+1,i] = n[l,z,i] + 1/v[l] * dz[l] * (k_G * (S-1)**Ng) * (n[l,z,i-1]/dx[i-1]) # št./m^3
elif c_raz[l,z] <= c_s[l]:
n[l,z+1,:] = n[l,z,:]
else:
n[l,z+1,:] = n[l,z,:]
if l > 0 and c_raz[l,0] < c_s[l]:
while c_raz[l,0] < c_s[l]:
initial_guesses[l] -= 1
if initial_guesses[l] <= 0:
initial_guesses[l] = 0
#from mL/min to m3/s
fi_v_prot_mlmin[l] = initial_guesses[l]
fi_v_prot[l] = fi_v_prot_mlmin[l] * 1e-6 / 60 # m^3/s
# total flow rate and velocity after each injection point
fi_v_prot_sum[l] = sum(fi_v_prot[0:l+1])
pretok_sum[l] = fi_v_top + fi_v_prot_sum[l]
v[l] = pretok_sum[l] / A_reaktor # m/s
# solvent/antisolvnt ratios
V_del_prot[l] = (fi_v_top * (1-del_top_vtok) + (fi_v_prot_sum[l] * (del_prot_pritok))) / pretok_sum[l] #1 - V_del_top # V_del_prot = sum(fi_v_prot[0:k+1] * del_prot_pritok) / pretok_sum # kontrola
V_del_top[l] = 1 - V_del_prot[l]
# solubility kg/m^3, mol/m^3
c_s_0[l] = interp(V_del_prot[l],V_del_nasic,c_nasic) # kg/m^3
c_s[l] = c_s_0[l] / M_m # mol/m^3
# metastable limit
c_nuk_0[l] = interp(V_del_prot[l],V_del_izpad,c_izpad) # kg/m^3
c_nuk[l] = c_nuk_0[l] / M_m # mol/m^3
# density and viscosity of solution
rho_mix[l] = V_del_prot[l] * rho_voda + V_del_top[l] * rho_etoh
mi_mix[l] = V_del_prot[l] * mi_voda + V_del_top[l] * mi_etoh
# Reynold' s and Dean's numbers
Re[l] = rho_mix[l] * v[l] * D_reaktor / mi_mix[l]
De[l] = Re[l] * 0.1**0.5
# composition before and after of injection points
m_top_i[l] = rho_etoh * V_del_top[l-1] * pretok_sum[l-1]
m_prot_i[l] = rho_voda * V_del_prot[l-1] * pretok_sum[l-1]
m_prot[l] = V_del_prot[l] * pretok_sum[l] * rho_voda
m_top[l] = V_del_top[l] * pretok_sum[l] * rho_etoh
# Euler's method
for z in range(0,Nz):
# mass balance
m_krist_na_s[l,z] = sum(n[l,z,:] * Vi) * fi_v_top * rho_krist # kg/s
# concentration balance
c_raz[l,z] = (n_0 - (m_krist_na_s[l,z] / M_m)) / pretok_sum[l] # mol/m^3
y_raz[l,z] = c_raz[l,z]*M_m
# supersaturation
S = c_raz[l,z] / c_s[l]
A = (m_prot[l]/m_top[l]) - (m_prot_i[l]/m_top_i[l])
if c_raz[l,z] > c_s[l] and c_raz[l,z] < c_nuk[l]*0.333:
for i in range(Nx):
if i == 0:
n[l,z+1,i] = n[l,z,i] + 1/v[l] * dz[l] * (- (k_G * (S-1)**Ng) * n[l,z,i]/dx[i]) # št./m^3
elif i < Nx-1:
n[l,z+1,i] = n[l,z,i] + 1/v[l] * dz[l] * ((k_G * (S-1)**Ng) * n[l,z,i-1]/dx[i-1] - (k_G * (S-1)**Ng) * n[l,z,i]/dx[i]) # št./m^3
else:
n[l,z+1,i] = n[l,z,i] + 1/v[l] * dz[l] * (k_G * (S-1)**Ng) * (n[l,z,i-1]/dx[i-1]) # št./m^3
elif c_raz[l,z] >= c_nuk[l]*0.333:
B1 = k_b1 * (A)**Nb1 # nuclei/m^3 s
for i in range(Nx):
if i == 0:
n[l,z+1,i] = n[l,z,i] + 1/v[l] * dz[l] * (- (k_G * (S-1)**Ng) * n[l,z,i]/dx[i] + B1) # št./m^3
elif i < Nx-1:
n[l,z+1,i] = n[l,z,i] + 1/v[l] * dz[l] * ((k_G * (S-1)**Ng) * n[l,z,i-1]/dx[i-1] - (k_G * (S-1)**Ng) * n[l,z,i]/dx[i]) # št./m^3
else:
n[l,z+1,i] = n[l,z,i] + 1/v[l] * dz[l] * (k_G * (S-1)**Ng) * (n[l,z,i-1]/dx[i-1]) # št./m^3
else:
n[l,z+1,:] = n[l,z,:]
break
else:
# from mL/min to m3/s
fi_v_prot_mlmin[l] = initial_guesses[l]
fi_v_prot[l] = fi_v_prot_mlmin[l] * 1e-6 / 60 # m^3/s
if l > 0:
if c_raz[l-1,-2] >= c_nuk[l-1] * 0.333:
initial_guesses[l] = 0
fi_v_prot_mlmin[l] = 0
fi_v_prot[l] = 0
# total flow rate and velocity after each injection point
fi_v_prot_sum[l] = sum(fi_v_prot[0:l+1])
pretok_sum[l] = fi_v_top + fi_v_prot_sum[l]
v[l] = pretok_sum[l] / A_reaktor # m/s
# solvent/antisolvnt ratios
V_del_prot[l] = (fi_v_top * (1-del_top_vtok) + (fi_v_prot_sum[l] * (del_prot_pritok))) / pretok_sum[l] #1 - V_del_top # V_del_prot = sum(fi_v_prot[0:k+1] * del_prot_pritok) / pretok_sum # kontrola
V_del_top[l] = 1 - V_del_prot[l]
# solubility kg/m^3, mol/m^3
c_s_0[l] = interp(V_del_prot[l],V_del_nasic,c_nasic) # kg/m^3
c_s[l] = c_s_0[l] / M_m # mol/m^3
# metastable limit
c_nuk_0[l] = interp(V_del_prot[l],V_del_izpad,c_izpad) # kg/m^3
c_nuk[l] = c_nuk_0[l] / M_m # mol/m^3
# density and viscosity of solution
rho_mix[l] = V_del_prot[l] * rho_voda + V_del_top[l] * rho_etoh
mi_mix[l] = V_del_prot[l] * mi_voda + V_del_top[l] * mi_etoh
# Reynold' s and Dean's numbers
Re[l] = rho_mix[l] * v[l] * D_reaktor / mi_mix[l]
De[l] = Re[l] * 0.1**0.5
# composition before and after of injection points
m_top_i[l] = rho_etoh * V_del_top[l-1] * pretok_sum[l-1]
m_prot_i[l] = rho_voda * V_del_prot[l-1] * pretok_sum[l-1]
m_prot[l] = V_del_prot[l] * pretok_sum[l] * rho_voda
m_top[l] = V_del_top[l] * pretok_sum[l] * rho_etoh
# izračun Euler
for z in range(0,Nz):
# mass balance
m_krist_na_s[l,z] = sum(n[l,z,:] * Vi) * fi_v_top * rho_krist # kg/s
# concentration balance
c_raz[l,z] = (n_0 - (m_krist_na_s[l,z] / M_m)) / pretok_sum[l] # mol/m^3
y_raz[l,z] = c_raz[l,z]*M_m
# supersaturation
S = c_raz[l,z] / c_s[l]
A = (m_prot[l]/m_top[l]) - (m_prot_i[l]/m_top_i[l])
if c_raz[l,z] > c_s[l] and c_raz[l,z] < c_nuk[l]*0.333:
for i in range(Nx):
if i == 0:
n[l,z+1,i] = n[l,z,i] + 1/v[l] * dz[l] * (- (k_G * (S-1)**Ng) * n[l,z,i]/dx[i]) # št./m^3
elif i < Nx-1:
n[l,z+1,i] = n[l,z,i] + 1/v[l] * dz[l] * ((k_G * (S-1)**Ng) * n[l,z,i-1]/dx[i-1] - (k_G * (S-1)**Ng) * n[l,z,i]/dx[i]) # št./m^3
else:
n[l,z+1,i] = n[l,z,i] + 1/v[l] * dz[l] * (k_G * (S-1)**Ng) * (n[l,z,i-1]/dx[i-1]) # št./m^3
elif c_raz[l,z] >= c_nuk[l]*0.333:
B1 = k_b1 * (A)**Nb1 # nuclei/m^3 s
for i in range(Nx):
if i == 0:
n[l,z+1,i] = n[l,z,i] + 1/v[l] * dz[l] * (- (k_G * (S-1)**Ng) * n[l,z,i]/dx[i] + B1) # št./m^3
elif i < Nx-1:
n[l,z+1,i] = n[l,z,i] + 1/v[l] * dz[l] * ((k_G * (S-1)**Ng) * n[l,z,i-1]/dx[i-1] - (k_G * (S-1)**Ng) * n[l,z,i]/dx[i]) # št./m^3
else:
n[l,z+1,i] = n[l,z,i] + 1/v[l] * dz[l] * (k_G * (S-1)**Ng) * (n[l,z,i-1]/dx[i-1]) # št./m^3
else:
n[l,z+1,:] = n[l,z,:]
# normalized comulative crystal number
n_kon[l] = n[l,-1,:]
n_kom = zeros([Nfi,Nx+1])
for i in range(Nx):
n_kom[l,i+1] = (n_kom[l,i]+n_kon[l,i])
n_kom_del[l] = n_kom[l,0:Nx+1]/sum(n_kon[l])
d10[l] = interp(0.1,n_kom_del[l,0:Nx],velikosti_um)
d50[l] = interp(0.5,n_kom_del[l,0:Nx],velikosti_um)
d90[l] = interp(0.9,n_kom_del[l,0:Nx],velikosti_um)
# target function
desired_size = 50 # um
objective = (desired_size - d50[-1])**2
return objective
# Solution
parametri = minimize(ode_enacbe, initial_guesses, method='Nelder-Mead', bounds =((25,200), (0,40),
(0,45), (0,40), (0,40)))
executionTime = (time() - startTime)/60
print('Execution time in minutes: ' +str(executionTime))
python
optimization
minimize
scipy-optimize-minimize
0 Answers
Your Answer