<jupyter_start><jupyter_text>Sensitivity of Concentration to Rate Constant and Initial Conditionjbr book reference: [Figure A.3](https://sites.engineering.ucsb.edu/~jbraw/chemreacfun/fig-html/appendix/fig-A-3.html)Reaction is $A \xrightarrow[]{k_1} B$Analytical solution:$$c_A(t) = c_{A0}e^{-k_1t}$$Sensitivity:$$\begin{eqnarray}\frac{dc_A}{dk_1} &=& -tc_{A0}e^{-k_1t} \\\frac{dc_A}{dc_{A0}} &=& e^{-k_1t}\end{eqnarray}$$In this example, $c_{A0} = 2\;mol/L$, $k_1 = 1\;min^{-1}$<jupyter_code>import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import pyomo.environ as pyo
from kipet import KipetModel
from kipet.library.core_methods.EstimabilityAnalyzer import *
# Though the sensitivity calculation doesn't need data, in KIPET sensitivity is part of Estimability routine 
# which needs concentration data. So the concentration data is simulated first prior to getting sensitivity
kipet_model = KipetModel()

r1 = kipet_model.new_reaction('rxn1')
    
r1.add_component('A', state = 'concentration', init = 2.0)
r1.add_parameter('k1', init = 1.0)
    
def rule_odes(m,t):
    exprs = dict()
    exprs['A'] = -m.P['k1']*m.Z[t,'A']
    return exprs

r1.add_equations(rule_odes)
r1.set_times(0.0, 10.0)
r1.settings.collocation.nfe = 60
r1.settings.collocation.ncp = 1
r1.settings.collocation.scheme = 'LAGRANGE-RADAU'
r1.settings.solver.linear_solver = 'ma27'
r1.simulate()
# data
C_data = r1.results.Z
C_data
fig, ax = plt.subplots()
ax.scatter(C_data.index, C_data['A'])
# build the model again for getting sensitivity
kipet_model = KipetModel()

r1 = kipet_model.new_reaction('rxn1')
    
r1.add_component('A', state = 'concentration', init = 2.0)

r1.add_parameter('k1', init = 1.0, bounds = (0.1, 10)) # rate constant
r1.add_parameter('ca0', init = 2.0, bounds = (0.0, 5.0)) # initial condition listed as parameters since we want sensitivity w.r.t ca0

def rule_odes(m,t):
    exprs = dict()
    exprs['A'] = -m.P['k1']*m.Z[t,'A']
    return exprs

r1.add_equations(rule_odes)
r1.set_times(0.0, 10.0)

r1.add_dataset('C_data', category = 'concentration', data = C_data)
r1.create_pyomo_model()
# specify values of parameters at which sensitivity is desired
r1.model.P['k1'].value = 1.0
r1.model.P['ca0'].value = 2.0

# we need an explicit constraint for conc(A) at time 0 = ca0
r1.model.ca0_cons = pyo.Constraint(expr = r1.model.Z[0, 'A'] == r1.model.P['ca0'])
# since we impose a constraint above for conc(A) at time 0, we need to relax the automatic initial condition 
# constraint on conc(A) at time 0 = 2 from the model
r1.model.init_conditions_c['A'].deactivate()
r1.model.init_conditions_c.pprint()
e_analyzer = EstimabilityAnalyzer(r1.model)
e_analyzer.apply_discretization('dae.collocation', nfe = 60, ncp = 2, scheme = 'LAGRANGE-RADAU')
sigmas = {'A': 1} # we can just keep this as one since we want unscaled sensitivity
param_uncertainties = {'k1':1,'ca0':1} # we can just keep this as one since we want unscaled sensitivity
dsdp, idx_to_param = e_analyzer.get_sensitivities_for_params(tee=True, sigmasq=sigmas)
# get the sensitivities and compare with analytical sensitivities
tlist = np.array(e_analyzer.model.allmeas_times)
sens = dsdp.copy()
dca_dk1 = sens[:, 0]
dca_dca0 = sens[:, 1]
# analytical sensitivity
dca_dk1_act = -tlist * 2 * np.exp(-tlist) # sensitivity dca/dk1 = -t*ca0*exp(-k1*t), ca0 = 2, k1 = 1
dca_dca0_act = np.exp(-tlist) # sensitivity dca/dca0 = exp(-k1*t), k1 = 1
fig, ax = plt.subplots()
ax.plot(tlist, dca_dk1, label = 'kipet')
ax.plot(tlist, dca_dk1_act,'--r', label = 'analytical')
ax.set_xlabel('time (min)')
ax.set_ylabel('dca/dk1')
ax.legend()
fig, ax = plt.subplots()
ax.plot(tlist, dca_dca0, label = 'kipet')
ax.plot(tlist, dca_dca0_act,'--r', label = 'analytical')
ax.set_xlabel('time (min)')
ax.set_ylabel('dca/dca0')
ax.legend()<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.1534640522875817, 'fraction_numerical': 0.03372549019607843, 'mean_word_length': 5.131410256410256, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 5, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '130346', 'license': 'MIT', 'n_tokens_mistral': 1579, 'n_tokens_neox': 1499, 'n_words': 349, 'path': 'kipet_paresto/jbr_Sfirstorder.ipynb', 'repo_name': 'notesofdabbler/learn_kipet', 'size': '3825'}