diff --git a/src/calculate/optimize/fit_coeffs.py b/src/calculate/optimize/fit_coeffs.py index 7f5d914..b3a2711 100644 --- a/src/calculate/optimize/fit_coeffs.py +++ b/src/calculate/optimize/fit_coeffs.py @@ -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}, @@ -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 @@ -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']) @@ -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: @@ -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 @@ -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) diff --git a/src/calculate/optimize/mech_optimize.py b/src/calculate/optimize/mech_optimize.py index 82d77b2..89ce9bb 100644 --- a/src/calculate/optimize/mech_optimize.py +++ b/src/calculate/optimize/mech_optimize.py @@ -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: