Skip to content

Commit

Permalink
Maybe working?
Browse files Browse the repository at this point in the history
  • Loading branch information
@tsikes
tsikes committed Aug 9, 2021
1 parent 5e69b69 commit 0d89b68
Show file tree
Hide file tree
Showing 2 changed files with 62 additions and 34 deletions.
94 changes: 61 additions & 33 deletions src/calculate/optimize/fit_coeffs.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -566,7 +566,7 @@ def LPL_HPL_Fcent(self): # CURRENTLY A KNOWN BUG IF FULL LIMIT IS HELD CONSTA
arrhenius_coeffs, T, ln_k = self.fit_arrhenius_rates()

# set nlopt opt parameters
options = [{'algorithm': nlopt.GN_DIRECT, 'xtol_rel': 1E-4, 'ftol_rel': 1E-4, # GN_DIRECT_NOSCAL GN_DIRECT
options = [{'algorithm': nlopt.GN_CRS2_LM, 'xtol_rel': 1E-4, 'ftol_rel': 1E-4, # GN_DIRECT_NOSCAL GN_DIRECT
'initial_step': 1, 'max_eval': 2000},
#{'algorithm': nlopt.LN_COBYLA, 'xtol_rel': 1E-2, 'ftol_rel': 1E-2,
# 'initial_step': 1E-3, 'max_eval': 500},
Expand Down Expand Up @@ -597,12 +597,16 @@ def LPL_HPL_Fcent(self): # CURRENTLY A KNOWN BUG IF FULL LIMIT IS HELD CONSTA
M = np.array(M)

# set default p0, s, and p_bnds
p0 = np.array([39.7, 39.1, 38.6, 31.5, 31.5, 31.5])
p0 = np.concatenate((p0, -0.51*np.ones((1,6)).flatten())) # use 5th order polynomial for Fcent
p0 = np.array([3.6E7, 137.0, -7.0, 1.4E7, 34.5, 0.4])
#p0 = np.concatenate((p0, -0.51*np.ones((1,6)).flatten())) # use 5th order polynomial for Fcent
p0 = np.concatenate((p0, [-0.65, -0.63, -0.57, -0.49, -0.43, -0.36])) # use 5th order polynomial for Fcent
s = np.ones_like(p0)
#s = np.array([4184, 1.0, 1E-2, 4184, 1.0, 1E-2, *(np.ones((1,6)).flatten()*1E-1)])

p_bnds = np.ones((2, 6))*np.log([[1E-14], [1E80]])
p_bnds = np.subtract(p_bnds, p0[:6])/s[:6]
p_bnds = set_arrhenius_bnds(p0[0:3], default_arrhenius_coefNames)
p_bnds = np.concatenate((p_bnds, set_arrhenius_bnds(p0[3:6], default_arrhenius_coefNames)), axis=1)
p_bnds[:,1] = np.log(p_bnds[:,1]) # setting bnds to ln(A), need to do this better
p_bnds[:,4] = np.log(p_bnds[:,4]) # setting bnds to ln(A), need to do this better
p_bnds = np.concatenate((p_bnds, np.ones((2, 6))*np.log([[1E-3], [1.0]])), axis=1)

# calculate k_0, k_inf, Fcent
Expand Down Expand Up @@ -669,9 +673,10 @@ def LPL_HPL_Fcent(self): # CURRENTLY A KNOWN BUG IF FULL LIMIT IS HELD CONSTA

self.obj_vars = {'fcn': fit_fcn, 'T': T, 'T_arrhenius': T_arrhenius, 'T_Fcent': T_Fcent,
'M': M, 'ln_k': ln_k, 'coefs_0': p0, 's': s}
#s = self.obj_vars['s'] = self.calc_s(np.zeros_like(p0))

s = self.obj_vars['s'] = self.calc_s(np.zeros_like(p0))

self.opt = opt = nlopt.opt(opt_options['algorithm'], len(p0)) # AUGLAG
self.opt = opt = nlopt.opt(nlopt.AUGLAG, len(p0)) # AUGLAG

opt.set_min_objective(self.objective)
opt.set_maxeval(opt_options['max_eval'])
Expand All @@ -680,60 +685,63 @@ def LPL_HPL_Fcent(self): # CURRENTLY A KNOWN BUG IF FULL LIMIT IS HELD CONSTA
opt.set_xtol_rel(opt_options['xtol_rel'])
opt.set_ftol_rel(opt_options['ftol_rel'])

opt.set_lower_bounds(p_bnds[0])
opt.set_upper_bounds(p_bnds[1])
opt.set_lower_bounds((p_bnds[0] - p0)/s)
opt.set_upper_bounds((p_bnds[1] - p0)/s)

opt.add_inequality_constraint(self.constraint, 1E-8)

opt.set_initial_step(opt_options['initial_step'])
#opt.set_population(int(np.rint(10*(len(idx)+1)*10)))
opt.set_population(int(np.rint(10*(len(p0)+1)*10)))

#sub_opt = nlopt.opt(opt_options['algorithm'], len(p0))
#sub_opt.set_initial_step(opt_options['initial_step'])
#sub_opt.set_xtol_rel(opt_options['xtol_rel'])
#sub_opt.set_ftol_rel(opt_options['ftol_rel'])
#opt.set_local_optimizer(sub_opt)
sub_opt = nlopt.opt(opt_options['algorithm'], len(p0))
sub_opt.set_initial_step(opt_options['initial_step'])
sub_opt.set_xtol_rel(opt_options['xtol_rel'])
sub_opt.set_ftol_rel(opt_options['ftol_rel'])
opt.set_local_optimizer(sub_opt)

p0_opt = np.zeros_like(p0)
x_fit = opt.optimize(p0_opt) # optimize!
x_fit = x_fit*s + p0

print(x_fit)

# evaluate fits at new temperatures
[k_0, k_inf, Fcent], fit_arrhenius_coeffs = self.arrhenius_fcent_fit(x_fit, T_res, return_coeffs=True)
res = np.array([T_res, k_0, k_inf, Fcent]) # res = [T, k_0, k_inf, Fcent]
[k_0, k_inf, Fcent] = self.arrhenius_fcent_fit(x_fit, T_res)
res = np.array([T_res, k_0, k_inf, Fcent]).T # res = [T, k_0, k_inf, Fcent]

# change ln_A to A
x_fit[1] = np.exp(x_fit[1])
x_fit[4] = np.exp(x_fit[4])

# Fit HPL and LPL Arrhenius parameters
for idx, arrhenius_type in enumerate(['low_rate', 'high_rate']):
x_idx = np.array(self.alter_idx[arrhenius_type])
if idx in idx_fit:
self.x[x_idx] = fit_arrhenius_coeffs[x_idx]
self.x[x_idx] = x_fit[x_idx]
else:
self.x[x_idx] = arrhenius_coeffs[idx][x_idx]

return res

def arrhenius_fcent_fit(self, x, T_eval, return_coeffs=False):
def arrhenius_fcent_fit(self, x, T_eval):
T_arrhenius = self.obj_vars['T_arrhenius']
T_Fcent = self.obj_vars['T_Fcent']

#Fcent = x[2]
#k_0, k_inf = np.exp(x[:2])
k_0 = np.exp(x[:3])
k_inf = np.exp(x[3:6])
Ea_0, ln_A_0, n_0 = x[:3]
Ea_inf, ln_A_inf, n_inf = x[3:6]
ln_Fcent = x[6:]

Ea_0, A_0, n_0 = fit_arrhenius(k_0, T_arrhenius)
Ea_inf, A_inf, n_inf = fit_arrhenius(k_inf, T_arrhenius)
Fcent_poly = np.polynomial.Polynomial.fit(T_Fcent, ln_Fcent, 5) # Could later impose constraints on this fit
ln_k_0 = ln_A_0 + n_0*np.log(T_eval) - Ea_0/(Ru*T_eval)
ln_k_inf = ln_A_inf + n_inf*np.log(T_eval) - Ea_inf/(Ru*T_eval)

ln_k_0[ln_k_0 > max_ln_val] = max_ln_val
ln_k_inf[ln_k_inf > max_ln_val] = max_ln_val
k_0, k_inf = np.exp(ln_k_0), np.exp(ln_k_inf)

k_0 = A_0*np.power(T_eval, n_0)*np.exp(-Ea_0/Ru/T_eval)
k_inf = A_inf*np.power(T_eval, n_inf)*np.exp(-Ea_inf/Ru/T_eval)
Fcent_poly = np.polynomial.Polynomial.fit(T_Fcent, ln_Fcent, 5) # Could later impose constraints on this fit
Fcent = np.exp(np.polyval(Fcent_poly.convert().coef[::-1], T_eval))

if return_coeffs:
return [k_0, k_inf, Fcent], np.array([Ea_0, A_0, n_0, Ea_inf, A_inf, n_inf])
else:
return [k_0, k_inf, Fcent]
return [k_0, k_inf, Fcent]

def ln_Troe(self, T, *x): # LPL, HPL, Fcent
if len(x) == 1:
Expand Down Expand Up @@ -772,7 +780,7 @@ def calc_s(self, x, grad=[]):

s = 1/y
#s = s/np.min(s)
s = s/np.max(s)
#s = s/np.max(s)

return s

Expand All @@ -798,6 +806,26 @@ def objective(self, coefs_in, grad=np.array([]), obj_type='obj_sum'):

return obj_val

def constraint(self, coefs_in, grad=np.array([])):
coefs_0, s, T = self.obj_vars['coefs_0'], self.obj_vars['s'], self.obj_vars['T']
coefs = coefs_in*s + coefs_0
T_max = T[-1]

Ea_0, ln_A_0, n_0 = coefs[:3]
Ea_inf, ln_A_inf, n_inf = coefs[3:6]

ln_k_0 = ln_A_0 + n_0*np.log(T_max) - Ea_0/(Ru*T_max)
ln_k_inf = ln_A_inf + n_inf*np.log(T_max) - Ea_inf/(Ru*T_max)

ln_k_max = np.max([ln_k_0, ln_k_inf])

con = ln_k_max - max_ln_val

if grad.size > 0:
grad[:] = self.constraint_gradient(x, numerical_gradient=True)

return con


def fit_generic(rates, T, P, X, rxnIdx, coefKeys, coefNames, is_falloff_limit, mech, bnds, mpPool=None):
rxn = mech.gas.reaction(rxnIdx)
Expand Down
2 changes: 1 addition & 1 deletion src/calculate/optimize/mech_optimize.py
View file
Original file line number Diff line number Diff line change
Expand Up @@ -185,7 +185,7 @@ def _set_coef_opt(self):

return coef_opt

def _set_rxn_coef_opt(self, min_T_range=1500, min_P_range_factor=3):
def _set_rxn_coef_opt(self, min_T_range=800, min_P_range_factor=2):
mech = self.parent.mech
rxn_coef_opt = []
for coef in self.coef_opt:
Expand Down

0 comments on commit 0d89b68

Please sign in to comment.