Asked by StarSign1997
on 22 Apr 2019

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!

Answer by David Wilson
on 22 Apr 2019

Edited by David Wilson
on 22 Apr 2019

Accepted Answer

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:

StarSign1997
on 22 Apr 2019

Thank you so much! Actually, I'll be doing a model discrimination among three models. One of them is Monod. Okay I understand now how the flow of passing the argument is supposed to be.

Thank you again!

StarSign1997
on 23 Apr 2019

StarSign1997
on 23 Apr 2019

I got it. I used the finite difference step size found in the example. The final parameter values are not anymore so sensitive to the initial guess and almost has the same value regardless of my initial guess.

Thanks!

Sign in to comment.

Opportunities for recent engineering grads.

Apply Today
## 1 Comment

## Star Strider (view profile)

Direct link to this comment:https://de.mathworks.com/matlabcentral/answers/457752-fitting-monod-equation-with-ode45-to-data-using-lsqcurvefit-function#comment_696765

Sign in to comment.