Multistart - number of iterations
2 Ansichten (letzte 30 Tage)
Ältere Kommentare anzeigen
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
%Data points%
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]
%multistart optimisation%
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);
%objective function%
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
1 Kommentar
Alan Weiss
am 29 Jan. 2020
I am not sure about what is going on with your optimization, but do you have an error in this line of code toward the end of DifEq:
dcdt(1)=0;
This overrides the definition toward the beginning of the function
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));
which makes more sense as a reaction rate.
Alan Weiss
MATLAB mathematical toolbox documentation
Antworten (0)
Siehe auch
Kategorien
Mehr zu Financial Toolbox finden Sie in Help Center und File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!