Monod nonlinear regression solved with Ode45
2 visualizzazioni (ultimi 30 giorni)
Mostra commenti meno recenti
Hi everyone! I'm trying to estimate the parameters of the monod model solved with Ode45 through nlinfit, but I get the following errors:
Error in MONOD>ParameterJack (line 58)
mu = (Miumax*C(2))/(Ks+C(2));
Error in MONOD (line 10)
w=ParameterJack(ps,t);
Thanks in advance for the help!
The code is as follows
load dati.txt
y=dati;
t=y(:,1); X=y(:,2); S=y(:,3)
ps=[Miumax Ks Yxs];
w=ParameterJack(ps,t);
figure(1)
plot(t,w,'b.')
hold on
plot(t,yp,'r.')
pause
options=statset('Display','final','TolX',1e-12,'TolFun',1.e-12,'MaxIter',10000,'FunValCheck','off');
ip=0;
[pr,r,J] = nlinfit(t,yp,@ParameterJack,ps,options);
ci = nlparci(pr,r,J);
disp(ci);
w=ParameterJack(pr,t);
figure(2)
plot(t,w,'b.')
hold on
plot(t,yp,'.r')
function Substrate5
Rentang = [0:30];
C0 = [50 300]
Miumax=0.2;
Ks=10;
Yxs=0.5;
[t,C]=ode45(@(t,C)ParameterJack(t,C,ps),Rentang,C0);
figure
plot(t,C(:,1),'-g',t,C(:,2));
legend('acetogeni','idrogeno','metanogeni','acetato','metano', 'Location','best')
end
function dCdt=ParameterJack(t,C,Miumax,Ks,Yxs)
mu = (Miumax*C(2))/(Ks+C(2));
dCdt(1,:) = mu*C(1); %crescita acetogeni
dCdt(2,:) = -C(1)*(mu/Yxs) %consumo globale di idrogeno
end
1 Commento
Jan
il 8 Feb 2023
@jacopo ferretti: You have posted only the part of the error message, which explains where the error occurs, but not the part, wich mentions what the problem is. Please post the complete message.
You code would be much easier to read, if you avoid the pile of blank lines.
Risposte (4)
Star Strider
il 10 Feb 2023
Modificato: Star Strider
il 11 Feb 2023
It could possibly work with your data.
Note — As posted, ‘S’ has 21 elements, however ‘Rentang’ (that I assume is the associated time vector) has 31.
EDIT — (11 Feb 2023 at 16:49)
Using my Monod Kinetics code —
t=[0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20];
C=[50 52 55 58 61 64 67 70 71 71 71 71 71 71 71 71 71 71 71 71 7
370 327 281 234 184 132 78 23 0 0 0 0 0 0 0 0 0 0 0 0 0];
t = t(:);
C = C.';
B0 = [rand(4,1); C(1,:).'];
[B,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat] = lsqcurvefit(@MonodKinetics1,B0,t, C, zeros(size(B0)));
B % B(1:4): Kinetic Parameters, B(5:6): Initial Conditions
Cfit = MonodKinetics1(B, t);
figure
plot(t,C(:,1),'pb')
hold on
plot(t,C(:,2),'pr')
plot(t, Cfit(:,1), '-b')
plot(t, Cfit(:,2), '-r')
hold off
grid
function S = MonodKinetics1(B, t)
% MONODKINETICS1 codes the system of differential equations describing one
% version of Monod kinetics:
% dS/dt = (mmax * X * S)/(Y * (Ks + S));
% dX/dt = (mmax * X * S)/ (Ks + S) - (b * X);
% with:
% Variables: x(1) = S, x(2) = X
% Parameters: mmax = B(1), Y = B(2), Ks = B(3), b = B(4)
x0 = B(end-1:end);
[T,Sv] = ode45(@DifEq, t, x0);
function dS = DifEq(t, x)
xdot = zeros(2,1);
xdot(1) = (B(1) .* x(2) .* x(1)) ./ (B(2) .* (B(3) + x(1)));
xdot(2) = (B(1) .* x(2) .* x(1)) ./ (B(3) + x(1)) - (B(4) .* x(2));
dS = xdot;
end
S = Sv;
% S = Sv(:,1);
end
That is a reasonably decent approximation.
.
20 Commenti
Star Strider
il 4 Mar 2023
That is exactly what it should do.
You defined the time for the second integration as:
t1=7:14
so I incorporated that into my‘t’ matrix.
.
Jan
il 8 Feb 2023
A guess: You call the function ParameterJack with 2 inputs:
w=ParameterJack(ps,t);
But the function requires 5:
function dCdt=ParameterJack(t,C,Miumax,Ks,Yxs)
Then this lines fails:
mu = (Miumax*C(2))/(Ks+C(2));
because Miumax and Ks are not defined. The error message should reveal this already.
2 Commenti
Jan
il 10 Feb 2023
@jacopo ferretti: The error message is clear: "Vectors must be the same length." ParameterJack() uses the 1st 2 elements of C only. maybe you mean:
function dCdt=ParameterJack(t, C, Miumax, Ks, Yxs)
mu = (Miumax * C(2, :)) ./ (Ks + C(2, :));
dCdt(1,:) = mu * C(1, :);
dCdt(2,:) = -C(1, :) .* (mu / Yxs)
end
Vedere anche
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!