Fitting Monod Equation with ODE45 to data using lsqcurvefit function

106 visualizaciones (últimos 30 días)
Hello!
I am fitting Monod equation to a data containing substrate (s), biomass (x), and ethanol (p) concentration against time. The objective is to get the parameters: 1) umax, 2) ks, 3) Yxs, and 4)Yps that will best represent the data. The differential equations are:
Here is my initial code using assumed values of the four parameters:
umax = 0.5;
ks = 6.5;
Yxs = 0.2;
Yps = 1.2;
%a(1) = x
%a(2) = s
%a(3) = p
f = @(t,a) [umax*a(1)*a(2)/(ks + a(2)); -(1/Yxs)*umax*a(1)*a(2)/(ks + a(2)); (1/Yps)*umax*a(1)*a(2)/(ks + a(2))];
xt0 = [0.0904,9.0115,0.0151];
[tspan,a] = ode45(f,[0 25],xt0);
figure
plot(tspan,a(:,1),tspan,a(:,2),tspan,a(:,3))
Here is the code for trying to fit it into the actual data (script file):
function pos = paramfun1(x,tspan)
umax = x(1);
ks = x(2);
Yxs = x(3);
Yps = x(4);
xt0 = x(5:7);
f = @(t,a) [umax*a(1)*a(2)/(ks + a(2)); -(1/Yxs)*umax*a(1)*a(2)/(ks + a(2)); (1/Yps)*umax*a(1)*a(2)/(ks + a(2))];
[~,pos] = ode45(f,tspan,xt0);
Here is my call function (in the command window):
xt0 = zeros(1,7);
xt0(1) = umax;
xt0(2) = ks;
xt0(3) = Yxs;
xt0(4) = Yps;
data =[0 3 5 8 9.5 11.5 14 16 18 20 25 27, 0.0904 0.1503 0.2407 0.3864 0.5201 0.6667 0.8159 0.9979 1.0673 1.1224 1.1512 1.2093; 0 3 5 8 9.5 11.5 14 16 18 20 25 27, 9.0115 8.8088 7.9229 7.2668 5.3347 4.911 3.5354 1.4041 0 0 0 0; 0 3 5 8 9.5 11.5 14 16 18 20 25 27, 0.0151 0.0328 0.0621 0.1259 0.2949 0.3715 0.4199 0.522 0.5345 0.6081 0.07662 0.7869];
%time =[0 3 5 8 9.5 11.5 14 16 18 20 25 27];
[pbest,exitflag,output] = lsqcurvefit(@paramfun,xt0,tspan,data);
fprintf('New parameters: %f, %f, %f, %f',pbest(1:4));
The error is function value not equal to YDATA. Btw, this code was based from an example in MATLAB. (https://www.mathworks.com/help/optim/ug/fit-differential-equation-ode.html)
My data is:
time = [0 3 5 8 9.5 11.5 14 16 18 20 25 27]
x = 0.0904 0.1503 0.2407 0.3864 0.5201 0.6667 0.8159 0.9979 1.0673 1.1224 1.1512 1.2093
s = 9.0115 8.8088 7.9229 7.2668 5.3347 4.911 3.5354 1.4041 0 0 0 0
p = 0.0151 0.0328 0.0621 0.1259 0.2949 0.3715 0.4199 0.522 0.5345 0.6081 0.07662 0.7869
Please help! I do not know how to input my data into the lsqcurvefit function.
Thanks in advance!

Respuesta aceptada

David Wilson
David Wilson el 22 de Abr. de 2019
Editada: David Wilson el 22 de Abr. de 2019
Here's my solution.
First set up the measured data, plot it, and look to see that it is all OK.
t = [0 3 5 8 9.5 11.5 14 16 18 20 25 27]'
x = [0.0904 0.1503 0.2407 0.3864 0.5201 0.6667 0.8159 0.9979 1.0673 1.1224 1.1512 1.2093]';
s = [9.0115 8.8088 7.9229 7.2668 5.3347 4.911 3.5354 1.4041 0 0 0 0]';
p = [0.0151 0.0328 0.0621 0.1259 0.2949 0.3715 0.4199 0.522 0.5345 0.6081 0.07662 0.7869]';
plot(t,[x,s,p],'o-')
No problems there. Note that the measured data is in columns.
Now, we need a function to export the predicted measurements as a function of the unknown parameters and initial conditions. My function is:
function ypred = paramfun1(p,t)
umax = p(1); ks = p(2); Yxs = p(3); Yps = p(4); % use disperse here ...
xt0 = p(5:7); % initial conditions
f = @(t,a) [umax*a(1)*a(2)/(ks + a(2)); ...
-(1/Yxs)*umax*a(1)*a(2)/(ks + a(2)); ...
(1/Yps)*umax*a(1)*a(2)/(ks + a(2))];
[~,ypred] = ode45(f,t,xt0);
end
May I suggest that you think about naming your functions sensible names. Also it is more elegant to use disperse (from the user's group) here.
Now we are ready to test it, and make sure our output is the same shape as the actual data.
xt0 = [0.1;9;0.1]; % read off from data plot [x0 ,S0, p0]';
p0 = [umax; ks; Yxs; Yps; xt0]
Ypred = paramfun1(p0,t); % measurement predictions
Now we are ready for the fit. I'm using your initial conditions, and we can get reasonable start values from the data.
pfit = lsqcurvefit(@paramfun1,p0,t,[x,s,p])
Now plot the fitted curve and the actual data. They should line up.
umax = pfit(1); ks = pfit(2);Yxs = pfit(3); Yps = pfit(4);
f = @(t,a) [umax*a(1)*a(2)/(ks + a(2)); -(1/Yxs)*umax*a(1)*a(2)/(ks + a(2)); (1/Yps)*umax*a(1)*a(2)/(ks + a(2))];
xt0 = pfit(5:7);
[tspan,a] = ode45(f,[0 25],xt0);
plot(tspan,a(:,1),tspan,a(:,2),tspan,a(:,3), '-', ...
t,x,'o',t,s,'+',t,p,'s')
Plot shows the regressed fit is not too bad:
  7 comentarios
Vinoj Liyanaarachchi
Vinoj Liyanaarachchi el 19 de Oct. de 2019
Thank you StarSign1997 !!
Are there any techniques to improve the parameter estimation? Above shapes of the curves are not acceptable in my case (microalgae cultivation).
Thanks in advance!
Antara Mazumder
Antara Mazumder el 6 de Jun. de 2020
Hi this was really helpful, but I dont know I worte exactlr same code, but it showing me an error for using options saying " Empty keys are not allowed in this container.", dint understand this.

Iniciar sesión para comentar.

Más respuestas (2)

Alex Sha
Alex Sha el 21 de Oct. de 2019
The follow results obtained from 1stOpt may be used for comparison:
Code:
Variable t,x,s,p;
ODEFunction x'=umax*s*x/(ks + s);
s'=-(1/Yxs)*umax*s*x/(ks + s);
p'=(1/Yps)*umax*s*x/(ks + s);
Data;
t = [0 2 4 6 8 10 12 14 16];
x = [0.06 0.11 0.46 0.78 1.42 2.36 2.49 2.49 2.54];
s = [9.54 9.33 9.52 9.06 8.05 7.27 6.03 5.80 5.51];
p = [0.00 0.00 0.00 0.01 0.18 0.35 0.49 0.53 0.55];
Results:
Root of Mean Square Error (RMSE): 0.250520910244571
Sum of Squared Residual: 1.50625743527444
Correlation Coef. (R): 0.981509121545543
R-Square: 0.963360155677103
Adjusted R-Square: 0.914727231437843
Determination Coef. (DC): 0.942399473562309
F-Statistic: 14.7361249635671
Parameter Best Estimate
-------------------- -------------
umax 0.117013263727751
ks -5.83263810513223
yxs 0.699511367324087
yps 5.01252941213249
c215.jpg
c216.jpg
c217.jpg
  2 comentarios
Vinoj Liyanaarachchi
Vinoj Liyanaarachchi el 21 de Oct. de 2019
Thank you Alex Sha. Actually I'm new to matlab.
Can you share the code..
Thanks in advanse!
Rifat Hasan
Rifat Hasan el 2 de Nov. de 2022
Hi ! I have tried to run the code by David welson. I am getting an error.
Solver stopped prematurely.
lsqcurvefit stopped because it exceeded the function evaluation limit,
options.MaxFunctionEvaluations = 7.000000e+02.
How can I solve this issue?

Iniciar sesión para comentar.


Alex Sha
Alex Sha el 21 de Oct. de 2019
Hi, Vinoj Liyanaarachchi, the code have been provided above alreadly, but note, that code is not for Matlab, but for another package called 1stOpt. In Matlab, no matter lsqcurvefit or fminsearch functions, are all to adopt local optimization algorithms, it is why those functions depended haveily on initial start-values. Global optimization toolbox in Matlab (like GA) is not good enough to ensure the global solution.

Categorías

Más información sobre Interpolation en Help Center y File Exchange.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by