forked from jtr5395/regression
-
Notifications
You must be signed in to change notification settings - Fork 0
/
FittingExercise.py
154 lines (124 loc) · 5.71 KB
/
FittingExercise.py
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
# coding: utf-8
# # Fitting of experimental chemical kinetics data
# You perform some experiments in a batch reactor to determine the rate expression and thermochemistry for the reversible chemical reaction
# $\require{mhchem}$
# $$\ce{A <=> B}$$
#
# Recall from thermodynamics that
# $\Delta G = \Delta H - T \Delta S$
# and $\Delta G = R T \ln K_a$
# where $K_a$ is the activity-based equilibrium constant of the chemical reaction, $R$ is the molar gas constant (8.314 J/mol/K) and $T$ is the temperature in Kelvin.
#
# If we assume ideal solution (unit fugacity coefficients) then $K_a = K_c$ giving us the concentration-based equilibrium constant $K_c$.
# From kinetics you recall
# $K_c = \frac{k_f}{k_r}$
# where
# $k_f$ is the forward rate coefficient and $k_r$ is the reverse rate coefficient.
# i.e. the rate of the reaction $\ce{A->B}$ is $k_f \times C_A$
# and the reverse reaction $\ce{B->A}$ is $k_r \times C_B$
# where $C_A$ and $C_B$ are the concentrations of species A and B respectively.
# In a batch reactor $\frac{dN_A}{dt} = r_{A(net)} V$, so (dividing through by the reactor volume $V$) $\frac{dC_A}{dt} = r_{A(net)}$ where $r_{A(net)}$ is the net rate of formation of species A, i.e. $r_{A(net)} = k_r C_B - k_f C_A$.
# Assume the forward rate coefficient $k_f$ follows Arrhenius form, $k_f = A \exp\left(\frac{-E_A}{R T}\right)$ where $A$ is the "pre-exponential factor" and $E_A$ is the activation energy.
#
# Fortunately, in this case you have good reason to believe that species A and B have very similar temperature-dependent heat capacities, so that $\Delta H_{rxn}$ and $\Delta S_{rxn}$ are independent of temperature.
#
# You start the experiment with no B ($C_B=0$), and at time zero have some way to initiate the reaction, starting with a set concentration of $C_A$.
#
# You wish to determine the four paramaters:
# $log_{10} A$,
# $E_A$,
# $\Delta H_{rxn}$,
# $\Delta S_{rxn}$.
#
# Based on a literature search, quantum chemistry calculations, and prior experience, your current estimates are as follows:
# ```
# logA = 6. # base-ten logarithm of A in s^-1
# Ea = 45. # Ea in kJ/mol
# dH = -10. # ∆H in kJ/mol
# dS = -50. # ∆S in J/mol/K
# ```
#
# In[1]:
get_ipython().magic('matplotlib inline')
import numpy as np
import scipy.integrate
from matplotlib import pyplot as plt
import random
import SALib as sa
import SALib.sample
# from SALib.sample import morris as ms
# from SALib.analyze import morris as ma
# from SALib.plotting import morris as mp
# In[2]:
# This cell just tries to make graphs look nicer
try:
import seaborn as sns
except ImportError:
# This block will be run if there's an ImportError, i.e you don't have seaborn installed.
sns = False
print ("If you want to try different figure formatting, "
"type 'conda install seaborn' at an anaconda command prompt or terminal. "
"See https://stanford.edu/~mwaskom/software/seaborn/ for details")
# If not using seaborn, we can still control the size of the figures this way
from pylab import rcParams
rcParams['figure.figsize'] = 3, 3
else:
# This block will be run if there is no ImportError
sns.set_style("ticks")
sns.set_context("paper",rc={"figure.figsize": (2, 2)})
# We create a "named tuple" data type to store the exprimental data in.
# In[3]:
from collections import namedtuple
ExperimentData = namedtuple('ExperimentData', ['T', 'cA_start', 'times', 'cA'])
def plot_experiment(e):
"""
Plots the experimental data provided in 'e'
which should be of the type ExperimentData.
"""
plt.plot(0, e.cA_start, 'ko')
plt.plot(e.times, e.cA,':o', label="T={:.0f}K".format(e.T))
plt.ylim(0,)
plt.ylabel('$C_A$ (mol/L)')
plt.xlabel('time (s)')
plt.legend()
# Now here are the data from your three experiments:
# In[4]:
from numpy import array
experiments = [ExperimentData(T=298.15,
cA_start=10.0,
times=array([ 10, 20, 30, 40, 50, 60, 70, 80, 90, 100]),
cA=array([ 8.649, 7.441, 7.141, 6.366, 6.215, 5.990, 5.852, 5.615, 5.481 , 5.644])),
ExperimentData(T=308.15,
cA_start=10.0,
times=array([ 10, 20, 30, 40, 50, 60, 70, 80, 90, 100]),
cA=array([ 7.230, 6.073, 5.452, 5.317, 5.121, 4.998, 4.951, 4.978, 5.015, 5.036])),
ExperimentData(T=323.15,
cA_start=10.0,
times=array([ 10, 20, 30, 40, 50, 60, 70, 80, 90, 100]),
cA=array([ 5.137, 4.568, 4.548, 4.461, 4.382, 4.525, 4.483, 4.565, 4.459, 4.635])),
]
for i,e in enumerate(experiments):
print("Experiment {} was at T={}K and ran for {} seconds".format(i, e.T, e.times[-1]))
plot_experiment(e)
# In[5]:
ParameterSet = namedtuple('ParameterSet', ['logA', 'Ea', 'dH', 'dS'])
# This is a sensible starting guess for your fitting
starting_guess = ParameterSet(
logA = 6. , # base-ten logarithm of A in s^-1
Ea = 45. , # Ea in kJ/mol
dH = -10. , # ∆H in kJ/mol
dS = -50. # ∆S in J/mol/K
)
# This should end up with your optimized parameters
optimized_parameters = ParameterSet(0,0,0,0)
# This should end up with your standard errors (one sigma)
# for the uncertainties in the fitted parameters.
# i.e. there should be a 68% chance the true value is
# at least this close to your optimized parameter.
standard_errors = ParameterSet(0,0,0,0)
# Ok, now insert some cells to determine the optimized_parameters and their standard_errors.
# In[6]:
# Finish your notebook with this cell
print(starting_guess)
print(optimized_parameters)
print(standard_errors)