Parameter Estimation of 3 ODE systems with lsqcurvefit does not converge

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

Answer 1

[Star Strider]

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.

[gina_et_berg]

Thanks for your reply and corrections. I have a problem here through. I have my initial concentrations (aka c0). So maybe that is the problem. Do I have to change ‘kinetics’ ode function? Thanks

[Star Strider]

My pleasure.

There is no need to change ‘kinetics’ other than to re-define ‘c0’ and ‘alpha0’ as I already posted, so that the last three elements of ‘alpha0’ are ‘c0’, and ‘c0’ is re-defined in ‘kinetics’ to reflect that. A common problem with ODEs in parameter estimation problems failing to converge is that the initial conditions are fixed. If you let them become parameters to be estimated, the regression usually works.

[gina_et_berg]

I just check and I have Global Optimization Toolbox.

[gina_et_berg]

In that case my inital values for two elemens (alpha10 and 11) are 3 order of magnitute bigger than actual experimetal values.

[gina_et_berg]

Can you please post the (ga) version?

[Star Strider]

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]

thanks a million! really appreciate your help

[Star Strider]

As always, my pleasure!

I wish I could have gotten it to work.

.

Answer 2

[gina_et_berg]

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

[Star Strider]

As always, my pleasure!