Monod nonlinear regression solved with Ode45

Question

Jack95 2023-2-8

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1908940-monod-nonlinear-regression-solved-with-ode45

评论： Jack95 2023-3-14

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 个评论
显示 -1更早的评论隐藏 -1更早的评论

Jan 2023-2-8

@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.

请先登录，再进行评论。

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

Answer 1

Star Strider 2023-2-10

1
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1908940-monod-nonlinear-regression-solved-with-ode45#answer_1168540

编辑：Star Strider 2023-2-11

在 MATLAB Online 中打开

See: Monod kinetics and curve fitting .

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)));

Local minimum possible. lsqcurvefit stopped because the final change in the sum of squares relative to its initial value is less than the value of the function tolerance.

B % B(1:4): Kinetic Parameters, B(5:6): Initial Conditions

B = 6×1

0.1803 17.0729 0.0000 0.4445 50.4131 415.1698

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 个评论
显示 18更早的评论隐藏 18更早的评论

Star Strider 2023-2-13

在 MATLAB Online 中打开

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)

Start Time: 2023-02-13 17:03:14.8045

[theta,fval,exitflag,output,population,scores] = ga(ftns, Parms, [],[],[],[],zeros(Parms,1),Inf(Parms,1),[],[],optsAns);

Optimization terminated: average change in the fitness value less than options.FunctionTolerance.

t1 = clock;

fprintf('\nStop Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t1)

Stop Time: 2023-02-13 17:03:31.3783

GA_Time = etime(t1,t0)

GA_Time = 16.5738

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)

Elapsed Time: 1.657384000000000E+01 00:00:16.573

fprintf('Fitness value at convergence = %.4f\nGenerations \t\t\t\t = %d\n\n',fval,output.generations)

Fitness value at convergence = 124.6317 Generations = 226

fprintf(1,'\tRate Constants:\n')

Rate Constants:

for k1 = 1:length(theta)

fprintf(1, '\t\tTheta(%2d) = %8.5f\n', k1, theta(k1))

end

Theta( 1) = 0.16700 Theta( 2) = 9.16559 Theta( 3) = 0.00005 Theta( 4) = 0.43100 Theta( 5) = 48.00156 Theta( 6) = 415.04239

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 2023-2-24

在 MATLAB Online 中打开

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)));

Solver stopped prematurely. lsqcurvefit stopped because it exceeded the function evaluation limit, options.MaxFunctionEvaluations = 6.000000e+02.

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'})

ParametersCI = 6×3 table

Lower 95% CI Parmaeter Estimate Upper 95% CI _____________________ ____________________ ____________________ -2.31754458605142e+03 1.47144566342795e-01 2.31783887518410e+03 -1.04014031416879e+02 8.00239363558959e-02 1.04174079289590e+02 -1.84284312807881e+08 1.15972933777972e+04 1.84307507394637e+08 -7.29027652795074e-01 2.49700343777263e-01 1.22842834034960e+00 9.94662444226436e+00 5.25197386125868e+01 9.50928527829093e+01 3.40198632474659e+02 3.85534641673575e+02 4.30870650872491e+02

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 2023-3-2

在 MATLAB Online 中打开

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 2023-3-2

编辑：Star Strider 2023-3-3

在 MATLAB Online 中打开

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

Local minimum possible. lsqcurvefit stopped because the final change in the sum of squares relative to its initial value is less than the value of the function tolerance.

Cfit = 1×1 cell array

{8×2 double}

Solver stopped prematurely. lsqcurvefit stopped because it exceeded the function evaluation limit, options.MaxFunctionEvaluations = 1.200000e+03.

Cfit = 1×2 cell array

{8×2 double} {8×2 double}

% 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 2023-3-4

在 MATLAB Online 中打开

Yes, it works! now i'm trying to split the two C matrices in two times, in order to represent my data. i'm trying to simply plot the respective C's at defined times, but something is wrong, maybe it's related to ode45? this is a bit advanced for me, but if i solve this problem i can finally start my study,

Thanks again!

% 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;
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, 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{1}(:,1),'pb')
    plot(t1, Cc{2}(:,1),'pb')
     plot(t, Cc{1}(:,2),'pr')
    plot(t1, Cc{2}(:,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

Star Strider 2023-3-4

在 MATLAB Online 中打开

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

Solver stopped prematurely. lsqcurvefit stopped because it exceeded the function evaluation limit, options.MaxFunctionEvaluations = 1.200000e+03. Solver stopped prematurely. lsqcurvefit stopped because it exceeded the function evaluation limit, options.MaxFunctionEvaluations = 1.200000e+03.

% 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 2023-3-4

在 MATLAB Online 中打开

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 2023-3-4

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 2023-3-14

在 MATLAB Online 中打开

HI! I'm back to working on this code, I'm trying to add a third dataset (and then I'll have to add yet another). I generally work by loading data from txt files, but in this case I don't quite understand how to do it, because there is an error that I can't find. can you explain me quickly? thanks as always!

% 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}=[113 128 123 96 136 117 120 151 
   387 333 320 302 254 187 124 87
   0 26 75 81 87 91 167 227]; % I want to add this serie
Cc{2}=[170 182 195 208 211 224 257 290
   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)) * x(2);
          xdot(3) = 1/B(2) .*x(2)
          dS = xdot;
      end
S = Sv;
% S = Sv(:,1);
end

请先登录，再进行评论。

Answer 2

Jan 2023-2-8

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1908940-monod-nonlinear-regression-solved-with-ode45#answer_1166730

在 MATLAB Online 中打开

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 个评论
显示无隐藏无

Jack95 2023-2-10

在 MATLAB Online 中打开

HI! Thanks for the reply! I removed the spaces from the code and put the data directly in the code at least you can try it directly if you want. I followed your indication and now the problem seems to be another:

>> MONOD

dCdt =

9.7368

-19.4737

Error using plot

Vectors must be the same length.

Error in MONOD (line 12)

plot(t,w,'b.')

The code is:

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];
Miumax=0.2;
Ks=10; 
Yxs=0.5; 
ps=[Miumax Ks Yxs];
w=ParameterJack(t,C,Miumax,Ks,Yxs);
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,Miumax,Ks,Yxs),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);                            
    dCdt(2,:) = -C(1)*(mu/Yxs)      
end

Jan 2023-2-10

在 MATLAB Online 中打开

@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

请先登录，再进行评论。

Answer 3

Jack95 2023-2-24

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1908940-monod-nonlinear-regression-solved-with-ode45#answer_1179145

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!

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Answer 4

Jack95 2023-3-2

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1908940-monod-nonlinear-regression-solved-with-ode45#answer_1184010

在 MATLAB Online 中打开

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!

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

请先登录，再进行评论。

Monod nonlinear regression solved with Ode45

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

回答（4 个）

20 个评论
显示 18更早的评论隐藏 18更早的评论

2 个评论
显示无隐藏无

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Monod nonlinear regression solved with Ode45

1 个评论 显示 -1更早的评论隐藏 -1更早的评论

回答（4 个）

20 个评论 显示 18更早的评论隐藏 18更早的评论

2 个评论 显示 无隐藏 无

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

0 个评论 显示 -2更早的评论隐藏 -2更早的评论

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

1 个评论
显示 -1更早的评论隐藏 -1更早的评论

20 个评论
显示 18更早的评论隐藏 18更早的评论

2 个评论
显示无隐藏无

0 个评论
显示 -2更早的评论隐藏 -2更早的评论

0 个评论
显示 -2更早的评论隐藏 -2更早的评论