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
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
%
% jbr,  4/8/18
%


model = struct;
%model.nlp_solver_options.ipopt.linear_solver = 'ma27';
model.x = {'ca'};
model.p = {'k', 'n'};
model.d = {'m_ca'};
model.transcription = 'shooting';

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);
randn('seed',0);
noise = measstddev*randn(1,nts);

y_noisy = y_ac + noise;


% Initial guess, upper and lower bounds for the estimated parameters
small = 1e-3;
large  = 5;
np = 3;

theta0 = struct();
theta0.k = kac;
theta0.n = nac;
theta0.ca = ca0ac;

ubtheta = struct();
ubtheta.k = large;
ubtheta.n = large;
ubtheta.ca = large;

lbtheta = struct();
lbtheta.k = small;
lbtheta.n = small;
lbtheta.ca = small;

% % 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, 0.95);

conf = pe.confidence(est,  0.95);

disp('Estimated parameters')
disp(est.thetavec)
disp(est.theta)
disp('Bounding box intervals')
disp(conf.bbox)

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;

nthorder.dat

  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
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
# Created by Octave 6.2.0, Sat Oct 16 11:49:34 2021 PDT <jbraw@telemark>
# name: table
# type: matrix
# rows: 100
# columns: 4
 0 1.7846141099929809 1.89004579624681 2
 0.050505050505050504 1.7001177863100367 1.78125430277643 1.8687931136110623
 0.10101010101010101 1.7235388638246931 1.6868781030464555 1.7571908624638
 0.15151515151515152 1.7512545863511977 1.6040856320626484 1.6608821014287887
 0.20202020202020202 1.5821108909491211 1.5307577302173225 1.5767690768632798
 0.25252525252525254 1.5399928894273969 1.465274145448924 1.502552252673361
 0.30303030303030304 1.398246759497612 1.4063730475402922 1.4364909830922776
 0.35353535353535354 1.3946250408548502 1.3530558489059628 1.3772429972547677
 0.40404040404040403 1.3728188140907482 1.3045210663899967 1.3237497194328502
 0.45454545454545453 1.3141388589785745 1.2601173312970941 1.2751653576539208
 0.50505050505050508 1.1951418974549601 1.2193093421812653 1.2308067389876673
 0.55555555555555558 1.2188458056415998 1.1816527605392029 1.1901149095262968
 0.60606060606060608 1.0283218028231784 1.1467754081936539 1.1526277305766268
 0.65656565656565657 1.1422501860060568 1.1143629842116849 1.1179595617451543
 0.70707070707070707 1.1639587587025149 1.0841480752484765 1.0857856339123233
 0.75757575757575757 1.1232228475034778 1.0559016007853379 1.0558306413114611
 0.80808080808080807 0.97992433021289382 1.0294260824473422 1.0278600187085298
 0.85858585858585856 0.85081440321933277 1.004550296390287 1.0016715060902739
 0.90909090909090906 0.99630231476345044 0.98112498609487053 0.97708974695720652
 0.95959595959595956 0.83317345218130423 0.95901939654142665 0.95396250323721243
 1.0101010101010102 0.77521354624161432 0.93811845067607658 0.93215692707428643
 1.0606060606060606 0.78414055049612152 0.91832043256305573 0.91155654132558928
 1.1111111111111112 1.0347335696735354 0.89953507353269158 0.89205827718162256
 1.1616161616161615 0.88813728162187788 0.88168196131100285 0.87357072302240579
 1.2121212121212122 0.72395056716215422 0.86468920985860465 0.85601161948454196
 1.2626262626262625 0.85759768125160796 0.84849234106619664 0.83930840727432832
 1.3131313131313131 0.70030077494409204 0.83303333969365945 0.82339672840860012
 1.3636363636363635 0.91003458342495214 0.81825985081784836 0.80821852526607763
 1.4141414141414141 0.8274335198580034 0.80412449516228057 0.79372116555559713
 1.4646464646464645 0.67075863309388906 0.79058428245245937 0.77985692449098376
 1.5151515151515151 1.0085959777597946 0.77760010669217794 0.7665825710062546
 1.5656565656565655 0.81133892409616115 0.76513631022573891 0.75385894052796965
 1.6161616161616161 0.73172340036945716 0.75316030581614735 0.74165024103241339
 1.6666666666666667 0.80536054094074505 0.74164224786407407 0.72992376406429549
 1.7171717171717171 0.76810998685287835 0.73055474541877208 0.71864989168571836
 1.7676767676767677 0.69632857514452129 0.71987261086824095 0.70780178828071738
 1.8181818181818181 0.62165723830971353 0.70957263920207803 0.6973538914373838
 1.8686868686868687 0.75563359047516176 0.69963341356329389 0.68728348399742434
 1.9191919191919191 0.76841961758526089 0.69003513348142731 0.67756931202800985
 1.9696969696969697 0.73954010986395891 0.68075946273716326 0.66819150232382829
 2.0202020202020203 0.58437435880391631 0.67178939427094031 0.65913158908574609
 2.0707070707070705 0.60474180953807743 0.6631091299326346 0.65037291722602752
 2.1212121212121211 0.58330737236324726 0.65470397319061591 0.64189926269833031
 2.1717171717171717 0.6955960865764903 0.64656023318772771 0.63369638879588119
 2.2222222222222223 0.77653852484799646 0.63866513875822961 0.62575046799756306
 2.2727272727272725 0.72034966852638727 0.63100676121093258 0.61804891731712819
 2.3232323232323231 0.71553715973347032 0.62357394484563877 0.6105799963995775
 2.3737373737373737 0.3534458949800543 0.61635624430753089 0.60333269383869692
 2.4242424242424243 0.48792752268894685 0.60934386800133522 0.5962967062292911
 2.4747474747474749 0.75277109716605739 0.60252762688721329 0.58946235750388698
 2.5252525252525251 0.50217884473213892 0.59589888806614788 0.58282075814614043
 2.5757575757575757 0.58048439934090512 0.58944953263632371 0.57636325000957389
 2.6262626262626263 0.69482487849234409 0.58317191736551965 0.5700815956973917
 2.6767676767676769 0.38126331797289026 0.57705883977939487 0.56396833887743125
 2.7272727272727275 0.51402130674863367 0.57110350631305107 0.55801632952237634
 2.7777777777777777 0.69985346772595158 0.56529950321447076 0.5522187707655023
 2.8282828282828283 0.63360449514769646 0.55964076992429801 0.5465694232978735
 2.8787878787878789 0.35703588486526816 0.55412157468767609 0.54106235266540859
 2.9292929292929295 0.47778794080332632 0.54873649218116971 0.53569188624934072
 2.9797979797979797 0.52456338970706595 0.54348038296170764 0.53045275627420085
 3.0303030303030303 0.56301017813248255 0.53834837456544837 0.52533991448444939
 3.0808080808080809 0.45128212527948314 0.53333584410290702 0.52034854845720224
 3.1313131313131315 0.62333572043856567 0.52843840221291238 0.51547409429034174
 3.1818181818181817 0.34921027419274453 0.52365187825228077 0.51071225402062537
 3.2323232323232323 0.52944170476012642 0.51897230661074856 0.50605885178857257
 3.2828282828282829 0.58052590407146243 0.51439591405191409 0.50151003516925596
 3.333333333333333 0.43725855081896103 0.50991910799087214 0.49706205457071578
 3.3838383838383841 0.58348361141191507 0.50553846562804428 0.49271138197885533
 3.4343434343434343 0.50073293081884174 0.50125072386658065 0.48845466754798683
 3.4848484848484849 0.4875771034225006 0.49705276994767933 0.48428871219001007
 3.5353535353535355 0.69132565223606857 0.49294163274444114 0.48021045886906416
 3.5858585858585856 0.46250258043068349 0.48891447466043753 0.47621699973600806
 3.6363636363636367 0.42908497573971877 0.48496858408418086 0.47230555238366256
 3.6868686868686869 0.5858758354128718 0.48110136835515283 0.46847344874752755
 3.7373737373737375 0.37507004387966303 0.47731034720107712 0.46471818335643916
 3.7878787878787881 0.47560289624626634 0.47359314660971913 0.46103735688621994
 3.8383838383838382 0.61995154109007844 0.46994749310176481 0.45742861714369776
 3.8888888888888888 0.61379734581378298 0.46637120837424018 0.45388969486621217
 3.9393939393939394 0.59044475257347939 0.46286220428659203 0.45041837394188755
 3.9898989898989896 0.41789253231626311 0.45941847816392467 0.4470126974258975
 4.0404040404040407 0.21445804132038609 0.45603810839405473 0.44367070687824745
 4.0909090909090908 0.55829741986490555 0.4527192502969925 0.44039043219782181
 4.1414141414141419 0.53838327382662299 0.44946013224723486 0.43717005704500683
 4.191919191919192 0.43235632663326873 0.44625905203085897 0.43400791154908552
 4.2424242424242422 0.49444325496784625 0.44311437342087112 0.4309023015700954
 4.2929292929292933 0.52068752654382233 0.44002452295558026 0.42785164125749114
 4.3434343434343434 0.42723584525465047 0.43698798690598972 0.424854376131585
 4.3939393939393936 0.44594329541630079 0.43400330841928958 0.421909051594203
 4.4444444444444446 0.47158187423125575 0.43106908482654333 0.41901426110640833
 4.4949494949494948 0.4178217854695972 0.4281839651035752 0.41616862829500995
 4.5454545454545459 0.51289529078425256 0.42534664747490508 0.41337082856119955
 4.595959595959596 0.36476659723131111 0.42255587715134246 0.41061959215013438
 4.6464646464646462 0.34152588819065544 0.41981044419255287 0.40791369055309745
 4.6969696969696972 0.35256325914776776 0.41710918148655574 0.40525195314801188
 4.7474747474747474 0.4139349913547416 0.41445096283869687 0.40263321831325488
 4.7979797979797976 0.32805982035831593 0.41183470116318671 0.40005636614994194
 4.8484848484848486 0.48453268532476579 0.40925934677078374 0.39752037052831801
 4.8989898989898988 0.5473400362836881 0.40672388574667179 0.39502418125744299
 4.9494949494949498 0.3319781897017397 0.40422733841298952 0.39256679987401144
 5 0.30154472321212739 0.4017687578708608 0.39014726608932465