Fitting and parameter estimation of kinetic model ODE system
13 views (last 30 days)
Show older comments
Nadya Shalsabila salman
on 29 Nov 2021
Commented: Star Strider
on 9 Dec 2021
Hello, I am trying to fit experimental data with simultaneous ODE like shown below:
clc
global t c
t = [0
24
48
96
144
192
240
288];
c = [0.08 22.55530474 0.003
0.508650519 22.84424379 0.1349480969
1.089965398 23.13318284 0.9653979239
2.584775087 20.78555305 2.169550173
4.951557093 16.1986456 6.197231834
5.906574394 13.30925508 10.01730104
6.321799308 12.62302483 13.38062284
6.737024221 13.23702032 14.79238754];
theta0 = [0.05;1.66;110;76;3.66;53.26;0.5;0.0128;2.35;0.02;0.0017;63.95;19.1;72.5];
options = optimset('display','iter');
options.MaxFunEvals = 10e8;
options.TolX = 10e-8;
options.TolFun = 10e-7;
options.MaxIter = 10e6;
[theta] = fminsearch(@minfunction, theta0, options);
tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv, c(1,:));
figure
plot(t, c, 'p')
hold on
hlp = plot(tv, Cfit(:,[1 2 3]));
hold off
grid
xlabel('Waktu(jam)')
ylabel('Biomassa(g/L),Substrat(g/L),Astaksantin(mg/L)')
legend(hlp, 'C_1(t)', 'C_2(t)', 'C_3(t)', 'Location','N')
for k1 = 1:length(theta)
fprintf(1, '\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
end
function C=kinetics(theta, t, c0)
[T,Cv]=ode45(@DifEq,t, c0);
function dC=DifEq(t,c)
dcdt=zeros(3,1);
I=80;
mu = ((theta(1).*c(2))./(theta(2)+c(2)+(c(2).^2./(theta(3))))).*((1-c(3)./theta(4)).^theta(5)).*(I/(theta(6)+I));
dcdt(1) = mu.*c(1);
dcdt(2) = -c(1).*((mu./(theta(7))+theta(8)));
dcdt(3) = c(1).*((theta(9).*mu)+theta(10)).*(c(2)/(theta(11)+c(2)+(c(2)^2/theta(12)))).*(c(3)/(theta(13)+c(3))).*(I/(theta(14)+I));
dC=dcdt;
end
C=Cv;
end
function error = minfunction(theta)
global t c
c0 = c(1,:);
c_estim = kinetics(theta,t,c0);
w1 = 1;
w2 = 1;
w3 = 1;
error_temp = w1*abs(c(:,1)-c_estim(:,1)) + ...
w2*abs(c(:,2)-c_estim(:,2)) + ...
w3*abs(c(:,3)-c_estim(:,3));
error = sum(error_temp);
end
However, in that code, I have stated that "I=80", where actually I want to make the value of "I" is 35 (from t 0 to 96 hours) and 80 (t more than 96 hours). "I" is the light intensity.
Could you help me what should I do to change the code according to my request? Thank you for the help!
0 Comments
Accepted Answer
Star Strider
on 29 Nov 2021
Optimising 14 parameters is asking a lot of fminsearch! I switched to lsqcurvefit here because fminsearch is not likely to be useful beyond about 7 parameters. Also, don’t clear at the beginning, and definitely do not use global variables! Pass the other arguments as extra parameters, if necessary.
This works, although it will be necessary to let it run longer than I have here so that lsqcurvefit can converge on a better sollution. I added the initial conditions as the last 3 elements of ‘theta’ so they are now parameters-to-be-estimated. In my experience, this works better than fixing them at specific values in the code.
The equation system can reasonably be called ‘stiff’ so I switched from ode45 toode15s with a significant speed imprivement.
The numeric ODE solvers have problems with abrupt, non-differentialbe ‘step’ transitions, and there are two ways to deal with that. One is to interrupt the integration with the first value of ‘I’ then save the last values and time and re-start it with those as the intial conditions and re-start the integration. The other is to use the tanh function to privide a differentiable transition at the appropriate time. That’s what I did here.
I also updated the plot calls so that the colours match and the plot makes more sense.
% clc
% global t c % Don't Use 'global' Variables!
t = [0
24
48
96
144
192
240
288];
c = [0.08 22.55530474 0.003
0.508650519 22.84424379 0.1349480969
1.089965398 23.13318284 0.9653979239
2.584775087 20.78555305 2.169550173
4.951557093 16.1986456 6.197231834
5.906574394 13.30925508 10.01730104
6.321799308 12.62302483 13.38062284
6.737024221 13.23702032 14.79238754];
theta0 = [0.05;1.66;110;76;3.66;53.26;0.5;0.0128;2.35;0.02;0.0017;63.95;19.1;72.5];
theta0 = [theta0.*rand(size(theta0)); c(1,:).'];
% theta0 = rand(17,1);
% % options = optimset('display','iter');
% % options.MaxFunEvals = 10e8;
% % options.TolX = 10e-8;
% % options.TolFun = 10e-7;
% % options.MaxIter = 10e6;
[theta] = lsqcurvefit(@objfunction, theta0, t, c, zeros(size(theta0)))
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, compose('C_%d',1:size(c,2)), 'Location','best')
for k1 = 1:length(theta)
fprintf(1, '\t\tTheta(%2d) = %10.6f\n', k1, theta(k1))
end
function C=kinetics(theta, t)
c0 = theta(end-2:end); % Last Three 'theta' Elements Are Initial Conditions (To Be Estimated As Parameters)
I = @(t) 35+22.5*(1+tanh(t-96)); % Differentiable Function For 'I' That 'Steps' At The Appropriate Time!
[T,Cv]=ode15s(@DifEq,t, c0);
function dC=DifEq(t,c)
dcdt=zeros(3,1);
mu = ((theta(1).*c(2))./(theta(2)+c(2)+(c(2).^2./(theta(3))))).*((1-c(3)./theta(4)).^theta(5)).*(I(t)/(theta(6)+I(t)));
dcdt(1) = mu.*c(1);
dcdt(2) = -c(1).*((mu./(theta(7))+theta(8)));
dcdt(3) = c(1).*((theta(9).*mu)+theta(10)).*(c(2)/(theta(11)+c(2)+(c(2)^2/theta(12)))).*(c(3)/(theta(13)+c(3))).*(I(t)/(theta(14)+I(t)));
dC=dcdt;
end
C=Cv;
end
% function error = minfunction(theta,t,c)
% % global t c % Don't Use 'global' Variables!
% c0 = c(1,:);
% c_estim = kinetics(theta,t,c0);
% w1 = 1;
% w2 = 1;
% w3 = 1;
% error_temp = w1*abs(c(:,1)-c_estim(:,1)) + ...
% w2*abs(c(:,2)-c_estim(:,2)) + ...
% w3*abs(c(:,3)-c_estim(:,3));
% error = sum(error_temp);
function c_estim = objfunction(theta,t)
% global t c % Don't Use 'global' Variables!
% c0 = c(1,:);
c_estim = kinetics(theta,t);
% w1 = 1;
% w2 = 1;
% w3 = 1;
% error_temp = w1*abs(c(:,1)-c_estim(:,1)) + ...
% w2*abs(c(:,2)-c_estim(:,2)) + ...
% w3*abs(c(:,3)-c_estim(:,3));
% error = sum(error_temp);
end
Let it run a bit longer to get a better fit by increasing the funciton evaluation limit and other options as necessary. Another option is to sue the genetic algorithm to estimate the parameters, however it is likely that lsqcurvefit can estimate them if given a longer time to do so.
.
10 Comments
More Answers (0)
See Also
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!