"""Calculations for Example 11.5.6 from REB, The Book""" # import libraries import numpy as np import pandas as pd from reb_utils import lls_parameters from reb_utils import Arrhenius_parameters import matplotlib.pyplot as plt import os # global constants available to all functions # given V = 1.0 # L # known R = 8.314E-3 # kJ/mol/K # experimental data as arrays expt_df = pd.read_csv('example_11_5_6_data.csv') expt = expt_df['Experiment'].to_numpy() T = expt_df['T'].to_numpy() + 273.15 # K CA_0 = expt_df['CA0'].to_numpy() # mol /L t_meas = expt_df['t_meas'].to_numpy() # min CA_meas = expt_df['CA_meas'].to_numpy() # mol /L # deliverables function def deliverables(): # get the data block temperatures block_T_values = np.array(expt_df['T'].unique()) # allocate storage for the results am_k = np.ones_like(block_T_values)*float('nan') am_k_CI_lower = np.ones_like(block_T_values)*float('nan') am_k_CI_upper = np.ones_like(block_T_values)*float('nan') am_model_r_sq = np.ones_like(block_T_values)*float('nan') am_model_plot_data = [] # loop through the blocks for iBlock, blockT in enumerate(block_T_values): # extract the same-temperature data block block_df = expt_df[expt_df['T'] == blockT] # extract the data in the block expt_block = block_df['Experiment'].to_numpy() CA0_block = block_df['CA0'].to_numpy() t_meas_block = block_df['t_meas'].to_numpy() CA_meas_block = block_df['CA_meas'].to_numpy() # calculate x and add it to the data for this block x = -CA_meas_block am_model_plot_data.append(x) # calculate y using backward differences and save it y = np.ones_like(x)*float('nan') current_expt = -1 for iExpt, expt_number in enumerate(expt_block): # use the appropriate backward differences if expt_number != current_expt: # this is the first data point in the experiment y[iExpt] = (CA_meas_block[iExpt] \ - CA0_block[iExpt])/t_meas_block[iExpt] # reset the current experiment current_expt = expt_number else: # this is not the first data point in the experiment y[iExpt] = (CA_meas_block[iExpt] - CA_meas_block[iExpt-1])/(t_meas_block[iExpt] - t_meas_block[iExpt-1]) am_model_plot_data.append(y) # fit a straight line through the origin to the x-y data for this block param, param_ci, r_squared, y_pred = lls_parameters(y, x , model_has_intercept=False, use_rel_errors=False) # add y_pred to the data for this block am_model_plot_data.append(y_pred) # save the results for this block am_k[iBlock] = param[0] am_k_CI_lower[iBlock] = param_ci[0][0] am_k_CI_upper[iBlock] = param_ci[0][1] am_model_r_sq[iBlock] = r_squared # fit the Arrhenius expression to the T-k data T_K = block_T_values + 273.15 am_k0, am_k0_CI, am_E, am_E_CI, am_r_sq_Arr = Arrhenius_parameters(am_k, T_K , R) # calculate the predicted responses and residuals am_k_pred = am_k0*np.exp(-am_E/R/T_K) k = am_k0*np.exp(-am_E/R/T) am_CA_pred = CA_0*np.exp(-k*t_meas) am_epsilon_expt = CA_meas - am_CA_pred # calculate the overall r_squared CA_mean = np.mean(CA_meas) ss_res = np.sum(np.square(CA_meas - am_CA_pred)) ss_tot = np.sum(np.square(CA_meas - CA_mean)) am_r_sq = 1 - ss_res/ss_tot # make sure a results folder and a png folder exist if not os.path.isdir('results'): # create the folder os.makedirs('./results') if not os.path.isdir('png'): # create the folder os.makedirs('./png') # tabulate, show, and save the fitting function results results = [ ['k0', am_k0, am_k0_CI[0], am_k0_CI[1], '/min'] , ['E', am_E, am_E_CI[0], am_E_CI[1], 'kJ/mol'] , ['R^2 Arr', am_r_sq_Arr, float('nan'), float('nan'),''] , ['R^2 full', am_r_sq, float('nan'), float('nan'),''] ] results_df = pd.DataFrame(results) results_df.to_csv('results/am_results.csv', index=False) print('') print('Arrhenius fitting results:') print(results_df) print('') # tabulate and save the Arrhenius plot data results_df = pd.DataFrame({'block_T_values': block_T_values , 'am_k': am_k, 'am_k_pred': am_k_pred}) results_df.to_csv('results/am_Arr_plot_data.csv', index=False) # tabulate, show, and save the model plot results results_df = pd.DataFrame({'block_T_values': block_T_values , 'am_k': am_k, 'am_k_CI_lower': am_k_CI_lower , 'am_k_CI_upper': am_k_CI_upper, 'am_model_r_sq': am_model_r_sq}) results_df.to_csv('results/am_model_plot_results.csv', index=False) # tabulate and save the parity/residual plots data results_df = pd.DataFrame({'am_CA_pred': am_CA_pred , 'am_epsilon_expt': am_epsilon_expt}) results_df.to_csv('results/am_parity_plot_data.csv', index=False) # tabulate and save the model plot data plot_data = np.transpose(np.array(am_model_plot_data)) results_df = pd.DataFrame(plot_data) results_df.to_csv('results/am_model_plot_data.csv', index=False) # plot, show, and save the combined model plots plt.figure() plt.plot(plot_data[:,0], plot_data[:,1], marker='x', ls='', color='tab:blue' , label=f'{T_K[0] - 273.15:.0f} °C') plt.plot(plot_data[:,0], plot_data[:,2], ls='-', color='tab:blue') plt.plot(plot_data[:,3], plot_data[:,4], marker='x', ls='' , color='tab:orange', label=f'{T_K[1] - 273.15:.0f} °C') plt.plot(plot_data[:,3], plot_data[:,5], ls='-', color='tab:orange') plt.plot(plot_data[:,6], plot_data[:,7], marker='x', ls='' , color='tab:green', label=f'{T_K[2] - 273.15:.0f} °C') plt.plot(plot_data[:,6], plot_data[:,8], ls='-', color='tab:green') plt.plot(plot_data[:,9], plot_data[:,10], marker='x', ls='' , color='tab:purple', label=f'{T_K[3] - 273.15:.0f} °C') plt.plot(plot_data[:,9], plot_data[:,11], ls='-', color='tab:purple') plt.xlabel('x') plt.ylabel('y') plt.legend() plt.savefig('png/am_model_plots.png', dpi=300) plt.show(block=False) # plot, show, and save the Arrhenius plot plt.figure() plt.semilogy(1/T_K,am_k,color='tab:blue',marker='o', ls='none') plt.semilogy(1/T_K,am_k_pred,color='k') plt.xlabel('T$^{-1}$ (K$^{-1}$)') plt.ylabel('k (min$^{-1}$)') plt.xticks(rotation=25) plt.tight_layout() plt.savefig('png/am_Arrhenius_plot.png', dpi=300) plt.show(block=False) # generate, show, and save the overall parity plot plt.figure() plt.plot(CA_meas, am_CA_pred, marker='x', ls='', color='tab:blue' ,markerfacecolor='none', label='Data') plt.plot([min(CA_meas),max(CA_meas)],[min(CA_meas),max(CA_meas)] , color = 'k', ls = '-', label = 'Parity Line') plt.xlabel('$C_{A,meas}$ (M)') plt.ylabel('$C_{A,pred}$ (M)') plt.legend() plt.savefig('./png/am_parity.png', dpi=300) plt.show(block=False) # generate, show and save the residuals plots plt.figure() plt.plot(CA_0, am_epsilon_expt, markerfacecolor='none', marker='x' , color='tab:blue', ls='') plt.axhline(y=0, color='k') plt.xlabel('$C_{A,0}$ (M)') plt.ylabel('Residual (M)') plt.tight_layout() plt.savefig('./png/am_CA_residuals.png', dpi=300) plt.show(block=False) plt.figure() plt.plot(T-273.15, am_epsilon_expt, markerfacecolor='none', marker='x' , color='tab:blue', ls='') plt.axhline(y=0, color='k') plt.xlabel('T (°C)') plt.ylabel('Residual (M)') plt.tight_layout() plt.savefig('./png/am_T_residuals.png', dpi=300) plt.show(block=False) plt.figure() plt.plot(t_meas, am_epsilon_expt, markerfacecolor='none', marker='x' , color='tab:blue', ls='') plt.axhline(y=0, color='k') plt.xlabel('$t_{meas}$ (min)') plt.ylabel('Residual (M)') plt.tight_layout() plt.savefig('./png/am_time_residuals.png', dpi=300) plt.show() # execution command if __name__ == '__main__': deliverables()