Fitting kinetic model to estimate kinetic parameter

Question

Nadya Shalsabila salman 2021-11-4

1
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1579144-fitting-kinetic-model-to-estimate-kinetic-parameter

编辑： Alex 2021-11-16

trycurvefitlunaflores.m

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.

6 个评论
显示 4更早的评论隐藏 4更早的评论

Nadya Shalsabila salman 2021-11-8

Hi @Star Strider! Thanks for trying the problem I had. if you don't mind, may i see the code you have tried?

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

请先登录，再进行评论。

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

Answer 1

Alex 2021-11-12

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1579144-fitting-kinetic-model-to-estimate-kinetic-parameter#answer_829984

编辑：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更早的评论隐藏 4更早的评论

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 2021-11-16

terbaruLuna.m

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.

请先登录，再进行评论。

Answer 2

Alex 2021-11-5

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1579144-fitting-kinetic-model-to-estimate-kinetic-parameter#answer_824414

编辑：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更早的评论隐藏 7更早的评论

Nadya Shalsabila salman 2021-11-7

在 MATLAB Online 中打开

Hi again! So I have plot again the kinetic model.. In this code below I added code to create two y-axis, on the right and left. but the result (plot and kinetic parameter) is still very far from the desired. Could you help me to solve my problem, please?

function terbaru
  
function C=kinetics(theta,t) 
c0=[0.6178; 28.90; 0.7586]; 
[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
yyaxis left
plot(t,c(:,[1 2]),'p')
hold on
hlp=plot(tv,Cfit(:,[1 2]));
hold off
ylabel('Concentration (g/L)')
ylim([min(c(:,1)) max(ylim)])
yyaxis right
plot(t,c(:,3),'p')
hold on
hlpp=plot(tv,Cfit(:,3));
hold off
ylabel('Concentration (\mug/L)')
ylim([min(c(:,3)) max(ylim)])
grid
legend([hlp([1 2]); hlpp],'C_1','C_2','C_3', 'Location','')
end 

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?

Alex Sha 2021-11-13