Coefficient estimation for a system of coupled ODEs

Hello to you all, I'm a complete beginner to Matlab and I'm trying to solve a system of 4 ODEs that have unknown parameters. My intention is to estimate the parameters to obtain a solution then use experimental data to optimize the solution of the ODE system.

I modiffied the code provided by user Star Strider in the question: "Parameter Estimation for a System of Differential Equations" to use the ODEs and data I posses since the two problems are very similar. However the resulting fit the program found indicated to me that I am doing something wrong.

I will attach the .m with my code, inside is also contained the experimental data.

For reference the ODE system I'm trying to solve is the one displayed on the attached image.

[Wan Ji]

@Star Strider

Answer 1

[Star Strider]

I do not see that you are doing anything wrong, and the code runs without error.

The changes I made were to estimate the initial conditions as well as the other parameters, and use random values for ‘theta0’. I expanded the format for the ‘Constants:’ output in order to see the full precision, and changed the plot loop so that the fitted lines are the same colours ad the data points. (I made this change in other posts on this topic, and so am updating them here as well.) The fit appears to be reasonably good with these changes.

The notice is simply that lsqcurvefit exceeded the function evaluation limit before it was convinced that the model converged. You can increase the function evaluation limit (see the documentation section on options), however that may not be necessary. (I increased it to be 15000 and the iteration limit to 1500, and the fit appears to be even better, and the estimation converges.) Use the estimated parameters here for ‘theta0’ if you want to run it with changed options, since this appears to be reasonably good approximation.

modelacion

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.

	Constants:
		Theta( 1) =     1.225621526046
		Theta( 2) =     8.818523145614
		Theta( 3) =     0.832213488120
		Theta( 4) =  8864.918733082799
		Theta( 5) =     1.234266936058
		Theta( 6) =     0.810227130705
		Theta( 7) =     0.157698146317
		Theta( 8) =     0.023932541443
		Theta( 9) =     0.089951988393
		Theta(10) =    36.855992659497
		Theta(11) =     0.565552316098
		Theta(12) =     0.000134265635

function modelacion

function C=kinetics(theta,t)
% c0=[1;0;0;0];
c0 = theta(9:12);
[T,Cv]=ode45(@DifEq,t,c0);
%
    function dC=DifEq(t,c)
    dcdt=zeros(4,1);
    dcdt(1)= (theta(1).*(1-(c(1)./theta(2)))-theta(3)).*c(1);
    dcdt(2)= -((1./theta(4)).*dcdt(1))-(theta(5).*c(1));
    dcdt(3)= theta(6).*dcdt(1)+theta(7).*c(1);
    dcdt(4)= theta(8).*dcdt(1);
    dC=dcdt;
    end
C=Cv;
end



t=[0
0.928
1.984
2.992
3.952
5.008
5.968
6.976
7.98400000000
8.992
9.952
11.008];

c=[0.075	36.77419355	0.317906336	0
0.075	36.77419355	0.545179063	0.000434783
0.175	36.4516129	0.71184573	0.001304348
0.325	36.12903226	0.909090909	0.003478261
0.425	36.12903226	1.106060606	0.006086957
0.525	35.48387097	1.242424242	0.010869565
0.675	34.19354839	1.439393939	0.013913043
0.925	33.5483871	1.636363636	0.020869565
1.15	32.25806452	1.985123967	0.026521739
1.475	30.32258065	2.454820937	0.032608696
1.8	28.38709677	3.015702479	0.040869565
2.075	26.12903226	3.54600551	0.047826087]; 
% theta0=[1;1;1;1;1;1;1;1];
theta0 = rand(12,1);

opts = optimoptions('lsqcurvefit', 'MaxFunctionEvaluations', 15000, 'MaxIterations',1500);
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c,zeros(size(theta0)),Inf(size(theta0)), opts);

fprintf(1,'\tConstants:\n')
for k1 = 1:length(theta)
    fprintf(1, '\t\tTheta(%2d) = %18.12f\n', k1, theta(k1))
end

tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);

figure(1)
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, compose('C_{%d}(t)',1:size(c,2)), 'Location','best')
end

Make appropriate changes to get the result you want.

.

[pablo cardona]

Your help is very much appreciated!

[Star Strider]

As always, my pleasure!

.

[Alex Sha]

Hi, Pablo, according to the picture of ODE equations you posted, as above, in your code:

dcdt(2)= -((1./theta(4)).*dcdt(1))-(theta(5).*c(1));

be should:

dcdt(2)= -((1./theta(4)).*dcdt(3))-(theta(5).*c(1));

?