Fitting and parameter estimation of kinetic model ODE system
Show older comments
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
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
Nadya Shalsabila salman
on 30 Nov 2021
My pleasure!
As I mentioned in my original Answer:
‘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 provide a differentiable transition at the appropriate time. That’s what I did here.’
So I created the ‘I’ function to provide this solution —
t = linspace(0, 288, 500);
I = @(t) 35+22.5*(1+tanh(t-96));
figure
plot(t,I(t))
grid
xlabel('t')
ylabel('I(t)')
xlim([0 288])
ylim([30 90])
So ‘I’ is 35 until the function gets to 93 then shifts at 96 to become 80 after that. The time at which it shifts can be adjusted by changing the 96 in the tanh argument to something else. Everything else about the function remains the same.
For example, changing the function slightly to add an additional argument —
t = linspace(0, 288, 500);
I = @(t,tau) 35+22.5*(1+tanh(t-tau));
figure
plot(t,I(t,90),'-b')
hold on
plot(t, I(t,96),'-g')
plot(t, I(t,102),'-r')
hold off
grid
xlabel('t')
ylabel('I(t)')
legend('\tau = 90','\tau = 96','\tau = 102', 'Location','best')
xlim([0 288])
ylim([30 90])
Choose whatever value of τ works best (that is, wherever the transition should be, either before (
), at (
), or after (
) 96 hours). The differentiable function provides a smooth transition, so it should work regardless of the value of τ, providing that the value is within reason and that it meets the criteria set out in the assumptions of the differential equation system.
), at (), or after () 96 hours). The differentiable function provides a smooth transition, so it should work regardless of the value of τ, providing that the value is within reason and that it meets the criteria set out in the assumptions of the differential equation system. .
Nadya Shalsabila salman
on 30 Nov 2021
Star Strider
on 30 Nov 2021
As always, my pleasure!
.
Nadya Shalsabila salman
on 30 Nov 2021
Star Strider
on 30 Nov 2021
No worries.
Use this version of ‘I’ —
I = @(t,tau,lo,hi) lo+(hi-lo)/2*(1+tanh(t-tau)); % Arguments: 't'=Time Vector, 'tau'=Delay, 'lo'=Low Value (Before 'tau'), 'hi'=High Value (After 'tau')
Copy and paste the entire line (including the comment) so that the comment-documentation stays with it. It now permits varying all three parameters, the delay (switching time), low value, and high value. I tested it and it appears to work for every real value of every parameter, regardless of sign.
.
Nadya Shalsabila salman
on 9 Dec 2021
I got it to plot correctly, however I had to abandon the custom colours and the legend, partially because that tends to confuse the interpretation of the yyaxis plot. (The key to making it work with yyaxis is completely copying the plotting loops under both yyaxis calls.)
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
yyaxis left
hdl = plot(t, c(:,[1 2]), 'p');
% for k1 = 1:size(c,2)
% CV(k1,:) = hd(k1).Color;
% hd(k1).MarkerFaceColor = CV(k1,:);
% end
hold on
hlpl = plot(tv, Cfit(:,[1 2]));
% for k1 = 1:size(c,2)
% hlp(k1).Color = CV(k1,:);
% end
hold off
ylabel('Concentration C_1, C_2')
yyaxis right
hdr = plot(t, c(:,3), 'p');
% for k1 = 1:size(c,2)
% CV(k1,:) = hd(k1).Color;
% hd(k1).MarkerFaceColor = CV(k1,:);
% end
hold on
hlpr = plot(tv, Cfit(:,3));
% for k1 = 1:size(c,2)
% hlp(k1).Color = CV(k1,:);
% end
hold off
grid
xlabel('Time')
ylabel('Concentration C_3')
% 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
.
Nadya Shalsabila salman
on 9 Dec 2021
Star Strider
on 9 Dec 2021
As always, my pleasure!
.
More Answers (0)
Categories
Find more on Genetic Algorithm in Help Center and File Exchange
on 29 Nov 2021
on 9 Dec 2021
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
