Parameter Estimation of 3 ODE systems with lsqcurvefit does not converge

Question

gina_et_berg 2020-4-7

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/516120-parameter-estimation-of-3-ode-systems-with-lsqcurvefit-does-not-converge

评论： Star Strider 2020-4-8

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 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

请先登录，再回答此问题。

Answer 1

Star Strider 2020-4-7

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/516120-parameter-estimation-of-3-ode-systems-with-lsqcurvefit-does-not-converge#answer_424626

在 MATLAB Online 中打开

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 个评论
显示 6更早的评论隐藏 6更早的评论

Star Strider 2020-4-8

在 MATLAB Online 中打开

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.

gina_et_berg 2020-4-8