diff --git a/src/calculate/optimize/fit_coeffs.py b/src/calculate/optimize/fit_coeffs.py index da6537e..c7265da 100644 --- a/src/calculate/optimize/fit_coeffs.py +++ b/src/calculate/optimize/fit_coeffs.py @@ -8,6 +8,7 @@ import warnings from copy import deepcopy from scipy.optimize import curve_fit, minimize, root_scalar, OptimizeWarning, least_squares, approx_fprime +from scipy import stats from timeit import default_timer as timer import itertools @@ -17,31 +18,23 @@ Ru = ct.gas_constant # Ru = 1.98720425864083 -min_pos_system_value = (np.finfo(float).tiny*(1E20))**(1/2) -max_pos_system_value = (np.finfo(float).max*(1E-20))**(1/2) +min_pos_system_value = (np.finfo(float).tiny*(1E20))**(0.5) +max_pos_system_value = (np.finfo(float).max*(1E-20))**(0.5) min_ln_val = np.log(min_pos_system_value) max_ln_val = np.log(max_pos_system_value) min_log_val = np.log10(min_pos_system_value) max_log_val = np.log10(max_pos_system_value) +ln_k_max = np.log(1E80) # max_log_val default_arrhenius_coefNames = ['activation_energy', 'pre_exponential_factor', 'temperature_exponent'] default_Troe_coefNames = ['activation_energy_0', 'pre_exponential_factor_0', 'temperature_exponent_0', 'activation_energy_inf', 'pre_exponential_factor_inf', 'temperature_exponent_inf', 'A', 'T3', 'T1', 'T2'] -#falloff_bnds = np.array([[1E-12, 1E-12, 0.1], [1E30, 1E30, 1]]) -#falloff_bnds[:,0:2] = np.log(falloff_bnds[:,0:2]) # ln(k_0), ln(k_inf), Fcent -falloff_bnds = np.log([[1E-6, 1E-6, 0.1], [1E30, 1E30, 1]]) # ln(k_0), ln(k_inf), ln(Fcent) - -troe_falloff_0 = [[0.6, 200, 600, 1200], # (0, 0, 0) - # (0, 0, 1) - [0.05, 1000, -2000, 3000], # (0, 1, 0) - [-0.3, 200, -20000, -50], # (0, 1, 1) - [0.9, -2000, 500, 10000], # (1, 0, 0) - # (1, 0, 1) - # (1, 1, 0) - # (1, 1, 1) - ] +troe_falloff_0 = [[np.nan, np.nan, np.nan, 1500], # only T2 matters + [0.6, 200, 600, 1200], # (0, 0, 0) + [0.05, 1000, -2000, 3000], # (0, 1, 0) + [0.9, -2000, 500, 10000]] # (1, 0, 0) #troe_min_pos_system_value = (np.finfo(float).tiny*(1E20))**(1/3) #troe_max_pos_system_value = (np.finfo(float).max*(1E-20))**(1/3) @@ -57,10 +50,15 @@ # '+': [-T_max/troe_min_ln_val, troe_max_pos_system_value]}, # 'T2': {'-': [-T_max*troe_max_ln_val, -T_max*troe_min_ln_val], # '+': [-T_max*troe_max_ln_val, -T_max*troe_min_ln_val]}} -troe_all_bnds = {'A': {'-': [-1E2, 1.0], '+': [-1E2, 1.0]}, - 'T3': {'-': [-1E30, -30], '+': [30.0, 1E30]}, - 'T1': {'-': [-1E30, -30], '+': [30.0, 1E30]}, - 'T2': {'-': [-1E4, 1E30], '+': [-1E4, 1E30]}} +#troe_all_bnds = {'A': {'-': [-1E2, 1.0], '+': [-1E2, 1.0]}, +# 'T3': {'-': [-1E14, -50], '+': [50.0, 1E14]}, +# 'T1': {'-': [-1E14, -50], '+': [50.0, 1E14]}, +# 'T2': {'-': [-1E4, 1E14], '+': [-1E4, 1E14]}} + +troe_all_bnds = {'A': {'-': [-1E2, 1.0], '+': [-1E2, 1.0]}, + 'T3': {'-': [-1E8, -1E2], '+': [ 1.0, 1E9]}, + 'T1': {'-': [-1E8, -1E2], '+': [ 1.0, 1E9]}, + 'T2': {'-': [-1E4, 1E6], '+': [-1E4, 1E6]}} def fit_arrhenius(rates, T, x0=[], coefNames=default_arrhenius_coefNames, bnds=[], loss='linear'): def fit_fcn_decorator(x0, alter_idx, jac=False): @@ -135,16 +133,22 @@ def ln_arrhenius_jac(T, *args): bisymlog_C = 1/(np.exp(1)-1) class falloff_parameters: # based on ln_Fcent - def __init__(self, T, Fcent_target, x0=[0.6, 200, 600, 1200], use_scipy=False, + def __init__(self, T, M, ln_k, x0=[0.6, 200, 600, 1200], use_scipy=False, nlopt_algo=nlopt.LN_SBPLX, loss_fcn_par=[2, 1]): self.T = T - self.Fcent = Fcent_target + self.M = M + self.ln_k = ln_k self.x0 = x0 self.s = np.ones_like(x0) - self.Fcent_min = 1.1E-8 - self.Tmin = 100 - self.Tmax = 20000 + self.Fcent_min = 1E-6 + self.Tmin = 273 # 100 + self.Tmax = 6000 # 20000 + + if np.isnan(self.x0).any(): + self.T2_only = True + else: + self.T2_only = False self.use_scipy = use_scipy @@ -157,24 +161,45 @@ def x_bnds(self, x0): for n, coef in enumerate(['A', 'T3', 'T1', 'T2']): if x0[n] < 0: bnds.append(troe_all_bnds[coef]['-']) - else: + elif x0[n] > 0: bnds.append(troe_all_bnds[coef]['+']) + else: # if doesn't match either of the above then put nan as bounds + bnds.append([np.nan, np.nan]) return np.array(bnds).T def fit(self): - Fcent = self.Fcent T = self.T - x0 = self.x0 - p0 = self.convert(x0, 'base2opt') - self.p_bnds = self.convert(self.x_bnds(x0)) + + self.p0 = self.x0 + self.p0[-4:] = self.convert_Fcent(self.p0[-4:], 'base2opt') + + #s = np.array([4184, 1.0, 1E-2, 4184, 1.0, 1E-2, *(np.ones((1,4)).flatten()*1E-1)]) + + p_bnds = set_arrhenius_bnds(self.p0[0:3], default_arrhenius_coefNames) + p_bnds = np.concatenate((p_bnds, set_arrhenius_bnds(self.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 + + Fcent_bnds = self.convert_Fcent(self.x_bnds(self.p0[-4:])) + + self.p_bnds = np.concatenate((p_bnds, Fcent_bnds), axis=1) + + if len(self.p_bnds) > 0: + self.p0 = np.clip(self.p0, self.p_bnds[0, :], self.p_bnds[1, :]) + + self.p0 = self.p0[~np.isnan(self.p0)] + self.p_bnds = self.p_bnds[~np.isnan(self.p_bnds)] + self.p_bnds = np.reshape(self.p_bnds, (2,-1)) + + self.s = np.ones_like(self.p0) if self.use_scipy: with warnings.catch_warnings(): warnings.simplefilter('ignore', OptimizeWarning) - x_fit, _ = curve_fit(self.function, T, Fcent, p0=p0, method='trf', bounds=self.p_bnds, # dogbox - #jac=fit_func_jac, x_scale='jac', max_nfev=len(p0)*1000) - jac='2-point', x_scale='jac', max_nfev=len(p0)*1000, loss='huber') + x_fit, _ = curve_fit(self.ln_Troe, T, self.ln_k, p0=self.p0, method='trf', bounds=self.p_bnds, # dogbox + #jac=fit_func_jac, x_scale='jac', max_nfev=len(self.p0)*1000) + jac='2-point', x_scale='jac', max_nfev=len(self.p0)*1000, loss='huber') #print('scipy:', x_fit) #cmp = np.array([T, Fcent, np.exp(fit_func(T, *x_fit))]).T @@ -186,20 +211,20 @@ def fit(self): else: xtol_rel = 1E-8 ftol_rel = 1E-4 - initial_step = 1E-4 + initial_step = 1E-2 - self.opt = opt = nlopt.opt(nlopt.AUGLAG, 4) + self.opt = opt = nlopt.opt(nlopt.AUGLAG, len(self.p0)) opt.set_min_objective(self.objective) - opt.add_inequality_constraint(self.constraint, 1E-8) + opt.add_inequality_constraint(self.Fcent_constraint, 0.0) + opt.add_inequality_constraint(self.Arrhenius_constraint, 1E-8) opt.set_maxeval(10000) - opt.set_maxtime(10) + #opt.set_maxtime(10) opt.set_xtol_rel(xtol_rel) opt.set_ftol_rel(ftol_rel) - self.p0 = p0 - p0_opt = np.zeros_like(p0) + p0_opt = np.zeros_like(self.p0) self.s = self.calc_s(p0_opt) opt.set_lower_bounds((self.p_bnds[0]-self.p0)/self.s) @@ -208,14 +233,16 @@ def fit(self): opt.set_initial_step(initial_step) #opt.set_population(int(np.rint(10*(len(idx)+1)*10))) - sub_opt = nlopt.opt(self.opt_algo, 4) + sub_opt = nlopt.opt(self.opt_algo, len(self.p0)) sub_opt.set_initial_step(initial_step) sub_opt.set_xtol_rel(xtol_rel) sub_opt.set_ftol_rel(ftol_rel) opt.set_local_optimizer(sub_opt) x_fit = opt.optimize(p0_opt) # optimize! - x_fit = x_fit*self.s + self.p0 + #print('Fcent_constraint: ', self.Fcent_constraint(x_fit)) + #print('Arrhe_constraint: ', self.Arrhenius_constraint(x_fit)) + x_fit = self.set_x_from_opt(x_fit) #print('nlopt:', x_fit) #cmp = np.array([T, Fcent, self.function(T, *x_fit)]).T @@ -224,12 +251,34 @@ def fit(self): # print(*entry) #print('') - x = self.convert(x_fit, 'opt2base') - res = {'x': x, 'fval': opt.last_optimum_value(), 'nfev': opt.get_numevals()} + # change ln_A to A # TODO: Maybe move this to convert + x_fit[1] = np.exp(x_fit[1]) + x_fit[4] = np.exp(x_fit[4]) + + res = {'x': x_fit, 'fval': opt.last_optimum_value(), 'nfev': opt.get_numevals()} return res - def convert(self, x, conv_type='base2opt'): + def set_x_from_opt(self, x): + x = x*self.s + self.p0 + + if self.T2_only: + x = np.array([*x[:-1], 1.0, 1.0E-30, 1.0E-30, x[-1]]) + x[-1] = self.convert_Fcent([x[-1]], 'opt2base')[0] # A, T3, T1, T2 + + else: + x[-4:] = self.convert_Fcent(x[-4:], 'opt2base') # A, T3, T1, T2 + + for i in [-3, -2]: + if x[i] >= 0: + if x[i] < 10: + x[i] = 1E-30 + elif x[i] > 1E8: + x[i] = 1E30 + + return x + + def convert_Fcent(self, x, conv_type='base2opt'): #x = x*self.s + self.p0 y = np.array(x) C = bisymlog_C @@ -251,36 +300,81 @@ def convert(self, x, conv_type='base2opt'): if flatten: # if it's 1d it's the guess y = y.flatten() + else: # if it 2d it's the bnds, they need to be sorted y = np.sort(y, axis=0) return y - def function(self, T, *x): - [A, T3, T1, T2] = self.convert(x, 'opt2base') + def ln_Troe(self, T, *x): # LPL, HPL, Fcent + M = self.M + + Ea_0, ln_A_0, n_0 = x[:3] + Ea_inf, ln_A_inf, n_inf = x[3:6] + Fcent_coeffs = x[-4:] # A, T3, T1, T2 + + ln_k_0 = ln_A_0 + n_0*np.log(T) - Ea_0/(Ru*T) + ln_k_inf = ln_A_inf + n_inf*np.log(T) - Ea_inf/(Ru*T) + + 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) + + Fcent = self.Fcent_calc(T[0,:], *Fcent_coeffs) + Fcent[Fcent <= 0.0] = 10000 + + with np.errstate(all='raise'): + try: + P_r = k_0/k_inf*M + log_P_r = np.log10(P_r) - Fcent_fit = self.Fcent_calc(T, A, T3, T1, T2) - Fcent_fit[Fcent_fit <= 0.0] = 10000 + #log_P_r = np.log10(k_0) - np.log10(k_inf) + np.log10(M) + #log_P_r[log_P_r < min_log_val] = min_log_val + #log_P_r[log_P_r > max_log_val] = max_log_val + #P_r = np.power(10, log_P_r) + except: + return np.ones_like(M)*max_pos_system_value + + log_Fcent = np.log10(Fcent) + C = -0.4 - 0.67*log_Fcent + N = 0.75 - 1.27*log_Fcent + f1 = (log_P_r + C)/(N - 0.14*(log_P_r + C)) + + ln_F = np.log(Fcent)/(1 + f1**2) + with np.errstate(all='raise'): + try: + ln_k_calc = np.log(k_inf*P_r/(1 + P_r)) + ln_F + except: + return np.ones_like(M)*max_pos_system_value - return Fcent_fit + return ln_k_calc def Fcent_calc(self, T, A, T3, T1, T2): - exp_T3 = np.zeros_like(T) - exp_T3[:] = np.exp(-T/T3) + def exp_safe(num, den): + with np.errstate(divide='ignore'): + x = num/den - exp_T1 = np.zeros_like(T) - exp_T1[:] = np.exp(-T/T1) + res = np.zeros_like(x) # since initialized to zero, values below are already zero - T2_T = T2/T - exp_T2 = np.where(T2_T <= max_ln_val, np.exp(-T2_T), 0.0) + idx_valid = np.argwhere((x >= min_ln_val) & (x <= max_ln_val)) + res[idx_valid] = np.exp(x[idx_valid]) + + idx_over = np.argwhere(x > max_ln_val) + res[idx_over] = max_ln_val + + return res + + exp_T3 = exp_safe(-T, T3) + exp_T1 = exp_safe(-T, T1) + exp_T2 = exp_safe(-T2, T) Fcent = (1-A)*exp_T3 + A*exp_T1 + exp_T2 return Fcent - def jacobian(self, T, *x): + def jacobian(self, T, *x): # defunct [A, B, C, D] = x - [A, T3, T1, T2] = self.convert(x, 'opt2base') # A, B, C, D = x + [A, T3, T1, T2] = self.convert_Fcent(x, 'opt2base') # A, B, C, D = x bC = bisymlog_C jac = [] @@ -294,10 +388,10 @@ def jacobian(self, T, *x): return jac def objective(self, x_fit, grad=np.array([]), obj_type='obj_sum'): - x = x_fit*self.s + self.p0 + x = self.set_x_from_opt(x_fit) T = self.T - resid = self.function(T, *x) - self.Fcent + resid = self.ln_Troe(T, *x) - self.ln_k if obj_type == 'obj_sum': obj_val = generalized_loss_fcn(resid, a=self.loss_alpha, c=self.loss_scale).sum() elif obj_type == 'obj': @@ -307,7 +401,7 @@ def objective(self, x_fit, grad=np.array([]), obj_type='obj_sum'): #s[:] = np.abs(np.sum(loss*fit_func_jac(T, *x).T, axis=1)) if grad.size > 0: - grad[:] = self.objective_gradient(x, resid) + grad[:] = self.objective_gradient(x, resid, numerical_gradient=True) #else: # grad = self.objective_gradient(x, resid) @@ -341,7 +435,7 @@ def objective_gradient(self, x, resid=[], numerical_gradient=False): def calc_s(self, x, grad=[]): if len(grad) == 0: - grad = self.objective_gradient(x) + grad = self.objective_gradient(x, numerical_gradient=True) y = np.abs(grad) if (y < min_pos_system_value).all(): @@ -351,11 +445,11 @@ 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 - def constraint(self, x, grad=np.array([])): + def Fcent_constraint(self, x_fit, grad=np.array([])): def f_fp(T, A, T3, T1, T2, fprime=False, fprime2=False): # dFcent_dT f = T2/T**2*np.exp(-T2/T) - (1-A)/T3*np.exp(-T/T3) - A/T1*np.exp(-T/T1) @@ -365,9 +459,8 @@ def f_fp(T, A, T3, T1, T2, fprime=False, fprime2=False): # dFcent_dT fp = T2*(T2 - 2*T)/T**4*np.exp(-T2/T) + (1-A)/T3**2*np.exp(-T/T3) +A/T1**2*np.exp(-T/T1) return f, fp - x = x*self.s + self.p0 - [A, T3, T1, T2] = self.convert(x, 'opt2base') - Fcent_min = self.Fcent_min + x = self.set_x_from_opt(x_fit) + [A, T3, T1, T2] = x[-4:] Tmin = self.Tmin Tmax = self.Tmax @@ -382,16 +475,39 @@ def f_fp(T, A, T3, T1, T2, fprime=False, fprime2=False): # dFcent_dT print(T) Fcent = self.Fcent_calc(T, A, T3, T1, T2) #TODO: OVERFLOW WARNING HERE - con = np.min(Fcent_min - Fcent) + min_con = np.max(self.Fcent_min - Fcent) + max_con = np.max(Fcent - 1.0) + con = np.max([max_con, min_con])*1E8 if grad.size > 0: - grad[:] = self.constraint_gradient(x, numerical_gradient=True) + grad[:] = self.constraint_gradient(x, numerical_gradient=self.Fcent_constraint) return con - def constraint_gradient(self, x, const_eval=[], numerical_gradient=False): - if numerical_gradient: - grad = approx_fprime(x, self.constraint, 1E-10) + def Arrhenius_constraint(self, x_fit, grad=np.array([])): + x = self.set_x_from_opt(x_fit) + + T = self.T + T_max = T[-1] + + Ea_0, ln_A_0, n_0 = x[:3] + Ea_inf, ln_A_inf, n_inf = x[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_limit_max = np.max([ln_k_0, ln_k_inf]) + + con = ln_k_limit_max - ln_k_max + + if grad.size > 0: + grad[:] = self.constraint_gradient(x, numerical_gradient=self.Arrhenius_constraint) + + return con + + def constraint_gradient(self, x, const_eval=[], numerical_gradient=None): + if numerical_gradient is not None: + grad = approx_fprime(x, numerical_gradient, 1E-10) else: # I've not calculated the derivatives wrt coefficients for analytical if len(resid) == 0: @@ -410,14 +526,15 @@ def constraint_gradient(self, x, const_eval=[], numerical_gradient=False): def falloff_parameters_decorator(args_list): - T_falloff, Fcent, x0, use_scipy, nlopt_algo = args_list - falloff = falloff_parameters(T_falloff, Fcent, x0=x0, use_scipy=use_scipy, nlopt_algo=nlopt_algo) #GN_CRS2_LM LN_SBPLX + T, M, ln_k, x0, use_scipy, nlopt_algo = args_list + falloff = falloff_parameters(T, M, ln_k, x0=x0, use_scipy=use_scipy, nlopt_algo=nlopt_algo) return falloff.fit() +bisymlog_C = 1/(np.exp(1)-1) class Troe: def __init__(self, rates, T, P, M, x0=[], coefNames=default_Troe_coefNames, bnds=[], is_falloff_limit=None, - scipy_curvefit=False, robust=False, mpPool=None, nlopt_algo=nlopt.LN_SBPLX): + scipy_curvefit=False, robust=False, mpPool=None, nlopt_algo=nlopt.GN_CRS2_LM): # GN_DIRECT_L, GN_DIRECT, GN_DIRECT_NOSCAL GN_CRS2_LM, LN_COBYLA, LN_SBPLX, LD_MMA self.debug = True @@ -427,7 +544,6 @@ def __init__(self, rates, T, P, M, x0=[], coefNames=default_Troe_coefNames, bnds self.P = P self.M = M self.P_limit_constrained = [False, False] # is low pressure limit, high pressure limit - self.ln_Fcent_min = np.log(1E-3) self.x0 = x0 if len(self.x0) != 10: @@ -450,7 +566,7 @@ def __init__(self, rates, T, P, M, x0=[], coefNames=default_Troe_coefNames, bnds self.scipy_curvefit = scipy_curvefit self.is_falloff_limit = is_falloff_limit - self.nlopt_algo = nlopt_algo + self.nlopt_algo = nlopt_algo self.pool = mpPool if robust: @@ -469,14 +585,58 @@ def __init__(self, rates, T, P, M, x0=[], coefNames=default_Troe_coefNames, bnds else: self.alter_idx['falloff_parameters'].append(n) + def convert_Fcent(self, x, conv_type='base2opt'): + #x = x*self.s + self.p0 + y = np.array(x) + C = bisymlog_C + + flatten = False + if y.ndim == 1: + y = y[np.newaxis, :] + flatten = True + + if conv_type == 'base2opt': # y = [A, T3, T1, T2] + y[:,1:] = np.sign(y[:,1:])*np.log(np.abs(y[:,1:]/C) + 1) + + else: + y[:,1:] = np.sign(y[:,1:])*C*(np.exp(np.abs(y[:,1:])) - 1) + + #A = np.log(y[0]) + #T3, T1 = 1000/y[1], 1000/y[2] + #T2 = y[3]*100 + + if flatten: # if it's 1d it's the guess + y = y.flatten() + else: # if it 2d it's the bnds, they need to be sorted + y = np.sort(y, axis=0) + + return y + def fit(self): x = self.x start = timer() - # fit LPL, HPL, Fcent - res = self.LPL_HPL_Fcent() + arrhenius_coeffs, T, ln_k = self.fit_arrhenius_rates() # fit arrhenius expression for each P - # fit Troe falloff parameters + # set T, P, M as arrays + P = np.unique(self.P) + P[1:] = np.roll(P[1:], 1) + T, P = np.meshgrid(T, P) + + M = [] # calculate M + for i in range(0, np.shape(P)[0]): + M.append(self.M(T[i,:], P[i,:])) + + M = np.array(M) + + # calculate k_0, k_inf, Fcent + # set indices to fit [log_k_0, log_k_inf, log_Fcent] + idx_fit = [0, 1, 2] + for i in range(0, 2)[::-1]: + if self.P_limit_constrained[i]: + del idx_fit[i] + + # fit Troe parameters idx = self.alter_idx['falloff_parameters'] #if len(self.x0) > 0 and list(self.x0[idx]) not in troe_falloff_0: # x0 = np.array(troe_falloff_0) @@ -484,17 +644,16 @@ def fit(self): #else: # x0 = troe_falloff_0 x0 = troe_falloff_0 - - T_falloff, Fcent = res[:,0], res[:,3] if self.pool is not None: - args_list = ((T_falloff, Fcent, Troe_0, False, self.nlopt_algo) for Troe_0 in x0) + args_list = ((T, M, ln_k, np.array([3.6E7, 137.0, -7.0, 1.4E7, 34.5, 0.4, *Troe_0]), False, self.nlopt_algo) for Troe_0 in x0) falloff_output = self.pool.map(falloff_parameters_decorator, args_list) else: falloff_output = [] for i, Troe_0 in enumerate(x0): - falloff = falloff_parameters(T_falloff, Fcent, x0=Troe_0, use_scipy=False, - nlopt_algo=nlopt.LN_SBPLX, loss_fcn_par=self.loss_fcn_par) #GN_CRS2_LM LN_SBPLX + p0 = np.array([3.6E7, 137.0, -7.0, 1.4E7, 34.5, 0.4, *Troe_0]) + falloff = falloff_parameters(T, M, ln_k, x0=p0, use_scipy=False, + nlopt_algo=self.nlopt_algo, loss_fcn_par=self.loss_fcn_par) #GN_CRS2_LM LN_SBPLX falloff_output.append(falloff.fit()) HoF = {'obj_fcn': np.inf, 'coeffs': []} @@ -504,7 +663,28 @@ def fit(self): HoF['coeffs'] = res['x'] HoF['i'] = i - x[idx] = HoF['coeffs'] + # Run local optimizer + x = HoF['coeffs'] + + str_x = '\t'.join([f'{val:.3e}' for val in x]) + print(f'G {len(x)}\t{HoF["obj_fcn"]:.3e}\t', str_x) + + x0 = deepcopy(x) + x0[1], x0[4] = np.log(x0[1]), np.log(x0[4]) + + if np.all(x0[-3:-1] == [1.0E-30, 1.0E-30]): + x0[-4:-1] = [np.nan, np.nan, np.nan] + + falloff = falloff_parameters(T, M, ln_k, x0=x0, use_scipy=False, + nlopt_algo=nlopt.LN_SBPLX, loss_fcn_par=self.loss_fcn_par) + res = falloff.fit() + + if res['fval'] < HoF['obj_fcn'] and len(res['x']) == len(x): # this length check shouldn't be needed + x = res['x'] + + #print(f'{HoF["obj_fcn"]:.3e}\t{res["fval"]:.3e}') + str_x = '\t'.join([f'{val:.3e}' for val in x]) + print(f'L {len(x)}\t{res["fval"]:.3e}\t', str_x) #if self.debug: # T = self.T @@ -568,7 +748,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_CRS2_LM, '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 LN_COBYLA '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}, @@ -582,41 +762,7 @@ def LPL_HPL_Fcent(self): # CURRENTLY A KNOWN BUG IF FULL LIMIT IS HELD CONSTA else: loss = 'linear' - # set T, P, M as arrays - P = np.unique(self.P) - P[1:] = np.roll(P[1:], 1) - T_max = np.max(T) - T, P = np.meshgrid(T, P) - - T_arrhenius = np.linspace(T[0,0], T[0,-1], 3) # for fitting HPL/LPL in fit fcn - T_Fcent = np.linspace(T[0,0], T[0,-1], 6) - T_res = np.linspace(T[0,0], T[0,-1], 50) - - M = [] # calculate M - for i in range(0, np.shape(P)[0]): - M.append(self.M(T[i,:], P[i,:])) - - M = np.array(M) - # set default p0, s, and p_bnds - 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 = 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 - # set indices to fit [log_k_0, log_k_inf, log_Fcent] - idx_fit = [0, 1, 2] - for i in range(0, 2)[::-1]: - if self.P_limit_constrained[i]: - del idx_fit[i] res = [] for n, opt_options in enumerate(options): @@ -690,7 +836,8 @@ def LPL_HPL_Fcent(self): # CURRENTLY A KNOWN BUG IF FULL LIMIT IS HELD CONSTA 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.add_inequality_constraint(self.Arrhenius_constraint, 1E-8) + opt.add_inequality_constraint(self.Fcent_constraint, 1E-8) opt.set_initial_step(opt_options['initial_step']) opt.set_population(int(np.rint(10*(len(p0)+1)*10))) @@ -723,120 +870,6 @@ def LPL_HPL_Fcent(self): # CURRENTLY A KNOWN BUG IF FULL LIMIT IS HELD CONSTA return res - 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]) - Ea_0, ln_A_0, n_0 = x[:3] - Ea_inf, ln_A_inf, n_inf = x[3:6] - ln_Fcent = x[6:] - - 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) - - 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)) - - return [k_0, k_inf, Fcent] - - def ln_Troe(self, T, *x): # LPL, HPL, Fcent - if len(x) == 1: - x = x[0] - - M = self.obj_vars['M'] - - [k_0, k_inf, Fcent] = self.arrhenius_fcent_fit(x, T) - - with np.errstate(all='raise'): - try: - P_r = k_0/k_inf*M - log_P_r = np.log10(P_r) - - #log_P_r = np.log10(k_0) - np.log10(k_inf) + np.log10(M) - #log_P_r[log_P_r < min_log_val] = min_log_val - #log_P_r[log_P_r > max_log_val] = max_log_val - #P_r = np.power(10, log_P_r) - except: - return np.ones_like(M)*max_pos_system_value - - log_Fcent = np.log10(Fcent) - C = -0.4 - 0.67*log_Fcent - N = 0.75 - 1.27*log_Fcent - f1 = (log_P_r + C)/(N - 0.14*(log_P_r + C)) - - ln_F = np.log(Fcent)/(1 + f1**2) - with np.errstate(all='raise'): - try: - ln_k_calc = np.log(k_inf*P_r/(1 + P_r)) + ln_F - except: - return np.ones_like(M)*max_pos_system_value - - return ln_k_calc - - def calc_s(self, x, grad=[]): - if len(grad) == 0: - grad = approx_fprime(x, self.objective, 1E-10) - - y = np.abs(grad) - if (y < min_pos_system_value).all(): - y = np.ones_like(y)*1E-14 - else: - y[y < min_pos_system_value] = 10**(OoM(np.min(y[y>=min_pos_system_value])) - 1) # TODO: MAKE THIS BETTER running into problem when s is zero, this is a janky workaround - - s = 1/y - #s = s/np.min(s) - #s = s/np.max(s) - - return s - - def objective(self, coefs_in, grad=np.array([]), obj_type='obj_sum'): - fit_fcn, T, ln_k = self.obj_vars['fcn'], self.obj_vars['T'], self.obj_vars['ln_k'] - coefs_0, s = self.obj_vars['coefs_0'], self.obj_vars['s'] - loss_scale, loss_alpha = self.loss_fcn_par - - coefs = coefs_in*s + coefs_0 - - resid = fit_fcn(T, coefs) - ln_k - if obj_type == 'obj_sum': - obj_val = generalized_loss_fcn(resid, a=loss_alpha, c=loss_scale).sum() - elif obj_type == 'obj': - obj_val = generalized_loss_fcn(resid, a=loss_alpha, c=loss_scale) - elif obj_type == 'resid': - obj_val = resid - - if grad.size > 0: - grad[:] = self.objective_gradient(coefs, resid) - #else: - # grad = self.objective_gradient(coefs, resid) - - 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/misc_fcns.py b/src/calculate/optimize/misc_fcns.py index 365c4c2..2a431f1 100644 --- a/src/calculate/optimize/misc_fcns.py +++ b/src/calculate/optimize/misc_fcns.py @@ -8,8 +8,8 @@ Ru = ct.gas_constant -min_pos_system_value = np.finfo(float).tiny*1E20 -max_pos_system_value = np.finfo(float).max*1E-20 +min_pos_system_value = (np.finfo(float).tiny*(1E20))**(0.5) +max_pos_system_value = (np.finfo(float).max*(1E-20))**(0.5) min_neg_system_value = -max_pos_system_value T_min = 300 T_max = 6000 diff --git a/src/calculate/reactors.py b/src/calculate/reactors.py index 785869c..994e197 100644 --- a/src/calculate/reactors.py +++ b/src/calculate/reactors.py @@ -220,10 +220,20 @@ def incident_shock_reactor(self, gas, details, t_end, **kwargs): y0 = np.hstack((0.0, var['A1'], gas.density, var['u_reac'], gas.T, 0.0, gas.Y)) # Initial condition ode = shock_fcns.ReactorOde(gas, t_end, var['rho1'], var['L'], var['As'], var['A1'], False) - sol = integrate.solve_ivp(ode, [0, t_end], y0, method=var['ODE_solver'], - dense_output=True, rtol=var['rtol'], atol=var['atol']) + with np.errstate(over='raise', divide='raise'): + try: + sol = integrate.solve_ivp(ode, [0, t_end], y0, method=var['ODE_solver'], + dense_output=True, rtol=var['rtol'], atol=var['atol']) + sol_success = True + sol_message = sol.message + sol_t = sol.t + + except: + sol_success = False + sol_message = sys.exc_info()[0] + sol_t = sol.t - if sol.success: + if sol_success: self.ODE_success = True # this is passed to SIM to inform saving output function details['success'] = True else: @@ -235,18 +245,18 @@ def incident_shock_reactor(self, gas, details, t_end, **kwargs): checkRxns = self.checkRxnRates(gas) if len(checkRxns) > 0: explanation += '\nSuggested Reactions: ' + ', '.join([str(x) for x in checkRxns]) - details['message'] = '\nODE Error: {:s}\n{:s}\n'.format(sol.message, explanation) + details['message'] = '\nODE Error: {:s}\n{:s}\n'.format(sol_message, explanation) - if var['sim_int_f'] > np.shape(sol.t)[0]: # in case of integration failure - var['sim_int_f'] = np.shape(sol.t)[0] + if var['sim_int_f'] > np.shape(sol_t)[0]: # in case of integration failure + var['sim_int_f'] = np.shape(sol_t)[0] if var['sim_int_f'] == 1: - t_sim = sol.t + t_sim = sol_t else: # perform interpolation if integrator sample factor > 1 j = 0 - t_sim = np.zeros(var['sim_int_f']*(np.shape(sol.t)[0] - 1) + 1) # preallocate array - for i in range(np.shape(sol.t)[0]-1): - t_interp = np.interp(np.linspace(i, i+1, var['sim_int_f']+1), [i, i+1], sol.t[i:i+2]) + t_sim = np.zeros(var['sim_int_f']*(np.shape(sol_t)[0] - 1) + 1) # preallocate array + for i in range(np.shape(sol_t)[0]-1): + t_interp = np.interp(np.linspace(i, i+1, var['sim_int_f']+1), [i, i+1], sol_t[i:i+2]) t_sim[j:j+len(t_interp)] = t_interp j += len(t_interp) - 1