-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathPETpc_kinetic_model.m
More file actions
497 lines (420 loc) · 41.5 KB
/
PETpc_kinetic_model.m
File metadata and controls
497 lines (420 loc) · 41.5 KB
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
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
% Applied physical chemistry project
% Lorenzo Paggetta
% A.Y. 2022/23
clear all
close all
clc
options_ga= optimoptions('ga','ConstraintTolerance',1e-6,'FunctionTolerance', 1e-8,...
'MaxGeneration',100,'UseParallel',true,'PopulationSize',400,...
'PlotFcn', {@gaplotbestf,@gaplotstopping,@gaplotrange,@gaplotselection,@gaplotscorediversity});
options_ga2 = optimoptions('ga','HybridFcn',{@fmincon},'ConstraintTolerance',1e-6,'FunctionTolerance', 1e-8,...
'MaxGeneration',100,'UseParallel',true,'PopulationSize',400,...
'PlotFcn', {@gaplotbestf,@gaplotstopping,@gaplotrange,@gaplotselection,@gaplotgenealogy,@gaplotscorediversity});
options_lsq = optimoptions('lsqnonlin','MaxIter',500,'UseParallel',false,'FunctionTolerance',1e-9,...
'OptimalityTolerance',1e-9,'StepTolerance',1e-9);
%% Data
% figure 5
A0=[472.1519 500 593.6709 726.5823 807.5949 965.8228];
conv5=[20.9032258064516 32.9032258064516 55.9354838709677 73.1612903225806 79.9354838709677];
% figure 6
conv6_m_prod=[0.602409638554217 7.53012048192771 22.8915662650602 28.9156626506024 43.9759036144578 55.7228915662651 73.7951807228916 86.7469879518072 91.5662650602410];
mass_prod=[0 1.59774436090226 4.60526315789474 5.63909774436090 8.36466165413534 10.9022556390977 12.7819548872180 16.3533834586466 17.1052631578947];
conv6_m_theo=[0.903614457831325 7.53012048192771 22.8915662650602 28.6144578313253 44.2771084337349 55.7228915662651 74.0963855421687 86.4457831325301 91.2650602409639];
mass_theo=[0 1.78571428571429 4.88721804511278 6.10902255639098 9.21052631578947 11.7481203007519 15.4135338345865 18.0451127819549 18.9849624060150];
conv6_yield=[8.73493975903615 25 31.0240963855422 46.9879518072289 59.0361445783133 78.0120481927711 90.9638554216867 96.0843373493976];
yield=[93.8223938223938 95.7528957528958 91.8918918918919 89.9613899613900 93.0501930501931 88.0308880308880 89.9613899613900 89.9613899613900];
% figure 7
% model
time_f7_170C=[2.55319148936170 5.10638297872340 8.93617021276596 15.3191489361702 20.4255319148936 23.6170212765957 26.8085106382979 31.2765957446809 35.1063829787234 37.6595744680851 40.2127659574468 42.1276595744681 45.3191489361702 47.2340425531915 51.0638297872340 54.2553191489362 56.8085106382979 61.2765957446809 66.3829787234043 70.8510638297872 75.4609929078014 80.6382978723404 85.8156028368794 90.7801418439716 95.8156028368794 102.624113475177 107.163120567376 111.347517730496 115.106382978723 118.510638297872 121.560283687943 124.468085106383 127.092198581560 129.361702127660 131.347517730496 133.262411347518 135.319148936170 137.375886524823 139.219858156028 140.992907801418 142.695035460993 144.255319148936 145.886524822695 147.517730496454 149.148936170213 150.851063829787 152.482269503546 154.113475177305 155.815602836879 157.588652482270 159.503546099291 161.489361702128 163.333333333333 164.822695035461 166.524822695035 169.148936170213 172.198581560284 175.106382978723 178.226950354610 181.702127659574 185.602836879433 189.858156028369 194.397163120567 198.936170212766 203.617021276596 208.510638297872 213.617021276596 218.723404255319 223.829787234043 236.453900709220 241.276595744681 246.382978723404 251.560283687943 256.666666666667 263.758865248227 268.723404255319 273.900709219858 279.007092198582 284.113475177305 289.219858156028 ...
294.326241134752];
conv7_170C=[0 0.360537528679122 1.40937397574566 2.35988200589971 4.42477876106195 5.70304818092429 7.66961651917404 9.43952802359882 12.0943952802360 14.8803670927565 16.8797115699771 18.5512946574893 23.4021632251721 25.4998361193051 27.5975090134382 29.6951819075713 31.8256309406752 33.9233038348083 36.0537528679122 38.1842019010161 40.3802032120616 42.5434283841364 43.9528023598820 46.3126843657817 48.1153720091773 50.0491642084562 52.1468371025893 54.3100622746640 56.4405113077680 60.6686332350049 62.7663061291380 64.8639790232711 66.8960996394625 68.9282202556539 69.6165191740413 70.7964601769912 72.2713864306785 75.3851196329072 77.3189118321862 79.2199278924943 81.0881678138315 82.8908554572272 84.6607669616519 86.2667977712225 87.9056047197640 89.0855457227139 90.2654867256637 90.9865617830220 92.2648312028843 93.3792199278925 93.3464437889217 94.3952802359882 95.3130121271714 95.1491314323173 95.7718780727630 96.0668633235005 97.1812520485087 97.2140281874795 97.6728941330711 97.9023271058669 97.5090134382170 98.0989839396919 98.4595214683710 98.8856112749918 99.2133726647001 99.2461488036709 99.1805965257293 99.2461488036709 99.4755817764667 99.6066863323501 99.2789249426418 99.8688954441167 99.8361193051459 99.8361193051459 99.7377908882334 99.8033431661750 99.7377908882334 99.8361193051459 99.7377908882334 99.8688954441167 99.8688954441167 99.7050147492625 99.7050147492625 99.7050147492625 99.7050147492625 99.7050147492625 99.7050147492625 99.8688954441167 99.7377908882334 99.8033431661750 99.8688954441167 99.8688954441167 ...
99.8361193051459];
time_f7_175C=[1.27659574468085 3.19148936170213 5.74468085106383 8.29787234042553 12.1276595744681 14.6808510638298 17.2340425531915 19.4326241134752 23.7588652482270 28.7234042553192 33.9007092198582 38.7943262411348 43.6170212765957 48.4397163120567 53.0496453900709 57.5886524822695 61.9148936170213 66.0992907801418 70.0709219858156 73.7588652482270 77.2340425531915 80.5673758865248 83.9007092198582 87.0212765957447 89.9290780141844 92.6241134751773 95.3191489361702 98.0851063829787 100.638297872340 103.049645390071 105.460992907801 107.872340425532 110.141843971631 112.340425531915 114.609929078014 116.950354609929 119.219858156028 121.631205673759 124.042553191489 126.382978723404 128.723404255319 131.134751773050 133.617021276596 136.312056737589 139.078014184397 141.773049645390 144.609929078014 147.588652482270 150.780141843972 154.184397163121 157.659574468085 161.276595744681 165.248226950355 169.432624113475 173.687943262411 178.226950354610 182.907801418440 187.730496453901 192.553191489362 197.517730496454 202.482269503546 207.517730496454 212.553191489362 217.588652482270 222.624113475177 227.446808510638 232.553191489362 237.659574468085 242.765957446809 247.872340425532 252.978723404255 258.085106382979 263.191489361702 268.297872340426 273.404255319149 278.510638297872 283.617021276596 288.723404255319 ...
293.829787234043];
conv7_175C=[1.47492625368732 1.17994100294985 1.47492625368732 1.47492625368732 1.47492625368732 1.47492625368732 1.76991150442478 1.47492625368732 1.76991150442478 1.47492625368732 1.47492625368732 1.76991150442478 1.47492625368732 1.47492625368732 1.47492625368732 1.47492625368732 1.47492625368732 1.76991150442478 1.47492625368732 1.63880694854146 2.16322517207473 2.29432972795805 2.49098656178302 3.04818092428712 3.47427073090790 4.85086856768273 5.67027204195346 7.21075057358243 8.91510980006555 10.6850213044903 12.5532612258276 14.4870534251065 16.4863979023271 18.5185185185185 20.6161914126516 22.7794165847263 24.9098656178302 27.0730907899049 29.3018682399213 31.5306456899377 33.7594231399541 35.9882005899705 38.2497541789577 40.4785316289741 42.7073090789905 44.9033103900361 47.1320878400524 49.3608652900688 51.5568666011144 53.7856440511308 56.0144215011472 58.1448705342511 60.1442150114716 62.1107833497214 64.0773516879712 66.0439200262209 68.0104883644707 69.9442805637496 71.7797443461160 73.5168797115700 74.9918059652573 76.1717469682072 77.1222549983612 78.0399868895444 78.8593903638151 79.4165847263192 79.7443461160275 79.9410029498525 80.0721075057358 80.1048836447067 80.4654211733858 80.4981973123566 80.4654211733858 80.3998688954441 80.6293018682399 80.2687643395608 80.3670927564733 80.4654211733858 80.4981973123566 80.5309734513274 ...
80.5309734513274];
time_f7_180C=[-2.41134751773050 6.95035460992908 12.0567375886525 17.2340425531915 21.0638297872340 24.8936170212766 29.3617021276596 31.9148936170213 36.3829787234043 39.3617021276596 43.1914893617021 46.5248226950355 49.9290780141844 53.1205673758865 56.1702127659574 59.3617021276596 62.4113475177305 65.1063829787234 67.5177304964539 69.9290780141844 72.4113475177305 74.8936170212766 77.3049645390071 79.7163120567376 82.1985815602837 84.6099290780142 86.9503546099291 89.2198581560284 91.4184397163121 93.9007092198582 96.6666666666667 97.0212765957447 98.5815602836879 99.5744680851064 101.773049645390 103.475177304965 105.886524822695 108.510638297872 111.063829787234 113.758865248227 116.453900709220 119.007092198582 121.773049645390 124.539007092199 127.446808510638 130.709219858156 134.184397163121 137.730496453901 141.276595744681 145.035460992908 149.007092198582 153.191489361702 157.801418439716 162.056737588652 164.680851063830 170.709219858156 175.248226950355 180.070921985816 185.035460992908 190 195.319148936170 197.304964539007 204.255319148936 208.723404255319 217.021276595745 222.765957446809 226.595744680851 230.354609929078 234.113475177305 235.319148936170 239.219858156028 244.326241134752 249.432624113475 254.468085106383 262.978723404255 269.361702127660 274.468085106383 279.574468085106 284.680851063830 289.290780141844 293.617021276596...
298.581560283688 ];
conv7_180C=[2.65486725663717 2.65486725663717 2.65486725663717 2.35988200589971 2.65486725663717 2.65486725663717 2.65486725663717 2.81874795149131 3.63815142576205 3.96591281547034 4.29367420517863 4.88364470665356 5.57194362504097 6.35857096034087 7.30907899049492 8.45624385447394 9.70173713536545 11.0127826941986 12.4877089478859 14.1265158964274 15.8636512618814 17.6335627663061 19.4362504097017 21.3044903310390 23.2710586692888 25.2376270075385 27.2697476237299 29.2690921009505 31.3012127171419 33.3988856112750 35.4637823664372 37.5286791215995 39.6591281547034 41.7895771878073 43.8872500819404 45.9849229760734 48.0825958702065 50.1802687643396 52.2779416584726 54.3756145526057 56.4732874467388 58.5709603408719 60.6358570960341 62.6352015732547 64.6345460504753 66.6338905276958 68.6004588659456 70.5014749262537 72.3369387086201 74.1068502130449 75.8112094395280 77.4500163880695 78.9577187807276 80.3343166175025 81.6453621763356 82.7925270403147 83.7758112094395 84.5952146837103 85.2835136020977 85.8734841035726 86.3979023271059 86.8239921337267 87.1517535234349 87.4139626352016 87.7417240249099 88.3644706653556 88.4627990822681 88.3316945263848 88.2333661094723 88.2661422484431 88.3316945263848 88.3644706653556 88.3972468043264 88.4627990822681 88.5283513602098 88.5283513602098 88.5283513602098 88.5283513602098 ...
88.6266797771223];
time_f7_185C=[1.27659574468085 3.82978723404255 7.02127659574468 10.8510638297872 12.7659574468085 15.9574468085106 17.8723404255319 23.6170212765957 27.4468085106383 31.9148936170213 33.8297872340426 35.7446808510638 38.2978723404255 39.5744680851064 40.8510638297872 42.7659574468085 45.9574468085106 45.9574468085106 49.1489361702128 52.3404255319149 53.6170212765958 54.6099290780142 57.6595744680851 59.9290780141844 62.0567375886525 64.3262411347518 66.6666666666667 68.9361702127660 71.1347517730496 73.3333333333333 75.5319148936170 77.7304964539007 79.8581560283688 82.1276595744681 84.4680851063830 86.8085106382979 89.1489361702128 91.4893617021277 93.9007092198582 96.3829787234043 98.3687943262411 98.9361702127659 103.475177304965 103.900709219858 107.234042553192 110.070921985816 112.978723404255 116.099290780142 119.361702127660 122.836879432624 126.524822695035 130.709219858156 137.375886524823 141.914893617021 146.382978723404 150.354609929078 154.468085106383 158.936170212766 164.042553191489 167.872340425532 170.851063829787 175.886524822695 180.638297872340 192.624113475177 197.517730496454 202.624113475177 207.730496453901 212.765957446809 217.730496453901 222.695035460993 227.588652482270 232.411347517731 237.163120567376 242.127659574468 247.234042553191 252.340425531915 257.375886524823 262.340425531915 267.304964539007 272.411347517731 277.517730496454 282.624113475177 287.730496453901 292.836879432624 ...
297.943262411348];
conv7_185C=[0 1.60603080957063 1.96656833824975 2.94985250737463 4.12979351032448 5.60471976401180 6.48967551622419 7.96460176991150 9.14454277286136 10.4883644706654 12.2254998361193 14.1592920353982 15.8964274008522 17.7318911832186 19.6329072435267 21.5011471648640 23.3693870862012 25.3031792854802 27.3680760406424 29.4657489347755 31.5306456899377 33.5955424451000 35.6932153392330 37.7908882333661 39.8885611274992 41.9534578826614 44.0511307767945 46.1488036709276 48.2792527040315 50.4097017371354 52.4745984922976 53.2284496886267 54.5722713864307 55.7522123893805 56.8010488364471 58.9970501474926 61.0619469026549 63.0940675188463 65.1589642740085 67.1910848901999 69.1576532284497 71.1897738446411 73.1563421828909 75.1556866601114 77.0567027204195 78.7938380858735 80.4654211733858 82.1697803998689 83.8413634873812 85.4146181579810 86.8567682726975 88.2333661094723 89.4460832513930 90.6260242543428 91.4454277286136 92.1992789249427 93.4447722058341 94.2641756801049 94.8541461815798 95.4441166830547 95.7063257948214 96.1979678793838 96.4601769911504 96.7551622418879 97.0501474926254 97.6401179941003 97.6401179941003 97.8695509668961 97.7712225499836 98.0006555227794 97.9678793838086 97.7056702720420 97.7056702720420 98.0989839396919 98.2300884955752 98.2300884955752 98.2300884955752 98.1973123566044 98.1645362176336 98.6561783021960 98.4922976073419 ...
98.0662078007211];
time_f7_190C=[0.638297872340426 2.55319148936170 7.09219858156028 10.2127659574468 14.6808510638298 18.0141843971631 22.9787234042553 25.5319148936170 29.3617021276596 33.1205673758865 35.8865248226950 39.2198581560284 44.9645390070922 46.5248226950355 49.1489361702128 51.7730496453901 54.0425531914894 56.2411347517730 58.5106382978723 60.7801418439716 62.7659574468085 64.1843971631206 67.0212765957447 68.9361702127660 69.9290780141844 72.7659574468085 75.0354609929078 77.1631205673759 79.3617021276596 83.9007092198582 86.1702127659575 88.4397163120567 90.9929078014184 93.5460992907802 94.4680851063830 95.7446808510638 97.6595744680851 102.269503546099 105.248226950355 108.226950354610 111.347517730496 114.609929078014 118.085106382979 121.843971631206 125.460992907801 129.858156028369 132.127659574468 135.177304964539 139.432624113475 143.971631205674 144.751773049645 149.361702127660 154.113475177305 155.886524822695 160 160.070921985816 169.148936170213 173.617021276596 174.113475177305 176.382978723404 177.659574468085 181.489361702128 186.595744680851 191.489361702128 196.241134751773 200.992907801418 206.099290780142 211.205673758865 216.312056737589 221.418439716312 224.397163120567 228.936170212766 233.758865248227 236.312056737589 238.865248226950 241.418439716312 243.971631205674 246.524822695035 249.078014184397 251.631205673759 254.184397163121 254.680851063830 255.957446808511 257.872340425532 260.425531914894 262.978723404255 265.531914893617 267.943262411348 273.049645390071 278.156028368794 283.262411347518 288.368794326241 ...
293.475177304965];
conv7_190C=[0 0.294985250737463 1.17994100294985 2.06489675516224 3.24483775811209 4.12979351032448 5.01474926253687 7.37463126843658 9.43952802359882 11.5044247787611 12.6843657817109 14.1592920353982 15.6342182890855 16.8141592920354 17.6991150442478 18.8790560471976 20.3539823008850 21.8289085545723 23.8938053097345 25.0737463126844 27.1386430678466 28.3841363487381 30.0557194362504 32.1861684693543 34.4149459193707 36.5453949524746 38.6430678466077 40.7735168797116 42.9039659128155 45.0344149459194 47.1648639790233 49.2953130121272 51.4585381842019 53.5562110783350 55.6538839724680 57.7515568666011 59.8492297607342 61.9469026548673 64.0445755490003 66.1094723041626 67.7810553916749 68.1415929203540 71.8780727630285 72.4024909865618 74.1068502130449 76.0734185512947 78.0399868895444 79.9082268108817 81.7109144542773 83.4480498197312 85.0540806293019 86.4306784660767 88.8233366109472 90.1343821697804 91.4126515896427 92.7564732874467 93.5103244837758 94.1002949852507 94.9852507374631 95.2802359882006 95.7063257948214 95.7718780727630 96.5912815470338 97.1812520485087 97.0173713536545 97.3451327433628 97.6401179941003 97.6073418551295 97.4434611602753 97.4762372992462 97.8695509668961 98.2628646345461 98.2300884955752 97.9023271058669 97.7056702720420 97.8039986889544 98.0989839396919 98.3611930514585 98.3611930514585 98.2300884955752 98.1645362176336 98.1973123566044 98.2300884955752 98.3611930514585 ...
98.3939691904294];
% exp data
time_f7_170C_exp=[33.2 65.7 98.3 131 163 196 229 262 294];
conv7_170C_exp=[0.885 3.2400 5.0100 19.50 61.4 76.4 79.1 81.1 82.3];
time_f7_175C_exp=[32.6 98.3 163 229 294 ];
conv7_175C_exp=[7.08 28.9 79.1 87 89.1];
time_f7_180C_exp=[32.6 65.7 98.3 131 163 197 229 261 295];
conv7_180C_exp=[9.14 24.8 56 78.2 91.2 95.6 97.1 97.9 99.1];
time_f7_185C_exp=[33.2 98.3 163 229 294];
conv7_185C_exp=[11.8 68.1 94.1 97.9 99.1];
time_f7_190C_exp=[65.7 98.3 131 163 196 229 262 294];
conv7_190C_exp=[46 71.4 87.9 98.8 99.1 99.4 99.4 99.7];
% fig 8
oneoverT_exp=[0.00215977443609022 0.00218334586466165 0.00220748120300751 0.00223206766917293 ...
0.00225733082706766];
lnK_exp=[-4.48248587570621 -4.58418079096045 -4.72118644067796 -5.06581920903954 ...
-5.31723163841807];
T_span_high=448:1:473;
x_span_high=T_span_high.^(-1);
y_model_high=@(x_high) -5013.2*x_high+6.346;
T_span_low=431:1:456;
x_span_low=T_span_low.^(-1);
y_model_low=@(x_low) -11977*x_low+21.69;
% fig 9a
% model
time_f9a_180C=[0 0 3.14095013741657 7.45975657636435 10.0117785630153 15.6065959952886 20.2198665096192 24.0478994895956 28.2685512367491 30.9187279151943 32.6855123674912 35.0412249705536 39.0655673341186 42.5009815469179 45.9363957597173 49.4699646643110 53.1016882606989 56.6352571652925 60.0706713780919 63.4079308990970 66.6470357283078 69.8861405575187 73.2234000785238 76.5606595995289 79.7016097369454 81.4683941892423 83.0388692579505 85.0019630938359 88.4373773066353 92.0691009030232 95.8971338829996 99.7251668629761 103.553199842953 107.675696898312 112.288967412642 117.196702002356 122.300745975658 127.502944640754 132.901452689439 138.496270121712 144.189242245779 150.176678445230 156.360424028269 162.544169611307 169.022379269729 175.696898311739 182.469572045544 189.045936395760 193.462897526502 195.916764821358 200.431880643895 205.241460541814 212.014134275618 218.590498625834 225.265017667845 229.681978798587 232.332155477032 234.098939929329 235.865724381625 241.166077738516 246.466431095406 250 254.416961130742 260.600706713781 263.250883392226 266.784452296820 269.434628975265 ...
271.201413427562];
conv9a_180C=[0 0 0 0 0 0.728597449908925 3.27868852459016 6.01092896174863 8.65209471766849 10.2459016393443 11.4754098360656 13.2969034608379 15.8014571948998 18.7158469945355 21.6302367941712 24.4535519125683 27.2768670309654 30.1001821493625 32.9690346083789 35.8378870673953 38.7522768670310 41.6666666666667 44.5810564663024 47.4954462659381 50.4553734061931 51.9581056466302 53.3242258652095 54.8269581056466 57.6958105646630 60.4735883424408 63.2513661202186 66.0291438979964 68.7613843351548 71.4025500910747 73.8615664845173 76.2295081967213 78.5063752276867 80.6921675774135 82.8324225865210 84.8360655737705 86.7941712204007 88.5701275045537 90.2094717668488 91.7577413479053 93.1238615664845 94.2622950819672 95.1730418943534 95.9016393442623 96.7668488160291 96.9489981785064 97.2677595628415 97.4043715846995 97.5409836065574 97.7686703096539 97.9508196721312 97.9508196721312 97.9508196721312 97.9508196721312 97.9508196721312 97.9508196721312 97.9508196721312 97.9508196721312 97.9508196721312 97.9508196721312 97.9508196721312 97.9508196721312 97.9508196721312 ...
97.9508196721312];
time_f9a_185C=[0 3.04279544562230 8.93207695327837 14.8213584609344 19.7290930506478 23.7534354142128 26.9925402434236 29.9371809972517 33.2744405182568 35.7283078131135 39.0655673341186 48.8810365135454 50.7459756576364 52.0219866509619 55.2610914801728 58.7946603847664 61.7393011385944 64.4876325088339 67.6285826462505 71.4566156262269 74.5975657636435 77.7385159010601 81.1739301138594 84.8056537102474 88.4373773066353 92.1672555948174 96.0934432665881 100.019630938359 104.142127993718 108.559089124460 111.209265802905 113.172359638791 115.822536317236 121.122889674126 126.619552414605 131.134668237142 132.312524538673 136.729485669415 142.913231252454 149.391440910876 155.673341185709 156.065959952886 159.206910090302 166.077738515901 172.163329407146 172.948566941500 181.978798586572 189.045936395760 195.916764821358 202.983902630546 204.358068315666 205.241460541814 210.836277974087 217.314487632509 224.087161366313 231.056144483706 237.926972909305 244.601491951315 251.472320376914 258.539458186101 265.606595995289 ...
272.673733804476];
conv9a_185C=[0 0 0 0.409836065573771 3.00546448087432 5.73770491803279 8.56102003642987 11.3387978142077 13.9799635701275 16.0291438979964 18.8524590163934 27.3679417122040 29.5992714025501 31.3296903460838 34.7449908925319 37.5683060109290 40.5737704918033 43.5792349726776 47.1766848816029 49.8178506375228 52.8233151183971 55.8287795992714 58.6520947176685 61.4754098360656 64.2987249544627 67.0309653916211 69.7632058287796 72.4954462659381 75.1366120218579 77.6867030965392 79.0528233151184 80.1457194899818 81.5573770491803 83.6976320582878 85.7468123861567 87.2950819672131 87.7049180327869 89.3442622950820 91.0746812386157 92.3497267759563 92.8506375227687 93.4426229508197 93.8524590163934 94.6721311475410 95.3096539162113 95.4918032786885 96.3114754098361 96.4480874316940 96.9489981785064 97.5409836065574 97.7231329690346 97.4043715846995 98.2240437158470 98.5883424408015 98.6338797814208 98.8160291438980 99.0892531876139 98.8615664845173 98.4517304189435 98.3606557377049 98.4972677595629 ...
98.8160291438980];
time_f9a_190C=[0 0 2.65017667844523 3.92618767177071 10.7970160973695 15.9992147624656 20.8087946603848 24.8331370239497 27.4833137023950 28.7593246957205 30.8205732234001 31.9984295249313 34.7467608951708 34.7467608951708 37.5932469572046 39.9489595602670 42.0102080879466 44.8566941499804 47.7031802120141 48.1939536709855 50.3533568904594 51.1385944248135 53.0035335689046 53.8869257950530 55.6537102473498 56.6352571652925 58.4020416175893 63.9968590498626 65.0765606595995 67.9230467216333 71.1621515508441 74.3031016882607 77.6403612092658 81.6647035728308 84.8056537102474 89.3207695327837 92.7561837455830 103.356890459364 107.872006281900 113.074204946996 118.669022379270 124.460149195132 130.447585394582 136.729485669415 143.207695327837 149.882214369847 156.654888103651 163.525716529250 170.494699646643 177.365528072242 188.064389477817 194.935217903416 208.186101295642 214.762465645858 221.829603455045 228.896741264232 235.865724381625 242.932862190813 250 257.165292500982...
264.134275618375];
conv9a_190C=[0 0 0 0 0 1.22950819672131 3.59744990892532 6.32969034608379 8.69763205828780 9.06193078324226 11.5209471766849 12.0673952641166 15.0728597449909 14.1621129326047 18.0783242258652 21.2204007285975 24.2714025500911 27.3224043715847 30.3734061930783 30.7377049180328 33.4244080145720 34.4717668488160 36.4754098360656 37.5227686703097 39.5264116575592 40.5737704918033 42.5774134790528 48.2240437158470 51.2295081967213 54.2349726775956 57.1493624772313 60.1092896174863 63.0236794171220 65.7103825136612 68.4881602914390 70.9927140255009 73.8615664845173 80.6921675774135 83.1967213114754 85.3825136612022 87.3861566484517 89.2531876138434 91.0291438979964 92.5774134790528 93.8979963570127 95.0364298724955 95.9927140255009 96.7668488160291 97.1311475409836 97.8142076502732 98.4972677595629 98.1785063752277 98.3151183970856 98.7704918032787 98.6794171220401 98.4972677595629 98.5883424408015 98.5428051001822 98.3151183970856 98.1329690346084 ...
98.5428051001822];
% exp data
time_f9a_180C_exp=[0 30.0353356890459 60.0706713780919 90.1060070671378 120.141342756184 150.176678445230 180.212014134276 210.247349823322 240.282685512367 ...
270.318021201413];
conv9a_180C_exp=[0 8.60655737704918 24.5901639344262 55.7377049180328 77.4590163934426 90.9836065573770 95.4918032786885 96.7213114754098 97.5409836065574 ...
98.3606557377049];
time_f9a_185C_exp=[0 30.0353356890459 90.1060070671378 150.176678445230 210.247349823322 ...
270.318021201413];
conv9a_185C_exp=[0 11.8852459016393 68.0327868852459 93.8524590163934 97.5409836065574 ...
98.7704918032787];
time_f9a_190C_exp=[0 30.0353356890459 60.0706713780919 90.1060070671378 120.141342756184 150.176678445230 180.212014134276 210.247349823322 240.282685512367 ...
270.318021201413];
conv9a_190C_exp=[0 11.0655737704918 45.9016393442623 71.3114754098361 87.7049180327869 98.3606557377049 98.7704918032787 99.1803278688525 99.5901639344262 ...
99.5901639344262];
% figure 9b
time_9b_170C=[4.48343079922027 10.9161793372320 18.1286549707602 24.9512670565302 30.7017543859649 37.3294346978557 61.3060428849903 66.7641325536062 72.0272904483431 76.3157894736842 80.7017543859649 84.1130604288499 89.4736842105263 93.8596491228070 97.3684210526316 100.877192982456 105.847953216374 108.869395711501 110.623781676413 113.937621832359 118.713450292398 123.684210526316 128.947368421053 132.456140350877 133.918128654971 137.231968810916 137.816764132554 142.690058479532 147.368421052632 152.534113060429 157.407407407407 162.280701754386 167.641325536062 172.807017543860 174.561403508772 179.922027290448 185.282651072125 190.740740740741 196.296296296296 201.754385964912 207.017543859649 212.280701754386 217.543859649123 221.929824561404 228.070175438596 234.015594541910 240.155945419103 246.198830409357 251.754385964912 257.894736842105 264.035087719298...
270.467836257310];
conv_9b_170C=[1.82149362477232 0.546448087431692 0.318761384335154 1.00182149362477 2.68670309653916 2.64116575591986 9.33515482695810 11.3843351548270 13.5245901639344 15.1639344262295 17.2131147540984 18.7613843351548 20.9016393442623 22.9508196721311 24.5901639344262 26.2295081967213 28.5519125683060 30.3278688524590 30.9198542805100 32.5591985428051 34.9726775956284 37.2950819672131 39.7540983606557 41.3934426229508 42.0765027322404 43.8069216757741 43.9890710382514 46.3114754098361 48.7704918032787 50.9562841530055 53.3242258652095 55.6921675774135 57.7868852459016 59.8360655737705 61.0200364298725 63.1147540983607 65.2094717668488 67.2131147540984 69.2167577413479 71.3114754098361 72.9508196721312 74.5901639344262 76.2295081967213 77.8688524590164 79.5081967213115 81.2386156648452 82.7868852459016 84.4262295081967 85.6557377049180 87.2950819672131 88.8888888888889 ...
90.2094717668488];
time_9b_175C=[3.50877192982456 7.01754385964912 12.2807017543860 14.9122807017544 19.2982456140351 26.2183235867446 33.1384015594542 36.1598440545809 41.2280701754386 46.4912280701754 51.1695906432749 55.6530214424951 61.4035087719298 67.5438596491228 70.1754385964912 76.3157894736842 76.9980506822612 80.7992202729045 85.9649122807018 90.3508771929825 91.0331384015595 95.0292397660819 99.1228070175439 103.508771929825 106.140350877193 109.649122807018 110.916179337232 114.912280701754 117.543859649123 119.298245614035 121.929824561404 126.023391812866 130.019493177388 134.113060428850 138.596491228070 142.982456140351 147.173489278752 151.754385964912 157.017543859649 162.183235867446 163.157894736842 167.543859649123 168.518518518519 171.052631578947 174.561403508772 176.315789473684 180.409356725146 185.964912280702 192.105263157895 198.245614035088 199.025341130604 205.263157894737 211.500974658869 217.836257309942 224.269005847953 230.799220272904 237.426900584795 244.152046783626 251.072124756335 257.894736842105 261.695906432749 264.717348927875 268.518518518519 ...
275.438596491228];
conv_9b_175C=[0.364298724954466 0.409836065573771 0.409836065573771 0.409836065573771 0.546448087431692 1.00182149362477 1.82149362477232 2.59562841530055 4.82695810564663 6.96721311475410 9.38069216757742 11.9307832422587 15.5737704918033 19.6721311475410 21.3114754098361 25.4098360655738 25.7741347905282 28.5519125683060 31.9672131147541 34.4262295081967 35.1092896174863 37.7959927140255 40.4371584699454 43.0327868852459 44.6721311475410 47.1311475409836 48.0874316939891 50.8196721311475 52.5045537340619 53.3697632058288 55.0546448087432 57.6958105646630 60.3825136612022 63.0236794171220 65.5737704918033 68.1238615664845 70.7194899817851 73.1785063752277 75.8196721311475 78.2786885245902 78.6885245901639 80.3734061930783 80.8287795992714 81.9672131147541 83.1967213114754 84.0163934426230 85.3369763205829 87.2950819672131 88.9344262295082 90.5737704918033 90.7103825136612 92.2131147540984 93.7613843351548 95.2185792349727 96.4936247723133 97.6775956284153 98.6794171220401 99.4080145719490 100 100 100 100 100 ...
100];
time_9b_170C_exp=[29.8 59.6 90.4 120 150 180 211 240 270];
conv9b_170C_exp=[1.23 3.28 5.33 19.7 61.5 76.6 79.1 81.1 82.4];
time_9b_175C_exp=[0 29.8 90.4 150 211 270];
conv9b_175C_exp=[0.41 6.97 29.5 79.5 87.3 89.3];
% figure 10
T_f10=[170 174.971751412429 179.971751412429 184.971751412429 190];
time_f10=[103.976377952755 55.0391476489167 29.4092263890742 21.4171004048222 ...
18.4645669291338];
%% PRELIMINARY CALCULATIONS
% General data
T_span=423.15:1:483.15; % [K] vector of temperature to evaluate k
d_paper=0.32299;
tspan=0:300; % [min] integration time span
y0=[1 0 0]; % Initial condition y_S0(0)=1 y_S1(0)=0 X(0)=0
k_paper_span=zeros(1,length(T_span)); % Pre allocating vector of k as a function of T
% Filling the matrix of X_exp, it is possibile to have a different number
% of experimental data depending on the temperature, the matrix is filled
% with zeros in order to produce a square matrix, since at t=0, X=0 the
% calculations are not compromised
t_exp=NaN(5,10); % Pre allocating the matrix of t_exp
t_exp(1,2:end)=time_9b_170C_exp;
t_exp(2,5:end)=time_9b_175C_exp;
t_exp(3,:)=time_f9a_180C_exp;
t_exp(4,5:end)=time_f9a_185C_exp;
t_exp(5,:)=time_f9a_190C_exp;
t_exp=t_exp';
conv_exp=NaN(5,10); % Pre allocating the matrix of X_exp
conv_exp(1,2:end)=conv9b_170C_exp;
conv_exp(2,5:end)=conv9b_175C_exp;
conv_exp(3,:)=conv9a_180C_exp;
conv_exp(4,5:end)=conv9a_185C_exp;
conv_exp(5,:)=conv9a_190C_exp;
conv_exp=conv_exp'/100;
% evaluation of the k given by the paper as a function of temperature
for i=1:length(T_span)
if T_span(i)>=452
k_paper_span(i)=572.5*exp(-5013.2/T_span(i));
else
k_paper_span(i)=2.637*1e9*exp(-11977/T_span(i));
end
end
% finding the values of k for relevent temperatures
indx=[find(T_span==443.15) find(T_span==448.15) find(T_span==453.15) find(T_span==458.15) find(T_span==463.15)];
k_paper=k_paper_span(indx);
T_paper=T_span(indx);
%% DATA FITTING AND NUMERICAL INTEGRATION OF ODEs SYSTEM
% Evaluation of new values for k in order to minimize the error and fit
% the model to the exp data, then the ode system is integrated using the k
% value just found; then the variation of conversion as a function of time
% with different k values is plotted at all the relevant temperatures.
for i=1:length(k_paper)
cc=jet(length(k_paper)); % creating colormap
% Results from paper
[t_paper,y_paper]=ode15s(@(t,y) Reactions(t,y,k_paper(i),d_paper),tspan,y0);
yS0_paper(:,i)=y_paper(:,1);
yS1_paper(:,i)=y_paper(:,2);
X_paper(:,i)=y_paper(:,3);
% Results from model lsqnonlin
[k_model1(i)]=lsqnonlin(@(k) optimization(k,d_paper,t_exp(:,i),conv_exp(:,i)),k_paper(i),[],[]);
[t_model1,y_model1]=ode15s(@(t,y) Reactions(t,y,k_model1(i),d_paper),tspan,y0);
yS0_model1(:,i)=y_model1(:,1);
yS1_model1(:,i)=y_model1(:,2);
X_model1(:,i)=y_model1(:,3);
% Results from ga
[k_model2(i) ]=ga(@(k) optimization_ga(k,d_paper,t_exp(:,i),conv_exp(:,i)),1,[],[],[],[],.000000001,[],[],options_ga);
[t_model2,y_model2]=ode15s(@(t,y) Reactions(t,y,k_model2(i),d_paper),tspan,y0);
yS0_model2(:,i)=y_model2(:,1);
yS1_model2(:,i)=y_model2(:,2);
X_model2(:,i)=y_model2(:,3);
% Results from patternsearch
[k_model3(i) ]=patternsearch(@(k) optimization_ga(k,d_paper,t_exp(:,i),conv_exp(:,i)),k_paper(i),[],[],[],[],[],[],[]);
[t_model3,y_model3]=ode15s(@(t,y) Reactions(t,y,k_model3(i),d_paper),tspan,y0);
yS0_model3(:,i)=y_model3(:,1);
yS1_model3(:,i)=y_model3(:,2);
X_model3(:,i)=y_model3(:,3);
% Creating legendinfo for figure 4
legendinfo1{i}=strcat(['T=' num2str(T_paper(i)-273.15) ' °C']);
% Creating 2 figures, one for low T and one for high T
if i<=2
figure(2)
hold on
plot(t_paper,X_paper(:,i)*100,'LineStyle','-','LineWidth',1.5,'Color',cc(i,:));
plot(t_model1,X_model1(:,i)*100,'LineStyle','--','LineWidth',1.5,'Color',cc(i,:));
plot(t_model2,X_model2(:,i)*100,'LineStyle',':','LineWidth',1.5,'Color',cc(i,:));
plot(t_model3,X_model3(:,i)*100,'LineStyle','-.','LineWidth',1.5,'Color',cc(i,:));
drawnow
xlabel('Time [min]','FontSize',12,'FontWeight','bold')
ylabel('Conversion [%]','FontSize',12,'FontWeight','bold')
title('Conversion of PET-pc as a function of temperature and time. Low T');
else
figure(3)
hold on
plot(t_paper,X_paper(:,i)*100,'LineStyle','-','LineWidth',1.5,'Color',cc(i,:));
plot(t_model1,X_model1(:,i)*100,'LineStyle','--','LineWidth',1.5,'Color',cc(i,:));
plot(t_model2,X_model2(:,i)*100,'LineStyle',':','LineWidth',1.5,'Color',cc(i,:));
plot(t_model3,X_model3(:,i)*100,'LineStyle','-.','LineWidth',1.5,'Color',cc(i,:));
drawnow
xlabel('Time [min]','FontSize',12,'FontWeight','bold')
ylabel('Conversion [%]','FontSize',12,'FontWeight','bold')
title('Conversion of PET-pc as a function of temperature and time. High T');
end
end
% plotting the experimental data
figure(2)
hold on
plot(time_9b_170C_exp,conv9b_170C_exp,'o','Color',cc(1,:),'LineWidth',1.5,'MarkerSize',10)
plot(time_9b_175C_exp,conv9b_175C_exp,'o','Color',cc(2,:),'LineWidth',1.5,'MarkerSize',10)
figure(3)
hold on
plot(time_f9a_180C_exp,conv9a_180C_exp,'o','Color',cc(3,:),'LineWidth',1.5,'MarkerSize',10)
plot(time_f9a_185C_exp,conv9a_185C_exp,'o','Color',cc(4,:),'LineWidth',1.5,'MarkerSize',10,'MarkerEdgeColor','k')
plot(time_f9a_190C_exp,conv9a_190C_exp,'o','Color',cc(5,:),'LineWidth',1.5,'MarkerSize',10)
% plotting a ghost graphical object in order to customize the legend for
% figure 2 and 3
figure(2)
LEGEND1=plot(0,0,'Color',cc(1,:));
hold on
LEGEND2=plot(0,0,'Color',cc(2,:));
legend([LEGEND1 LEGEND2],{'T=170 °C','T=175 °C'},'Location','best','FontSize',12,'FontWeight','bold')
figure(3)
LEGEND3=plot(0,0,'Color',cc(3,:));
hold on
LEGEND4=plot(0,0,'Color',cc(4,:));
LEGEND5=plot(0,0,'Color',cc(5,:));
legend([LEGEND3 LEGEND4 LEGEND5],{'T=180 °C','T=185 °C','T=190 °C'},'Location','best','FontSize',12,'FontWeight','bold')
hold off
%% OPTIMIZATION WITH MORE THAN ONE PARAMETER
% This other model tries to find the best fit by using two parameters: k1
% and k2, k1 is referred to y_S0 and k2 to y_S1
for i=1:length(k_paper)
cc=jet(length(k_paper)); % creating colormap
% Results from paper
[t_paper_new,y_paper_new]=ode15s(@(t,y) Reactions_new(t,y,k_paper(i),d_paper),tspan,y0);
yS0_paper_new(:,i)=y_paper_new(:,1);
yS1_paper_new(:,i)=y_paper_new(:,2);
X_paper_new(:,i)=y_paper_new(:,3);
% Results from model lsqnonlin
[parameters_1]=lsqnonlin(@(par) optimization_new(par,t_exp(:,i),conv_exp(:,i)),[k_paper(i) k_paper(i)],[1 1]*1e-9,[.1 .1]);
k1_model1(i)=parameters_1(1);
k2_model1(i)=parameters_1(2);
[t_model1_new,y_model1_new]=ode15s(@(t,y) Reactions_new(t,y,k1_model1(i),k2_model1(i)),tspan,y0);
yS0_model1_new(:,i)=y_model1_new(:,1);
yS1_model1_new(:,i)=y_model1_new(:,2);
X_model1_new(:,i)=y_model1_new(:,3);
% Results from ga
parameters_2=ga(@(par) optimization_ga_new(par,t_exp(:,i),conv_exp(:,i)),2,[],[],[],[],[1 1]*1e-9,[.1 .1],[],options_ga);
k1_model2(i)=parameters_2(1);
k2_model2(i)=parameters_2(2);
[t_model2_new,y_model2_new]=ode15s(@(t,y) Reactions_new(t,y,k1_model2(i),k2_model2(i)),tspan,y0);
yS0_model2_new(:,i)=y_model2_new(:,1);
yS1_model2_new(:,i)=y_model2_new(:,2);
X_model2_new(:,i)=y_model2_new(:,3);
% Results from patternsearch
parameters_3=patternsearch(@(par) optimization_ga_new(par,t_exp(:,i),conv_exp(:,i)),[k_paper(i) k_paper(i)],[],[],[],[],[1 1]*1e-9,[.1 .1],[]);
k1_model3(i)=parameters_3(1);
k2_model3(i)=parameters_3(2);
[t_model3_new,y_model3_new]=ode15s(@(t,y) Reactions_new(t,y,k1_model3(i),k2_model3(i)),tspan,y0);
yS0_model3_new(:,i)=y_model3_new(:,1);
yS1_model3_new(:,i)=y_model3_new(:,2);
X_model3_new(:,i)=y_model3_new(:,3);
% Creating legendinfo
legendinfo1{i}=strcat(['T=' num2str(T_paper(i)-273.15) ' °C']);
% Creating 2 figures, one for low T and one for high T
if i<=2
figure(5)
hold on
plot(t_paper_new,X_paper_new(:,i)*100,'LineStyle','-','LineWidth',1.5,'Color',cc(i,:));
plot(t_model1_new,X_model1_new(:,i)*100,'LineStyle','--','LineWidth',1.5,'Color',cc(i,:));
plot(t_model2_new,X_model2_new(:,i)*100,'LineStyle',':','LineWidth',1.5,'Color',cc(i,:));
plot(t_model3_new,X_model3_new(:,i)*100,'LineStyle','-.','LineWidth',1.5,'Color',cc(i,:));
drawnow
xlabel('Time [min]','FontSize',12,'FontWeight','bold')
ylabel('Conversion [%]','FontSize',12,'FontWeight','bold')
title('Conversion of PET-pc as a function of temperature and time. Low T');
else
figure(6)
hold on
plot(t_paper_new,X_paper_new(:,i)*100,'LineStyle','-','LineWidth',1.5,'Color',cc(i,:));
plot(t_model1_new,X_model1_new(:,i)*100,'LineStyle','--','LineWidth',1.5,'Color',cc(i,:));
plot(t_model2_new,X_model2_new(:,i)*100,'LineStyle',':','LineWidth',1.5,'Color',cc(i,:));
plot(t_model3_new,X_model3_new(:,i)*100,'LineStyle','-.','LineWidth',1.5,'Color',cc(i,:));
drawnow
xlabel('Time [min]','FontSize',12,'FontWeight','bold')
ylabel('Conversion [%]','FontSize',12,'FontWeight','bold')
title('Conversion of PET-pc as a function of temperature and time. High T');
end
end
% plotting the experimental data
figure(5)
hold on
plot(time_9b_170C_exp,conv9b_170C_exp,'o','Color',cc(1,:),'LineWidth',1.5,'MarkerSize',10)
plot(time_9b_175C_exp,conv9b_175C_exp,'o','Color',cc(2,:),'LineWidth',1.5,'MarkerSize',10)
figure(6)
hold on
plot(time_f9a_180C_exp,conv9a_180C_exp,'o','Color',cc(3,:),'LineWidth',1.5,'MarkerSize',10)
plot(time_f9a_185C_exp,conv9a_185C_exp,'o','Color',cc(4,:),'LineWidth',1.5,'MarkerSize',10,'MarkerEdgeColor','k')
plot(time_f9a_190C_exp,conv9a_190C_exp,'o','Color',cc(5,:),'LineWidth',1.5,'MarkerSize',10)
% plotting a ghost graphical object in order to customize the legend for figures 5 and 6
figure(5)
LEGEND1=plot(0,0,'Color',cc(1,:));
hold on
LEGEND2=plot(0,0,'Color',cc(2,:));
legend([LEGEND1 LEGEND2],{'T=170 °C','T=175 °C'},'Location','best','FontSize',12,'FontWeight','bold')
figure(6)
LEGEND3=plot(0,0,'Color',cc(3,:));
hold on
LEGEND4=plot(0,0,'Color',cc(4,:));
LEGEND5=plot(0,0,'Color',cc(5,:));
legend([LEGEND3 LEGEND4 LEGEND5],{'T=180 °C','T=185 °C','T=190 °C'},'Location','best','FontSize',12,'FontWeight','bold')
function dydt=Reactions_paper(t,y,k,d) % WARNING -->
% this function contains the ODE system given by the paper (with errors)
y_S0=y(1); y_S1=y(2); X=y(3);
dy_S0dt=k*y_S0;
dy_S1dt=4*k*y_S0-0*k*y_S1;
dXdt=(k*y_S0-d*k*y_S1)*(1-X);
dydt=[dy_S0dt dy_S1dt dXdt]';
end
function [dydt]=Reactions(t,y,k,d) % ODE system with corrections, k is determined by the optimization function
y_S0=y(1); y_S1=y(2); X=y(3);
dys0dt=-k*y_S0; % added a - since y_S0 should decrease
dys1dt=4*k*y_S0-1/3*d*k*y_S1; % the consumption term of S1 specie is missing, the stechiometry is unknown, with some trial and error it is possibile to estimate the stoichiometric coefficient for the consumption term
dXdt=-(1/4*k*y_S0-3*d*k*y_S1)*(1-X); % added a -, the derivation of this equation int the paper was wrong, added coefficients to better match the exp data
dydt=[dys0dt dys1dt dXdt];
dydt=dydt';
end
function F=optimization(k,d,t_exp,X_exp) % Optimization function used for function that minimize vector like lsqnonlin
t_exp(isnan(t_exp))=[];
X_exp(isnan(X_exp))=[];
y0=[1 0 0];
[t,y]=ode15s(@(t,y)Reactions(t,y,k,d),0:300,y0);
yS0=y(:,1);
yS1=y(:,2);
X=y(:,3);
for i=1:length(t_exp)
indx(i)=find(abs(t-t_exp(i))<.5); % finding the value of the model that corrisponds to the exp data with .5 tolerance
end
yS0_model=yS0(indx);
yS1_model=yS1(indx);
X_model=X(indx);
err=X_model-X_exp;
F=err;
end
function F=optimization_ga(k,d,t_exp,X_exp) % Optimization function used for function that minimize scalar like ga
t_exp(isnan(t_exp))=[];
X_exp(isnan(X_exp))=[];
y0=[1 0 0];
[t,y]=ode15s(@(t,y)Reactions(t,y,k,d),0:300,y0);
yS0=y(:,1);
yS1=y(:,2);
X=y(:,3);
for i=1:length(t_exp)
indx(i)=find(abs(t-t_exp(i))<.5); % finding the value of the model that corrisponds to the exp data with .5 tolerance
end
yS0_model=yS0(indx);
yS1_model=yS1(indx);
X_model=X(indx);
% evaluating the error, normalizing, deleting NaN elements that can appear
% at the start of the integration
err=(X_model-X_exp)./X_exp;
err=err.^2;
err=err.^.5;
err(isnan(err))=[];
F=sum(err);
end
function [dydt]=Reactions_new(t,y,k1,k2) % ODE system with corrections, k is determined by the optimization function
y_S0=y(1); y_S1=y(2); X=y(3);
dys0dt=-k1*y_S0; % added a - since y_S0 should decrease
dys1dt=4*k1*y_S0-k2*y_S1; % the parameter k2 considers also the stoichometry of the system
dXdt=-(k1*y_S0-k2*y_S1)*(1-X); % added a -, the derivation of this equation int the paper was wrong
dydt=[dys0dt dys1dt dXdt];
dydt=dydt';
end
function F=optimization_new(par,t_exp,X_exp) % Optimization function used for function that minimize vector like lsqnonlin
t_exp(isnan(t_exp))=[];
X_exp(isnan(X_exp))=[];
y0=[1 0 0];
k1=par(1);
k2=par(2);
[t,y]=ode15s(@(t,y)Reactions_new(t,y,k1,k2),0:300,y0);
yS0=y(:,1);
yS1=y(:,2);
X=y(:,3);
for i=1:length(t_exp)
indx(i)=find(abs(t-t_exp(i))<.5); % finding the value of the model that corrisponds to the exp data with .5 tolerance
end
yS0_model=yS0(indx);
yS1_model=yS1(indx);
X_model=X(indx);
err=X_model-X_exp;
F=err;
end
function F=optimization_ga_new(par,t_exp,X_exp) % Optimization function used for function that minimize scalar like ga
t_exp(isnan(t_exp))=[];
X_exp(isnan(X_exp))=[];
y0=[1 0 0];
k1=par(1);
k2=par(2);
[t,y]=ode15s(@(t,y)Reactions_new(t,y,k1,k2),0:300,y0);
yS0=y(:,1);
yS1=y(:,2);
X=y(:,3);
for i=1:length(t_exp)
indx(i)=find(abs(t-t_exp(i))<.5); % finding the value of the model that corrisponds to the exp data with .5 tolerance
end
yS0_model=yS0(indx);
yS1_model=yS1(indx);
X_model=X(indx);
% evaluating the error, normalizing, deleting NaN elements that can appear
% at the start of the integration
err=(X_model-X_exp)./X_exp;
err=err.^2;
err=err.^.5;
err(isnan(err))=[];
F=sum(err);
end