MATLAB Answers

Parameter Estimation of 3 ODE systems with lsqcurvefit does not converge

11 views (last 30 days)
Hello
I have atttached my codes. I have used similar codes as Star Strider 's Monod kinetics and curve fitting techniques. But it does not converge. Can you please help me figure out what I am doing wrong here?
Thanks

  0 Comments

Sign in to comment.

Accepted Answer

Star Strider
Star Strider on 7 Apr 2020
First, re-define ‘alpha0’ to include the initial conditions::
alpha0 = [alpha0; 80.2;4.451e-3;4.264e-3];
and change ‘kinetics’ correspondingly so that:
c0 = alpha(end-2:end);
With those changes, it converged very quickly for me (R2020a) producing:
Local minimum possible.
lsqcurvefit stopped because the final change in the sum of squares relative to
its initial value is less than the value of the function tolerance.
<stopping criteria details>
kinetic parameters:
alpha(1) = 1160.81907
alpha(2) = 0.60175
alpha(3) = 2.38292
alpha(4) = 1.71922
alpha(5) = 1482.53287
alpha(6) = 0.21982
alpha(7) = 2.88669
alpha(8) = 2.53375
alpha(9) = 77.58045
alpha(10) = 2.62417
alpha(11) = 2.62398
If it does not produce the parameter values you expect, it would be best to use a global search, such as the genetic algorithm (ga) function. If you have the Global Optimization Toolbox, I will post a version of that code that uses ga.

  8 Comments

Show 5 older comments
Star Strider
Star Strider on 8 Apr 2020
Unfortunately, ga is not able to do anything with your system. I let it run for several hours this afternoon with no results at all. I would just go with lsqcurvefit.
The ga code runs as expected with the original code, however not with your differential equation system. With your system, it displays a few generations, then stalls, and even including 'StallTimeLimit' in the options structure does not solve the problem. It does not throw errors or fail to integrate, since it calls the ‘kinetics’ function appropriately and it produces the appropriate output, however it completely fails to produce any results.
The code that I ran (for you to experiment with) is:
function ParaFun
function C=kinetics(alpha,t)
c0=[80.2;4.451e-3;4.264e-3];
% c0 = alpha(end-2:end);
[T,Cv]=ode45(@DifEq,t,c0);
%
function dC=DifEq(t,c)
dcdt=zeros(3,1);
dcdt(1)=-alpha(1).*alpha(2).* (c(2)^alpha(3))*(c(1)^alpha(4))-...
alpha(5).*alpha(6).*(c(3)^alpha(7))*(c(1)^alpha(8));
dcdt(2)= -alpha(1).*alpha(2).* (c(2)^alpha(3))*(c(1)^alpha(4));
dcdt(3)= -alpha(5).*alpha(6).*(c(3)^alpha(7))*(c(1)^alpha(8));
dC=dcdt;
end
C=Cv;
end
t_obs = [0 1 2 3 4 5 6 7 11 14 21 35 56];
c_obs = [80.2000 0.0045 0.0043
78.5000 0.0026 0.0026
77.6700 0.0022 0.0021
75.6000 0.0019 0.0018
74.2000 0.0016 0.0017
72.2000 0.0015 0.0016
71.1000 0.0014 0.0014
71.4000 0.0013 0.0013
70.5000 0.0012 0.0012
69.5000 0.0010 0.0011
69.0000 0.0009 0.0010
68.0000 0.0008 0.0009
67.7000 0.0007 0.0008];
alpha0 = [1.18E+03;7.00E-01;2.39E+00;1.70E+00;1.50E+03;1.80E-01;2.89E+00;2.59E+00];
c = c_obs;
t = t_obs;
ftns = @(theta) norm(c-kinetics(theta,t));
PopSz = 500;
Parms = 8;
opts = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',(randi(1E+4,PopSz,Parms)*1E-3)+alpha0', 'MaxGenerations',2E3, 'PlotFcn',@gaplotbestf, 'PlotInterval',1, 'StallTimeLimit',30);
t0 = clock;
fprintf('\nStart Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t0)
[theta,fval,exitflag,output] = ga(ftns, Parms, [],[],[],[],zeros(Parms,1),Inf(Parms,1),[],[],opts)
t1 = clock;
fprintf('\nStop Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t1)
GA_Time = etime(t1,t0)
hmsdv = seconds(GA_Time);
hmsdv.Format = 'hh:mm:ss'
fprintf('\nElapsed Time: %23.15E\t\t%s\n', GA_Time, hmsdv)
fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta)
fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
end
tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);
figure
hd = plot(t, c, 'p');
for k1 = 1:size(c,2)
CV(k1,:) = hd(k1).Color;
hd(k1).MarkerFaceColor = CV(k1,:);
end
hold on
hlp = plot(tv, Cfit);
for k1 = 1:size(c,2)
hlp(k1).Color = CV(k1,:);
end
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, 'C_1(t)', 'C_2(t)', 'C_3(t)', 'C_4(t)', 'Location','N')
end
I experimented with several different approaches for the 'InitialPopulationMatrix', none of whiich worked. I do not see any errors in any of the code.
I could possibly let it run for a few more hours, however I doubt that would be worthwhile. If I do and if it produces meaningful results, I will post them back here as an edit to this Comment.

Sign in to comment.

More Answers (1)

gina_et_berg
gina_et_berg on 8 Apr 2020
That was a huge help. Now at least I know whatever I used manually is a better fit to my data. Thanks a lot.

Sign in to answer this question.