Fitting kinetic model to estimate kinetic parameter

115 次查看（过去 30 天）
Nadya Shalsabila salman
评论： Star Strider 2024-11-24，19:49
Hello, I want to ask a question. I want to curve fit the data with the kinetic model. I am able to run the code (written and attached below), but the result is still far from my expectation.
function trycurvefitlunaflores
function C=kinetics(theta,t)
c0=[0.6178; 28.90; 0.4589];
[T,Cv]=ode45(@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));
dcdt(1) = mu.*c(1);
dcdt(2) = -c(1).*((mu./(theta(6))+theta(7)));
dcdt(3) = c(1).*((theta(8).*mu)+theta(9));
dC=dcdt;
end
C=Cv;
end
t = [0
12
24
36
48
60
72];
c = [0.6178489703 28.90160183 0.7586206897
0.823798627 28.83295195 4.045977011
2.402745995 21.62471396 6.827586207
5.972540046 13.04347826 18.5862069
8.306636156 6.178489703 34.01149425
9.885583524 2.402745995 44.88505747
9.542334096 2.059496568 50.44827586];
theta0 = [1;1;1;1;1;1;1;1;1];
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);
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(1)
plot(t, c, 'p')
hold on
hlp = plot(tv, Cfit);
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, 'C_1(t)', 'C_2(t)', 'C_3(t)', 'Location','N')
end
In the normal circumstance, the plot line (from kinetic model) should follow the data's tren, like this figure below.
Could you help me to fix the plot to make it more similar to the data? Any little advice or help will be much appreciated. Thank you.
  7 个评论
显示 5更早的评论
Nadya Shalsabila salman
Nadya Shalsabila salman 2021-11-8
Is it possible that this difference is due to different units of the two plots? because if I get a suitable plot, the parameter value is very far from what it should be
the substrate, product, and biomass data should be as follows
but in the plotting, I used the initial product (karotenoid) = 0.75, not 75.86 mu g/L
and the parameter should be like this
and ms 0.005 g g-1 h-1
I'm still confused about this unit difference

请先登录，再进行评论。

请先登录，再回答此问题。

采纳的回答

Alex
Alex 2021-11-12
编辑：Alex 2021-11-16
Using fminsearch, defining a custom cost-function and using the initial values of Alex Sha finally did the trick.
I seperated the program into seperate files. The main file:
clear
close all
clc
global t c
t = [0
12
24
36
48
60
72];
c = [0.6178489703 28.90160183 0.7586206897
0.823798627 28.83295195 4.045977011
2.402745995 21.62471396 6.827586207
5.972540046 13.04347826 18.5862069
8.306636156 6.178489703 34.01149425
9.885583524 2.402745995 44.88505747
9.542334096 2.059496568 50.44827586];
theta0 = [-0.024;-12.55;-28.80;50.63;0.168;0.3501;0.003146;1.9802;0.0954];
options = optimset('display','iter');
options.MaxFunEvals = 10e8;
options.TolX = 10e-5;
options.TolFun = 10e-5;
options.MaxIter = 10e3;
[theta] = fminsearch(@minfunction, theta0, options);
tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv, c(1,:));
figure(1)
plot(t, c, 'p')
hold on
hlp = plot(tv, Cfit);
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, 'C_1(t)', 'C_2(t)', 'C_3(t)', 'Location','N')
the kinetics function:
function C=kinetics(theta, t, c0)
[T,Cv]=ode45(@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));
dcdt(1) = mu.*c(1);
dcdt(2) = -c(1).*((mu./(theta(6))+theta(7)));
dcdt(3) = c(1).*((theta(8).*mu)+theta(9));
dC=dcdt;
end
C=Cv;
end
and the cost function:
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
In the cost function one could use different weights for the error of c1, c2 and c3. This could be of interest if one would worry more about the error for example of c1 than c2.
fminsearch now tries to find parameters for the differential equation which minimizes:
Hope this helps.
  6 个评论
显示 4更早的评论
Alex
Alex 2021-11-16
编辑：Alex 2021-11-16
The idea behind the weights was to give you the option to say: "i worry more about the error of c3, than c2 and therefore want to give the error a higher weight".
Looking again at the cost function makes me think about an error I made there. Maybe it is better to define the cost function like:
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
The error is now defined the following way:
(Before it would have been possible that errors cancel out each other, which could have been the reason why I had to add the weights. This is of course not the prefered way.). As you can see I now have set the weights to 1, since it looks as if they are not really necessary anymore. But if you want to play around with the weights, you are of course free to go if it improves the results.
I will update the original answer.
Nadya Shalsabila salman
Nadya Shalsabila salman 2021-11-16
Okay.. thank you very much @Alex. Hereby I attach the latest code I have tried, this code is using the latest error code you provided and the result is pretty decent.

请先登录，再进行评论。

更多回答（2 个）

Alex
Alex 2021-11-5
编辑：Alex 2021-11-5
Your basic structure seems to be ok. Here are some hints you could consider:
  • Try playing around with the solver options, like using different algorithms, different initial points and tolerances:
theta0 = [1;1;1;1;1;1;1;1;1];
options = optimoptions('lsqcurvefit', 'Algorithm','trust-region-reflective');
options.MaxFunctionEvaluations = 10e5;
options.StepTolerance = 10e-7;
options.FunctionTolerance = 10e-7;
options.MaxIterations = 10e5;
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c, [],[], options);
  • Use the accurate intitial condition for your differential equation:
c0=[0.6178489703 28.90160183 0.7586206897];
  • Watch out for local minimum warnings. Matlab provides further thoughts on how to avoid local minimum convergence when clicking on the link in the matlab console (or directly via this link):
>> trycurvefitlunaflores
Local minimum possible.
After playing around a bit I was able to produce the following result:
I am sure that with a little tweaking you will be able to also fit C_1 correctly. Let me know if you found a solution.
  9 个评论
显示 7更早的评论
Nadya Shalsabila salman
Nadya Shalsabila salman 2021-11-13
Hi @Alex Sha, thank you very much for your answer! Anyway, I want to try 1stOpt Software to estimate the parameter. But, when I request the link, the software informer tell me "The download link for requested software cannot be found. We will update our database soon and let you know when 1stOpt is available for download"
Can you tell me how to download the 1st Opt?

请先登录，再进行评论。

Manish
Manish 2024-11-24，17:42
function C=kinetics(theta,t)
c0=[0.6178; 28.90; 0.4589];
[T,Cv]=ode45(@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));
dcdt(1) = mu.*c(1);
dcdt(2) = -c(1).*((mu./(theta(6))+theta(7)));
dcdt(3) = c(1).*((theta(8).*mu)+theta(9));
dC=dcdt;
end
C=Cv;
end
t = [0
12
24
36
48
60
72];
c = [0.6178489703 28.90160183 0.7586206897
0.823798627 28.83295195 4.045977011
2.402745995 21.62471396 6.827586207
5.972540046 13.04347826 18.5862069
8.306636156 6.178489703 34.01149425
9.885583524 2.402745995 44.88505747
9.542334096 2.059496568 50.44827586];
theta0 = [1;1;1;1;1;1;1;1;1];
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);
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(1)
plot(t, c, 'p')
hold on
hlp = plot(tv, Cfit);
hold off
grid
xlabel('Time')
ylabel('Concentration')
legend(hlp, 'C_1(t)', 'C_2(t)', 'C_3(t)', 'Location','N')
end

类别

在 Help Center 和 File Exchange 中查找有关 Eigenvalue Problems 的更多信息

标签

Community Treasure Hunt

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

Start Hunting!

Translated by