Finding Unknown Parameters from Monod Equations
Show older comments
Hi, new to Matlab and still learning.
I am using a modified code provided by Star Strider.Thanks Star Strider! I am making this as a new post so I can accept this one as well =).
The code is shown below. I got the data from another Matlab post. The equation is as shown in the screenshot: I converted all constants into theta. When I run the code below, I get huge constants that doesn't makes sense. Anything I'm doing wrong? Is it my equations? C0 or theta0? Thank you!
Function Igor_Mour
function C=kinetics(theta,t)
c0=[0.1;9;0.1];
[T,Cv]=ode45(@DifEq,t,c0);
%
function dC=DifEq(t,c)
dcdt=zeros(3,1);
dcdt(1)=-theta(1).*c(1).*c(2)/(theta(2)+c(1));
dcdt(2)= -(1/theta(3)).*theta(1).*c(1).*c(2)/(theta(2)+c(1));
dcdt(3)= (1/theta(4)).*theta(1).*c(1).*c(2)/(theta(2)+c(1));
dC=dcdt;
end
C=Cv;
end
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]';
%Note: x,s and p are combined (in separate columns) to form c.
theta0=[1;1;1;1;1;1];
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c);
fprintf(1, '\tRate Constants:\n')
for k1 = 1:length(theta)
fprintf(1, '9\t\tTheta(%d) = %8.5f\n', k1, theta(k1))
end
tv = linspace(min(t), max(t));
Cfit = kinetics(theta, tv);
figure(1)
plot(t, c, 'p')
hold on
hlp = plot(tv, Cfit);
hold off
grid
xlabel(Time)
ylabel(Concentration)
legend(hlp, 'C_1(t)','C_2(t)','C_3(t)','C_4(t)', 'Location','N')
end
1 Comment
Alex Sha
on 16 Oct 2022
Hi, 0.07662, one of data for p, should be 0.7662? if so, see the result below:
Sum Squared Error (SSE): 1.43925599360388
Root of Mean Square Error (RMSE): 0.208839215637285
Correlation Coef. (R): 0.982548927303979
R-Square: 0.9654023945462
Parameter Best Estimate
--------- -------------
theta1 0.0214068471653595
theta2 -1.28632596117882
theta3 -0.113732223600425
theta4 -1.55421915737194
Accepted Answer
The estimated value for ‘theta(1)’ definitely does not appear to be as expected for a rate constant, however if it is 
than it could be appropriate. I am not certain what the other parameters are by comparing the code to the symbolic expressions It would help to have them clarified.
than it could be appropriate. I am not certain what the other parameters are by comparing the code to the symbolic expressions It would help to have them clarified. According to the posted symbolic differntial equations, the first differential equation should not be negated. Changing that gives a decent fit to the data. Beyond that, I suggest checking to be certain that the equations are coded correctly. Here, I have ‘x’ as ‘c(1)’, ‘s’ as ‘c(2)’, and ‘p’ as ‘c(3)’. Change those if necessary to be the correct columns.
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]';
%Note: x,s and p are combined (in separate columns) to form c.
c = [x s p];
theta0=rand(7,1);
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c,zeros(1,numel(theta0)));
fprintf(1,'\tRate Constants:\n')
for k1 = 1:numel(theta)
fprintf(1, '\t\tTheta(%2d) = %8.5f\n', k1, theta(k1))
end
tv = linspace(min(t), max(t));
Cfit = kinetics(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','best')
function C=kinetics(theta,t)
c0=theta(5:7);
[T,Cv]=ode45(@DifEq,t,c0);
%
function dC=DifEq(t,c)
dcdt=zeros(3,1);
dcdt(1)= theta(1).*c(1).*c(2)./(theta(2)+c(1));
dcdt(2)= -(1./theta(3)).*theta(1).*c(1).*c(2)./(theta(2)+c(1));
dcdt(3)= (1./theta(4)).*theta(1).*c(1).*c(2)./(theta(2)+c(1));
dC=dcdt;
end
C=Cv;
end
Please define the parameters in terms of the symbolic equations. With that information, it will be easier to determine if there are any inconsistencies.
.
12 Comments
Learning
on 14 Oct 2022
Thanks for your responses! not sure of the symbolic equations you are referring to but they were written as shown below.
OR:
dX/dt = umax*a(1)*a(2)/(ks + a(2));
dS/dt = -(1/Yxs)*umax*a(1)*a(2)/(ks + a(2));
dP/dt = (1/Yps)*umax*a(1)*a(2)/(ks + a(2))];
I made umax as K(1), ks as K(2), Yxs as K(3), and Yps as K(4), s as c(1), and x as c(2).
I found a typo in my diff equations in my question. The equations should be
dcdt=zeros(3,1);
dcdt(1)=-theta(1).*c(1).*c(2)/(theta(2)+c(1));
dcdt(2)= -(1/(theta(3))).*theta(1).*c(1).*c(2)/(theta(2)+c(1));
dcdt(3)= (1/(theta(4))).*theta(1).*c(1).*c(2)/(theta(2)+c(1));
dC=dcdt;
Few questions
- From the equations above, I have only s and x and hence I'm thinking only 2 concentrations, c(1) and c(2)? where does the c(3) arises from? from dP/dt?
- How do you choose the length of theta0, you have theta0 = [1;1;1;1;1;1;1]. How do you know you have to have 7 1s and not 4 1s?
- When you were plotting the data after estimating theta's you made c0=theta(5:7); is it wrong if I made c0 be the same as the original initial concentrations?,
Thanks!
The first ‘dcdt(1)’ should not be negated, according to the symbolic equations. I now understand that ‘theta(1)’ is 
however ‘c(1)=s’ and ‘c(2)=x’ changes the code so that it now does not work.
however ‘c(1)=s’ and ‘c(2)=x’ changes the code so that it now does not work. Copying and pasting your corrected equations and changed compartments, I get this —
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]';
%Note: x,s and p are combined (in separate columns) to form c.
c = [s x p];
theta0 = [rand(4,1); s(1); x(1); p(1)];
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c,zeros(1,numel(theta0)));
fprintf(1,'\tRate Constants:\n')
for k1 = 1:numel(theta)
fprintf(1, '\t\tTheta(%2d) = %8.5f\n', k1, theta(k1))
end
tv = linspace(min(t), max(t));
Cfit = kinetics(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','best')
function C=kinetics(theta,t)
c0=theta(5:7);
[T,Cv]=ode45(@DifEq,t,c0);
%
function dC=DifEq(t,c)
dcdt=zeros(3,1);
dcdt(1)= theta(1).*c(1).*c(2)./(theta(2)+c(1));
dcdt(2)= -(1./(theta(3))).*theta(1).*c(1).*c(2)./(theta(2)+c(1));
dcdt(3)= (1./(theta(4))).*theta(1).*c(1).*c(2)./(theta(2)+c(1));
dC=dcdt;
end
C=Cv;
end
I suspect that the new compartment assignments:
c = [s x p];
are incorrect.
Please experiment with it to get the differential equations to be correct, and the appropriate compartment assignments to be correct. The compartments in order are the column vectors in the ‘c’ matrix.
I will do my best to help you get this to work, however I am limited by not knowing the equations you want to code in the context of the compartments you want to model.
.
Star Strider
on 15 Oct 2022
As always, my pleasure!
‘Is 7 due to the fact that I have 4 unknowns parameters (theta(1), theta(2), theta(3) and theta(4)) and 3 concentrations in c? In this case it appears the last 3 theta values in the estimated theta are the estimated initial concentrations for the provided data c?’
Yes.
‘If I already know the initial concentrations c0, do I still have to estimate them? If I did not want to estimate them then does theta becomes rand(4,1)?’
For reasons not entirely clear to me, the code works best if the initial concentrations are parameters to be estimated. I generally use the initial concentrations as the initial estimates for those parameters. (I suspect that is because they are data like all the rest, and are subjkect to some measurement variation.)
‘So can I just use the actual unknown parameters theta(1) to theta(4) in the kinetic function instead if I already know my c0 and not the estimated c0’s?’
Yes, although I use the initial concentrations as the initial estimates for those parameters when I estimate them as parameters themselves. The estimated initial concentrations are generally close to the initial concentrations in the data.
.
Learning
on 15 Oct 2022
Star Strider
on 15 Oct 2022
Thank you!
1. It originally was:
c = [x s p];
I changed it because in your Comment you wrote: ‘I made umax as K(1), ks as K(2), Yxs as K(3), and Yps as K(4), s as c(1), and x as c(2).’
2. That is your call. If the new data matches the earlier data, then that would work. If not (for example if the temperature or other conditions were changed), the system parameters would likely need to be estimated from the new data.
3. That depends on whether you are estimating new data. If you are using the previously estimated parameters with the new data, and the system conditions are similar, yes.
4. The residuals can be calculated from the original data as:
resid = c - kinetics(theta, t)
That will produce a matrix of the residuals the same size as ‘c’ that you can then do some basic stitistics with.
You can use nlparci to get confidence intervals on the parameters of the‘kinetics’ function, however getting confidence intervals on the fit are (to the best of my knowledge) not possible with a function such as ‘kinetics’ because when I last looked (and just now checked the fitnlm documentation to be sure), the Statistics and Machine Learning Toolbox regression functions did not allow error estimation of functions with multiple outputs. I do not have the Curve Fitting Toolbox, however checking the documentation for fit just now, it apparently has the same restriction.
.
Learning
on 16 Oct 2022
Star Strider
on 16 Oct 2022
As always, my pleasure!
Learning
on 24 Nov 2022
Star Strider
on 24 Nov 2022
Happy Thanksgiving to you, too!
In all likelihood, the reason for the the code taking a long time is likely that the differential equations are ‘stiff’. One reason for this is the parameters differing by several orders of magnitude. Using ode15s (or other stiff solvers such as ode23s) should significantly shorten the time required to integrate them.
There could of course be other reasons (I’d have to have all the information — code and data — to be certain), however this is the most likely possibility in my experience.
.
Learning
on 3 Dec 2022
Star Strider
on 3 Dec 2022
As always, my pleasure!
Updating the parameter estimates in real time is not something I’ve considered. I’ve not done real time parameter estimation in decades.
I’m not sure that a real time application of this is possible, however it depends on what ‘real time’ actually means. If there is time between samples to numerically integrate the differential equations, and good estimates of the existing paramerers are known, then the current approach could be used without alteration, using the existing estimates for the initial estimates for parameter estimates with the additional data.
Otherwise, a Kalman filter approach could be possible, however I’ve little experience with Kalman filters so I don’t know if one could be used in this instance. There are a large number of Kalman filter functions in MATLAB that might be possible to use. One such is kalman, and there are more than 50 others in the online documentation that come up with a search for ‘kalman’.
More Answers (1)
Works for me.
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.7662 0.7869]';
c = [x s p];
theta0=[1;1;1;1;1;1];
options = optimset('MaxFunEvals',100000,'MaxIter',100000);
[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c,[],[],options);
theta
C = kinetics(theta,t)
hold on
plot(t,c,'o')
plot(t,C)
hold off
function C=kinetics(theta,t)
%c0=[0.1;9;0.1];
c0 = [0.0904;9.0115;0.0151];
[T,Cv]=ode45(@DifEq,t,c0);
%
function dC=DifEq(t,c)
dcdt=zeros(3,1);
dcdt(1)=-theta(1).*c(1).*c(2)/(theta(2)+c(1));
dcdt(2)= -(1/theta(3)).*theta(1).*c(1).*c(2)/(theta(2)+c(1));
dcdt(3)= (1/theta(4)).*theta(1).*c(1).*c(2)/(theta(2)+c(1));
dC=dcdt;
end
C=Cv;
end
Categories
Find more on Systems of Nonlinear Equations in Help Center and File Exchange
Community Treasure Hunt
Find the treasures in MATLAB Central and discover how the community can help you!
Start Hunting!
