Monod nonlinear regression solved with Ode45

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

[Jan]

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

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.

.

[Jack95]

HI! Thanks for the reply, it was really helpful! I have only one question, if the parameters are 4, why do I have 6 values of B? what does the values of B in position 5 and 6 correspond to?

[Star Strider]

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]

Ah, I understand your reasoning. Another thing I don't understand is that B0 values make more scientific sense than B values, does this happen to you too?

[Star Strider]

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]

I don't know if you noticed, but the values of B change every time you run. maybe there is something wrong with the code, what do you think?

[Star Strider]

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]

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

[Star Strider]

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]

you are the best! Thank you very much!

[Star Strider]

My pleasure!

If my Answer helped you solve your problem, please Accept it!

.

[Star Strider]

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]

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]

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]

Hi Star Strider, thanks as always! I'm trying to use and study this code, but I'm getting an error related to the brace indexing:

Unable to perform assignment because brace indexing is not supported for variables of this type.

Error in testULTIMO (line 15)
    Cfit{k} = MonodKinetics1(Bk(:,k), t)

[Star Strider]

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]

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]

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]

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]

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]

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]

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.

[Jack95]

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]

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

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!

Answer 4

[Jack95]

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!