Figure A.3:

Solution to first-order differential equation dc_A/dt=-kc_A, and sensitivities S_1=d c_A/d k and S_2=d c_A/d c_{A0}.

Code for Figure A.3

Text of the GNU GPL.

main.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
from paresto import Paresto
import numpy as np
import matplotlib.pyplot as plt
import os
plt.ion()

# Define parameters
ca0 = 2
k = 1
tfinal = 5
nts = 100
tout = np.linspace(0, tfinal, nts)

# Initialize the model
pe = Paresto()
pe.transcription = 'shooting'
pe.x = ['ca']
pe.p = ['k']
pe.d = ['m_ca']
pe.tout = tout

# Define the ODE function
pe.ode = lambda t, y, p: [-p.k * y.ca]

# Define the least squares function
pe.lsq = lambda t, y, p: [y.ca - y.m_ca]

# Initialize the Paresto object
pe.initialize()

# Initialize state and parameters
x0 = ca0
y_ac, z_ac = pe.simulate({'k': k, 'ca': x0})

# Initialize parameters and bounds
pe.ic.k = k
pe.ic.ca = ca0

pe.lb_ub()
# lb and ub are copies of ic (all equal), so parameters are fixed

# Optimize the parameters
est, y, p = pe.optimize(y_ac)

conf = pe.confidence(est, 0.95)

table = np.column_stack((tout, y['ca'].reshape(nts), est.dca_dk.reshape(nts), est.dca_dca0.reshape(nts)))

print('Estimated parameters')
print(f'  ca_final = {table[-1, 1]:.6g}')
print(f'  dca_dk_final = {table[-1, 2]:.6g}')
print(f'  dca_dca0_final = {table[-1, 3]:.6g}')

# Save the data to a .dat file
np.savetxt('Sfirstorder.dat', table, fmt='%.6e')

# Conditional plotting
if not os.getenv('OMIT_PLOTS') == 'true':
    plt.plot(tout, y['ca'].reshape(nts), label='y.ca')
    plt.plot(tout, est.dca_dk.reshape(nts), label='est.dca_dk')
    plt.plot(tout, est.dca_dca0.reshape(nts), label='est.dca_dca0')
    plt.show()