Hi,
I've got a multistart optimisation set up to find optimal rate constants for a kinetic model. This works fine, however I wanted to check whether increasing the number of iterations would improve the result. In short, it doesn't seem to, but it also seems to make little sense to me - I actually get worse results? Maybe I'm being really naive, but I'm not quite following why I shouldn't at least get the same result (i.e. if I increase no. of iterations to 50 from 30, shouldn't I still at least get the same best minimum as for 30?)? If anyone has any insight, I'd be all ears. This isn't exactly a critical thing, but it's sort of bugging me!
I'm using the same rng stream each time. Should also explain - I initially ran it with 30 iterations, and saved the rng stream from that. If I re-run it with 30 I do always get the same result.
Below are the scripts for: 1) the data points; 2) the multistart optimisation; 3) the objective function. I've also attached the rng stream.
Any insight welcome!
Thanks a lot,
Matt
time_full = [0;1;2;3;4;5;6;7;8;9;10;24;48;96;144;192;240;288;336;384;432];
label_full = [95.898177 0 0 0 0 4 0.1 0.001055 0.000768
95.898177 0.116654171 0.000169089 0 0 3.883345829 0.099830911 0.001055 0.000768
95.898177 0.221311714 0.000215784 0 0 3.778688286 0.099784216 0.001055 0.000768
95.898177 0.361815956 0.001337332 0 0 3.638184044 0.098662668 0.001055 0.000768
95.898177 0.443043182 0.001897635 0 0 3.556956818 0.098102365 0.001055 0.000768
95.898177 0.511783075 0.002904908 2.73612E-06 0 3.488216925 0.097095092 0.001052264 0.000768
95.898177 0.596086415 0.003847949 2.19614E-05 0 3.403913585 0.096152051 0.001033039 0.000768
95.898177 0.70422927 0.004991988 3.41544E-05 0 3.29577073 0.095008012 0.001020846 0.000768
95.898177 0.736165548 0.005402772 4.34905E-05 0 3.263834452 0.094597228 0.00101151 0.000768
95.898177 0.863634039 0.006882534 4.92657E-05 0 3.136365961 0.093117466 0.001005734 0.000768
95.898177 0.961691531 0.007922809 8.34772E-05 0 3.038308469 0.092077191 0.000971523 0.000768
95.898177 1.694849679 0.02488041 0.000242074 0.000298399 2.305150321 0.07511959 0.000812926 0.000469601
95.898177 2.156329216 0.046290117 0.00048092 0.00018762 1.843670784 0.053709883 0.00057408 0.00058038
95.898177 2.28066375 0.058965709 0.000558424 0.000402402 1.71933625 0.041034291 0.000496576 0.000365598
95.898177 2.263872217 0.060281286 0.000637646 0.000436729 1.736127783 0.039718714 0.000417354 0.000331271
95.898177 2.280749095 0.059985812 0.000590469 0.000523726 1.719250905 0.040014188 0.000464531 0.000244274
95.898177 2.250775532 0.059775064 0.000604293 0.00049205 1.749224468 0.040224936 0.000450707 0.00027595
95.898177 2.248860841 0.059972783 0.000622058 0.000422935 1.751139159 0.040027217 0.000432942 0.000345065
95.898177 2.28310359 0.059096809 0.000547755 0.000441093 1.71689641 0.040903191 0.000507245 0.000326907
95.898177 2.324286746 0.058745232 0.000581056 0.000433505 1.675713254 0.041254768 0.000473944 0.000334495
95.898177 2.298780141 0.059541016 0.000586692 0.000465395 1.701219859 0.040458984 0.000468308 0.000302605]
rng(stream);
problem = createOptimProblem ('lsqcurvefit','objective', @kinetics_full, 'x0', randn(1,7), ...
'lb', [0;0;0;0;0;0;0.03], 'ub', [10;10;10;10;10;10;0.08],'xdata',time_full,'ydata',label_full);
ms = MultiStart('PlotFcns',@gsplotbestf);
[k,fval] = run(ms,problem,30);
function C=kinetics_full(k,t)
c0=[95.898177 0 0 0 0 4 0.1 0.001055 0.000768];
[T,C]=ode45(@(t,c) DifEq(t,c,k),t,c0);
function dC=DifEq(~,c,k)
dcdt=zeros(9,1);
dcdt(1)=-k(1).*c(1)+k(7).*(c(2)+c(3)+c(4)+c(5)+c(6)+c(7)+c(8)+c(9))+k(5).*(c(4)+c(8))+k(6).*(c(5)+c(9));
dcdt(2)=0.55*k(1).*c(1)-k(2).*c(2)-k(7).*c(2);
dcdt(3)=k(2).*c(2)-k(3).*c(3)-k(7).*c(3);
dcdt(4)=k(3).*c(3)-k(4).*c(4)-k(5).*c(4)-k(7).*c(4);
dcdt(5)=k(4).*c(4)-k(6).*c(5)-k(7).*c(5);
dcdt(6)=0.45*k(1).*c(1)-k(2).*c(6)-k(7).*c(6);
dcdt(7)=k(2).*c(6)-k(3).*c(7)-k(7).*c(7);
dcdt(8)=k(3).*c(7)-k(4).*c(8)-k(5).*c(8)-k(7).*c(8);
dcdt(9)=k(4).*c(8)-k(6).*c(9)-k(7).*c(9);
dcdt(1)=0;
dC=dcdt;
end
end
