MATLAB Answers

How to obtain model parameters by fitting experimental data to the monod model. monod model.

2 ビュー (過去 30 日間)
O'Brien Ikart
O'Brien Ikart 2021 年 7 月 30 日
コメント済み: Alex Sha 2021 年 7 月 31 日
t= [0, 4, 8, 12, 16, 20, 24, 28, 32, 36, 40, 48]
x= [15, 18.3, 22, 23.4, 23.8, 24.1, 24.5, 24.8, 24.87, 24.9, 24.95, 25, 25]
S= [49.7, 28, 26.1, 18.1, 13.7, 10.1, 4.1, 1.2, 0.3, 0.08, 0.05, 0.05, 0.05 ]
P= [0, 5.2, 8.1, 9.3, 10.5, 12.3, 12.8, 13.2, 13.6, 13.8, 13.9, 13.95, 13.95]
I used initial guess values of umax, ks, y1 and y2 as 0.5, 55, 5 and 1.4 respectively. dx/dt= umax*s*x/(ks + s); ds/dt= -y1*umax*s*x/(ks + s) dp/dt= y2*umax*s*x/(ks +s). I used the code provided by star strider on a similar question but I keep getting errors, I also have to fit the data to other fermentation models. Also can simulink be used for the fitting?
  7 件のコメント
O'Brien Ikart
O'Brien Ikart 2021 年 7 月 30 日
Really new to Matlab and simulink

サインインしてコメントする。

採用された回答

Alex Sha
Alex Sha 2021 年 7 月 31 日
Refer to the results below:
Root of Mean Square Error (RMSE): 1.50437634361661
Sum of Squared Residual: 81.4733345963976
Correlation Coef. (R): 0.97672011839718
R-Square: 0.9539821896818
Parameter Best Estimate
-------------------- -------------
umax -0.0619295716977556
ks -88.5509341201568
y1 4.65073776898894
y2 1.33836136350295
  2 件のコメント
Alex Sha
Alex Sha 2021 年 7 月 31 日
If you want all parameters to be positive without upper bound limition, the result will be a bit strange as below:
Root of Mean Square Error (RMSE): 1.70454457017103
Sum of Squared Residual: 104.596998901184
Correlation Coef. (R): 0.976253765588915
R-Square: 0.953071414826535
Parameter Best Estimate
-------------------- -------------
y1 4.78475504434034
y2 1.35947688117111
umax 515071175352863
ks 5.10885251418351E17

サインインしてコメントする。

その他の回答 (1 件)

Star Strider
Star Strider 2021 年 7 月 30 日
It should be relatively straightforward to adapt your system to existing code.
One example is Parameter Estimation for a System of Differential Equations that might be easier to understand than the original Monod kinetics Answer, since I was building it as I was writing it.
That code in its entirety is:
function Igor_Moura
% 2016 12 03
% NOTES:
%
% 1. The ‘theta’ (parameter) argument has to be first in your
% ‘kinetics’ funciton,
% 2. You need to return ALL the values from ‘DifEq’ since you are fitting
% all the values
function C=kinetics(theta,t)
c0=[1;0;0;0];
[T,Cv]=ode45(@DifEq,t,c0);
%
function dC=DifEq(t,c)
dcdt=zeros(4,1);
dcdt(1)=-theta(1).*c(1)-theta(2).*c(1);
dcdt(2)= theta(1).*c(1)+theta(4).*c(3)-theta(3).*c(2)-theta(5).*c(2);
dcdt(3)= theta(2).*c(1)+theta(3).*c(2)-theta(4).*c(3)+theta(6).*c(4);
dcdt(4)= theta(5).*c(2)-theta(6).*c(4);
dC=dcdt;
end
C=Cv;
end
t=[0.1
0.2
0.4
0.6
0.8
1
1.5
2
3
4
5
6];
c=[0.902 0.06997 0.02463 0.00218
0.8072 0.1353 0.0482 0.008192
0.6757 0.2123 0.0864 0.0289
0.5569 0.2789 0.1063 0.06233
0.4297 0.3292 0.1476 0.09756
0.3774 0.3457 0.1485 0.1255
0.2149 0.3486 0.1821 0.2526
0.141 0.3254 0.194 0.3401
0.04921 0.2445 0.1742 0.5277
0.0178 0.1728 0.1732 0.6323
0.006431 0.1091 0.1137 0.7702
0.002595 0.08301 0.08224 0.835];
theta0=[1;1;1;1;1;1];
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c,zeros(size(theta0)));
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
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, 'C_1(t)', 'C_2(t)', 'C_3(t)', 'C_4(t)', 'Location','N')
end
I subsequently updated the figure code so the data and the curves fitting them are all the same colours. I include that update here.
.

Community Treasure Hunt

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

Start Hunting!

Translated by