Monod nonlinear regression solved with Ode45
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
Jack95
il 13 Feb 2023
Star Strider
il 13 Feb 2023
My pleasure!
I intended to explain that in the comment with respect to ‘B’ —
B % B(1:4): Kinetic Parameters, B(5:6): Initial Conditions
When I fit differential equations to data, I always include the initial conditions as parameters to be estimated, and they are always the last values in the parameter vector. All the measurements are (in my opinion) subject to the same measurement errors, so fitting the initial conditions not only makes the fitting possible, it also improves the fit by estimating the actual initial conditions (beginning variable values).
In ‘MonodKinetics1’ I specify that as:
x0 = B(end-1:end);
where ‘B’ is the parameter vector and ‘x0’ are the initial conditions.
.
Jack95
il 13 Feb 2023
Star Strider
il 13 Feb 2023
I do not understand.
The ‘B0’ vector are just the initial parameter estimates. They are modified in the fitting process to fit the data. (The only constraint I use is a lower bound of zero for all parameters, since none of them should be negative in this instance.)
Jack95
il 13 Feb 2023
I noticed, however that is typical for nonlinear parameter estimation routines. The parameters will differ, depending on the initial parameter estimates in ‘B0’. One option is to use the Global Optimization Toolbox ga function, or another global optimiser, to search the parameter space for the best set of parameters (lowest residual norm, or norm of the residuals), however even that may not produce unique values for the parameters, since it is always possible for slight variations in the parameters to produce the same value for the lowest residual norm. (The ga function allows its parameter estimates to be ‘fine tuned’ using a gradient descent optimiser such as fmincon. See 'HybridFcn' in the options documentation section.) The lsqcurvefit function uses a gradient-descent approach, and even if there is a global minimum, the parameter values may not always be the same, since it stops iterating when the tolerances are acceptable. Changing the appropriate tolerances using optimoptions (or the appropriate options function for global optimisation functions) can produce more precise parameter estimates, however most of the time this is not necessary.
A genetic algorithm approach would be —
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 71
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.';
ftns = @(theta) norm(C-MonodKinetics1(theta,t));
PopSz = 250;
Parms = 6;
optsAns = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms).*[ones(1,4)*1E-3 1E-2 1E-1], 'MaxGenerations',5E3, 'FunctionTolerance',1E-10); % Options Structure For 'Answers' Problems
t0 = clock;
fprintf('\nStart Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t0)
[theta,fval,exitflag,output,population,scores] = ga(ftns, Parms, [],[],[],[],zeros(Parms,1),Inf(Parms,1),[],[],optsAns);
t1 = clock;
fprintf('\nStop Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t1)
GA_Time = etime(t1,t0)
DT_GA_Time = datetime([0 0 0 0 0 GA_Time], 'Format','HH:mm:ss.SSS');
fprintf('\nElapsed Time: %23.15E\t\t%s\n\n', GA_Time, DT_GA_Time)
fprintf('Fitness value at convergence = %.4f\nGenerations \t\t\t\t = %d\n\n',fval,output.generations)
fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(theta)
fprintf(1, '\t\tTheta(%2d) = %8.5f\n', k1, theta(k1))
end
tv = linspace(min(t), max(t));
Cfit = MonodKinetics1(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, compose('C_%d',1:size(C,2)), 'Location','N')
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
Experiment with it with different tolerances and other options variations.
Star Strider
il 24 Feb 2023
@jacopo ferretti Commented:
I think I can stay on this topic. How do you set the confidence interval (e.g. 95%) of the parameters? i'm trying with nlparci but it doesn't work with lsqcurvefit. Thanks again!
Yes, it does!
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)));
ci = nlparci(B,Rsd,'Jacobian',Jmat);
% B % B(1:4): Kinetic Parameters, B(5:6): Initial Conditions
fmt = format;
format longE
ParametersCI = table(ci(:,1), B, ci(:,2), 'VariableNames',{'Lower 95% CI','Parmaeter Estimate','Upper 95% CI'})
format short
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
Only parameters 5 and 6 are significantly different from zero at the 95% confidence limit.
.
Jack95
il 25 Feb 2023
Star Strider
il 25 Feb 2023
My pleasure!
If my Answer helped you solve your problem, please Accept it!
.
Star Strider
il 2 Mar 2023
Change these calls:
zeros(Tspan1,(theta0))
zeros(Tspan2,(theta0))
to:
zeros(numel(Tspan1),numel(theta0))
zeros(numel(Tspan2),numel(theta0))
If those are the only errors, the code should work.
Jack95
il 2 Mar 2023
it does not work. So maybe the problem is another but I don't know which one..
How could this be done? taking as reference this code for example, where in the substrate I changed to day 8 with a new addition:
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 370 327 281 250 230 200 180 170 140 100 40 30 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)));
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)
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
Star Strider
il 2 Mar 2023
Modificato: Star Strider
il 3 Mar 2023
All I did was to see if I could get your code to run, and found the previous error. My originally posted code works with the provided data.
You probably need to run this separately for each row of ‘C’ (using a for loop and saving the results each time) and then combine the results either in one plot (use the hold function) or different plots for each set of data (row of ‘C’).
EDIT — (3 Mar 2023 at 23:13)
The ‘C’ matrix needs to be re-written.
Try this —
% 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 71
% 370 327 281 234 184 132 78 23 370 327 281 250 230 200 180 170 140 100 40 30 0];
t = 0:7;
Cc{1}=[50 52 55 58 61 64 67 70
370 327 281 234 184 132 78 23];
Cc{2}=[50 52 55 58 61 64 67 70
370 327 281 250 230 200 180 170];
t = t(:);
% C = Cc.';
B0 = [rand(4,1); Cc{1}(1,:).'];
for k = 1:2
[B,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat] = lsqcurvefit(@MonodKinetics1,B0,t, Cc{k}.', zeros(size(B0)));
Bk(:,k) = B;
Cfit{k} = MonodKinetics1(Bk(:,k), t)
end
% Cfit = MonodKinetics1(B, t);
figure
% plot(t,Ck{k}(:,1),'pb')
hold on
% plot(t,C(:,2),'pr')
for k = 1:2
Cc{k} = Cc{k}.';
% Cfit{k}
plot(t, Cc{k}(:,1),'pb')
plot(t, Cc{k}(:,2),'pr')
plot(t, Cfit{k}(:,1), '-b')
plot(t, Cfit{k}(:,2), '-r')
end
hold off
grid
function S = MonodKinetics1(B, t)
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
.
Jack95
il 4 Mar 2023
Star Strider
il 4 Mar 2023
You need to copy and run my latest code exactly as I wrote it. If ‘Cfit’ was defined as something other than a cell array earlier, the code will throw that error. I only define it as a cell array here.
Jack95
il 4 Mar 2023
Create a second time vector and reference it appropriately —
% 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 71
% 370 327 281 234 184 132 78 23 370 327 281 250 230 200 180 170 140 100 40 30 0];
t = [0:7; 7:14]; % Concatenate Time Vectors
% t1=7:14
Cc{1}=[50 62 75 88 91 104 127 170
370 327 281 234 184 132 78 23];
Cc{2}=[170 182 195 208 211 224 257 190
420 307 261 240 220 180 140 110];
t = t.';
% C = Cc.';
B0 = [rand(4,1); Cc{1}(1,:).'];
for k = 1:2
[B,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat] = lsqcurvefit(@MonodKinetics1,B0,t(:,k), Cc{k}.', zeros(size(B0)));
Bk(:,k) = B;
Cfit{k} = MonodKinetics1(Bk(:,k), t(:,k));
end
% Cfit = MonodKinetics1(B, t);
figure
% plot(t,Ck{k}(:,1),'pb')
hold on
% plot(t,C(:,2),'pr')
for k = 1:2
Cc{k} = Cc{k}.';
tv = t(:,k);
plot(tv, Cc{k}(:,1),'pb')
plot(tv, Cc{k}(:,2),'pr')
plot(tv, Cfit{k}(:,1), '-b')
plot(tv, Cfit{k}(:,2), '-r')
end
hold off
grid
function S = MonodKinetics1(B, t)
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
Also, plot the ‘Cfit’ vectors.
.
Jack95
il 4 Mar 2023
wonderful, you are the best! However, one thing seems strange to me: at time zero I should have the same value of the corresponding matrix, instead it generates this:
t = [0:7; 7:14]; % Concatenate Time Vectors
% t1=7:14
Cc{1}=[43 52 72 118 137 150 165 190
302 290 273 244 214 177 150 124];
Cc{2}=[190 214 232 290 317 324 340 330
420 307 261 240 220 180 140 110];
t = t.';
% C = Cc.';
B0 = [rand(4,1); Cc{1}(1,:).'];
for k = 1:2
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat] = lsqcurvefit(@MonodKinetics1,B0,t(:,k), Cc{k}.', zeros(size(B0)));
Bk(:,k) = theta;
Cfit{k} = MonodKinetics1(Bk(:,k), t(:,k));
end
figure
hold on
for k = 1:2
Cc{k} = Cc{k}.';
tv = t(:,k);
plot(tv, Cc{k}(:,1),'pb')
plot(tv, Cc{k}(:,2),'pr')
plot(tv, Cfit{k}(:,1), '-b')
plot(tv, Cfit{k}(:,2), '-r')
end
hold off
grid
function S = MonodKinetics1(theta, t)
x0 = theta(end-1:end);
u=0.1
[T,Sv] = ode45(@DifEq, t, x0);
function dS = DifEq(t, x)
xdot = zeros(2,1);
xdot(1) = (u.* x(1) .* x(2)) ./ (theta(2) + x(2)) - (theta(3) .* x(1));
xdot(2) = -(u .* x(1) .* x(2)) ./ (theta(2) .* (theta(2) + x(2)));
dS = xdot;
end
S = Sv;
% S = Sv(:,1);
end
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.
.
Jack95
il 14 Mar 2023
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
Jack95
il 10 Feb 2023
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
Jack95
il 24 Feb 2023
0 voti
Jack95
il 2 Mar 2023
Hi everyone! based on this code:
Tspan1 = 0:20;
Tspan2 = 20:70;
C01 = [50 300 0 10 0];
Miumax=0.07;
Ks=10;
Yxs=0.1;
alfa=0.9;
betha=0;
y=0.2;
km=100;
ka=3;
a2=0.8;
a1=0.559;
u=0.1;
[t1,C1] = ode45(@(t,C)ParameterJack(t,C,betha),Tspan1,C01);
C02 = C1(end,:);
C02(1,2) = 300;
[t2,C2] = ode45(@(t,C)ParameterJack(t,C,betha),Tspan2,C02);
t = [t1; t2];
C = [C1; C2];
figure
hold on
yyaxis left
plot(t,C(:,1),'-g',t,C(:,2),'-b',t,C(:,4),'-m');
ylabel('acetogens (mgCOD/L), hidrogen (mgCOD/L),metanogens')
ylim([min(C(:,1)) max(ylim)])
ylim([min(C(:,2)) max(ylim)])
yyaxis right
plot(t,C(:,3),'r', t,C(:,5),'-o');
ylabel('acetate (mgCOD/L), methane (mgCOD/L)')
ylim([0 max(ylim)])
hold off
xlabel('tempo (giorni)')
grid
legend('acetogens','Hidrogen','metanogens','acetate','methane', 'Location','best')
function dCdt=ParameterJack(t,C,betha)
Miumax=0.07;
alfa=0.9;
u=0.1;
Ks = 80;
y=0.2;
km=100;
Yxs = 0.175;
mu = (Miumax*C(2))/(Ks+C(2));
if t<24
term2 = 0,
term3 = 0;
term4 = 0;
term5 = 0;
else
ka = 60;
a2=0.121;
km=10;
y=0.2
a3= (C(4)/C(5))*a2
mum = (u)*(C(2))/(km+C(2));
mo=(u)*(C(3))/(ka+C(3));
u= 0.1
term2 = (C(4)*(mum/(1-y)));
term3 = (C(4)*(mo/(1-a2)));
term4 = mum*C(4)+(u)*((C(3)/((C(3))+ka))*C(4));
term5 = term3+(C(2)*y);
end
dCdt(1,:) = mu*C(1);
dCdt(2,:) = -C(1)*(mu/Yxs)-term2;
dCdt(3,:) = (C(1)*((alfa*mu)+betha))-term3;
dCdt(4,:) = term4;
dCdt(5,:) = term5;
end
where I have two tspans which correspond to two moments in which I administered substrate, and the graph is:
I'm now trying to figure out how to bring the same concept back into the nonlinear regression code (considering for now only the terms hydrogen, biomass, acetate).
clear all
load data2.txt
y=data2;
t=y(:,1);
x=y(:,2);
s=y(:,3);
p=y(:,4);
c = [x s p];
Tspan1 = 0:20;
Tspan2 = 20:70;
options = optimset('MaxFunEvals',100000,'MaxIter',100000);
theta0=rand(3,1);
options = optimoptions('lsqcurvefit','Algorithm','levenberg-marquardt');
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c,zeros(Tspan1,(theta0)));
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c,zeros(Tspan2,(theta0)));
ci = nlparci(theta,Rsd,'Jacobian',Jmat);
fmt = format;
format longE
ParametersCI = table(ci(:,1), theta, ci(:,2), 'VariableNames',{'Lower 95% CI','Parmaeter Estimate','Upper 95% CI'})
C = kinetics(theta,t);
hold on
plot(t,c(:,1),'pb')
plot(t,c(:,2),'pr')
plot(t,c(:,3),'pg')
plot(t, C(:,1), '-b')
plot(t, C(:,2), '-r')
plot(t, C(:,3), '-g')
legend('biomassa','idrogeno','acetato')
hold off
function C=kinetics(theta,t)
%c0=[0.1;9;0.1];
c0 = [20;375;13];
u=0.1
[T1,C1]=ode45(@DifEq,t,c0); %FISSA Yxs e umax E TROVA Ks!
C02 = C1(end,:);
C02(1,2) = 300;
%
function dC=DifEq(t,c)
dcdt=zeros(3,1);
dcdt(1)=u*(c(2)/(theta(1)+c(2)))*c(1);
dcdt(2)= -(1/theta(2))*u*c(1)*c(2)/(theta(1)+c(2));
dcdt(3)= (1/theta(3))*u*c(1)*c(2)/(theta(1)+c(2));
dC=dcdt;
end
C=C;
ERROR:
>> testrefil
Error using zeros
Size inputs must be scalar.
Error in testrefil (line 15)
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c,zeros(Tspan1,(theta0)));
I'm trying to work this way, but I think there is some programming error! it's a bit advanced for me and I can't find the error, if you can help me I would be grateful!
Categorie
Scopri di più su Descriptive Statistics in Centro assistenza e File Exchange
il 8 Feb 2023
il 14 Mar 2023
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
Translated by 
