Figure 9.21:

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

Code for Figure 9.21

Text of the GNU GPL.

main.m


 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
%
% jbr,  4/8/18
%


model = struct;
model.x = {'ca'};
model.p = {'k', 'n'};
model.d = {'m_ca'};

model.ode = @(t, y, p) {-p.k * y.ca^p.n};

model.lsq = @(t, y, p) {y.ca - y.m_ca};

tfinal = 5;
nts    = 100;
tout = linspace(0,tfinal,nts);
model.tout = tout;

pe = paresto(model);

kac   = 0.5;
nac   = 2.5;
ca0ac = 2;
thetaac = [kac; nac; ca0ac];

x_ac = ca0ac;
p_ac = [kac; nac];

y_ac = pe.simulate(0, x_ac, p_ac);

% add measurement noise
measvar = 1e-2;
measstddev = sqrt(measvar);
rng(0);
noise = measstddev*randn(1,nts);

y_noisy = y_ac + noise;


% Initial guess, upper and lower bounds for the estimated parameters
theta0 = [0.5; 2.5; 2];
lbtheta = 0.5*[p_ac; x_ac(:)];
ubtheta = 1.5*[p_ac; x_ac(:)];

% Estimate parameters
[est, y, p] = pe.optimize(y_noisy, theta0, lbtheta, ubtheta);
est.d2f_dtheta2;

% Also calculate confidence intervals with 95 % confidence
theta_conf = pe.confidence(est, 1:numel(theta0), 0.95);

disp('Optimal parameters')
est.theta
disp('Confidence interval')
theta_conf

if (~ strcmp (getenv ('OMIT_PLOTS'), 'true')) % PLOTTING
% Plot simulated measurement values
plot(tout, y_noisy, 'ro', tout, y.ca)
% TITLE
end % PLOTTING

table  = [tout; y_noisy; y.ca; y_ac]';
save nthorder.dat table;