Figure 9.21

Figure 9.21:

Experimental measurement and best parameter fit for nth-order kinetic model, r=k c_A^n.

Code for Figure 9.21

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
63
64
65
66
67 # [makes] nthorder.dat
# nth-order kinetics parameter estimation; mirrors nthorder.m structure.
from paresto import Paresto
from misc import save_ascii
import numpy as np

# Model
pe = Paresto()
pe.print_level = 0
pe.transcription = 'shooting'
pe.x = ['ca']
pe.p = ['k', 'n']
pe.d = ['m_ca']
pe.ode = lambda t, y, p: [-p.k * y.ca**p.n]
pe.lsq = lambda t, y, p: [y.ca - y.m_ca]

tfinal = 5.0
nts    = 100
tout   = np.linspace(0, tfinal, nts)
pe.tout = tout

pe.initialize()

# True parameters used to generate the data
kac, nac, ca0ac = 0.5, 2.5, 2.0

# Noise-free trajectory
pin_ac = {'k': kac, 'n': nac, 'ca': ca0ac}
y_ac, _ = pe.simulate(pin_ac)
y_ac = np.array(y_ac).flatten()

# Add measurement noise
measvar    = 1e-2
measstddev = np.sqrt(measvar)
np.random.seed(0)
noise   = measstddev * np.random.randn(nts)
y_noisy = y_ac + noise

# Initial guess (true values, like nthorder.m) and bounds
small = 1e-3
large = 5.0
pe.ic.k  = kac
pe.ic.n  = nac
pe.ic.ca = ca0ac

pe.lb_ub()
pe.lb.k = small;  pe.ub.k = large
pe.lb.n = small;  pe.ub.n = large
pe.lb.ca = small; pe.ub.ca = large

# Parameter estimation
est, yy, pp = pe.optimize(y_noisy.reshape(1, nts, 1))

# 95% confidence intervals
conf = pe.confidence(est, 0.95)

print('Estimated parameters')
for k, v in est.theta.items():
    print(f'  {k}: {float(v):.4f}')
print('Bounding box intervals (95%)')
for k, v in conf.bbox.items():
    print(f'  {k}: {float(v):.4f}')

# Save table for gnuplot: time, noisy data, fit, true
y_fit = np.array(yy['ca']).flatten()
table = np.column_stack([tout, y_noisy, y_fit, y_ac])
save_ascii('nthorder.dat', table)