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

5 views (last 30 days)
O'Brien Ikart on 30 Jul 2021
Commented: Alex Sha on 31 Jul 2021
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?
O'Brien Ikart on 30 Jul 2021
Really new to Matlab and simulink

Alex Sha on 31 Jul 2021
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
Alex Sha on 31 Jul 2021
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

Star Strider on 30 Jul 2021
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.
.
Star Strider on 30 Jul 2021
My pleasure!
.