Fitting and parameter estimation of kinetic model ODE system

16 views (last 30 days)
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!

Accepted Answer

Star Strider
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)))
Solver stopped prematurely. lsqcurvefit stopped because it exceeded the function evaluation limit, options.MaxFunctionEvaluations = 1.700000e+03.
theta = 17×1
0.0329 0.6924 14.6130 42.2253 3.4208 0.6521 0.2545 0.0003 0.0337 0.3789
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
Theta( 1) = 0.032858 Theta( 2) = 0.692390 Theta( 3) = 14.612964 Theta( 4) = 42.225345 Theta( 5) = 3.420782 Theta( 6) = 0.652061 Theta( 7) = 0.254480 Theta( 8) = 0.000271 Theta( 9) = 0.033661 Theta(10) = 0.378941 Theta(11) = 49.187783 Theta(12) = 11.247759 Theta(13) = 0.845725 Theta(14) = 13.764339 Theta(15) = 0.216072 Theta(16) = 22.779756 Theta(17) = 0.189724
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

Sign in to comment.

More Answers (0)

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!