Figure A.10:

Fit of model to measurements using estimated parameters.

Code for Figure A.10

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
66
67
68
69
70
71
72
%
% Estimate parameters for the A->B->C model using paresto.m
% jbr,  4/2018
%

%Model
% This m-file loads data file 'ABC_data.dat'.
model = struct;
model.nlp_solver_options.ipopt.linear_solver = 'ma27';
model.x = {'ca', 'cb', 'cc'};
model.p = {'k1', 'k2'};
model.d = {'m_ca', 'm_cb', 'm_cc'};


model.ode = @massbal_ode;

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

tfinal = 6;
nplot = 75;
tplot = linspace(0, tfinal, nplot)';

k1 = 2; k2 = 1;
p_ac = [k1; k2];

ca0 = 1; cb0 = 0; cc0 = 0;
x_ac = [ca0; cb0; cc0];

thetaic   = [0.5; 3; x_ac];
thetalb   = [1e-4; 1e-4; 0.9*x_ac];
thetaub   = [10; 10; 1.1*x_ac];
est_ind = 1:2;

% load measurements from a file
tabledat = load ('ABC_data.dat');
tmeas = tabledat(:,1);
ymeas = tabledat(:,2:4);

[tout,~,ic] = unique([tmeas; tplot]);
% index of times at which measurement is made
meas_ind = ic(1:numel(tmeas));
model.tout = tout;
% interploate measurement onto new grid
y_noisy = interp1(tmeas, ymeas, tout, 'previous');
y_noisy(isnan(y_noisy)) = 0.;
model.lsq_ind = meas_ind'; % only use index of measurement times in objective function

% Create a paresto instance
pe = paresto(model);

% estimate the parameters
est = pe.optimize(y_noisy', thetaic, thetalb, thetaub);

% Also calculate confidence intervals with 95 % confidence
theta_conf = pe.confidence(est, est_ind, 0.95);

disp('Estimated Parameters and Bounding Box')
[est.theta(est_ind)  theta_conf]

%plot the model fit to the noisy measurements
figure(1);
plot(model.tout, est.x, tmeas, ymeas', 'o');

% if 2-d, plot confidence ellipse
% if (length(objective.estflag) == 2)
%   [xe, ye] = ellipse(estimates.Pinv, estimates.b, 100, estimates.parest);
%   figure(2);
%   plot(xe, ye)
% end

tableest = [model.tout, est.x'];
save ABC.dat tableest tabledat

ABC_data.dat


 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
   0.000   0.957  -0.031  -0.015
   0.263   0.557   0.330   0.044
   0.526   0.342   0.512   0.156
   0.789   0.224   0.499   0.310
   1.053   0.123   0.428   0.454
   1.316   0.079   0.396   0.556
   1.579   0.035   0.303   0.651
   1.842   0.029   0.287   0.658
   2.105   0.025   0.221   0.750
   2.368   0.017   0.148   0.854
   2.632  -0.002   0.182   0.845
   2.895   0.009   0.116   0.893
   3.158  -0.023   0.079   0.942
   3.421   0.006   0.078   0.899
   3.684   0.016   0.059   0.942
   3.947   0.014   0.036   0.991
   4.211  -0.009   0.014   0.988
   4.474  -0.030   0.036   0.941
   4.737   0.004   0.036   0.971
   5.000  -0.024   0.028   0.985

massbal_ode.m


1
2
3
4
function xdot = massbal_ode(t, x, p)
  r1 = p.k1*x.ca;
  r2 = p.k2*x.cb;
  xdot = {-r1, r1-r2, r2};