Finding Unknown Parameters from Monod Equations

Question

Learning 2022-10-14

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1826608-finding-unknown-parameters-from-monod-equations

评论： Star Strider 2022-12-3

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

Alex Sha 2022-10-16

在 MATLAB Online 中打开

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 

请先登录，再进行评论。

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

Answer 1

Star Strider 2022-10-14

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1826608-finding-unknown-parameters-from-monod-equations#answer_1075168

在 MATLAB Online 中打开

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.

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

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.

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

Rate Constants:

for k1 = 1:numel(theta)

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

end

Theta( 1) = 432.22494 Theta( 2) = 11061.84689 Theta( 3) = 0.13367 Theta( 4) = 2.38808 Theta( 5) = 0.02007 Theta( 6) = 8.88929 Theta( 7) = 0.05209

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

Learning 2022-10-14

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!

Star Strider 2022-10-14

在 MATLAB Online 中打开

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.

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

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.

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

Rate Constants:

for k1 = 1:numel(theta)

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

end

Theta( 1) = 0.00014 Theta( 2) = 973.54767 Theta( 3) = 752.06527 Theta( 4) = 0.00002 Theta( 5) = 4.01626 Theta( 6) = 0.70152 Theta( 7) = 0.03684

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 2022-10-15

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.

.

Star Strider 2022-10-15

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 2022-12-3

Thanks for your response! My final question is how do I implement this in real time? Say I identify the parameters for the series of odes, do I need to have another set of data to estimate the parameters? Or I can fix the parameters if I know the process doesn’t change that much? Real time in the sense that assuming I have live stream data from a hardware that I want to compare with the concentrations estimated by the ode at specific times.

Thank you!

Star Strider 2022-12-3

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

请先登录，再进行评论。

Answer 2

Torsten 2022-10-14

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1826608-finding-unknown-parameters-from-monod-equations#answer_1075163

在 MATLAB Online 中打开

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

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.

theta

theta = 6×1

-3.5510 133.8605 -0.1503 -1.9775 1.0000 1.0000

C = kinetics(theta,t)

C = 12×3

0.0904 9.0115 0.0151 0.1811 8.4079 0.0610 0.2782 7.7620 0.1101 0.4889 6.3604 0.2166 0.6186 5.4979 0.2822 0.8008 4.2856 0.3743 1.0126 2.8768 0.4814 1.1491 1.9684 0.5505 1.2507 1.2924 0.6019 1.3212 0.8237 0.6375

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

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

请先登录，再进行评论。

Finding Unknown Parameters from Monod Equations

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

采纳的回答

12 个评论
显示 10更早的评论隐藏 10更早的评论

更多回答（1 个）

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

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Finding Unknown Parameters from Monod Equations

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

采纳的回答

12 个评论 显示 10更早的评论隐藏 10更早的评论

更多回答（1 个）

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

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

WeChat

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

12 个评论
显示 10更早的评论隐藏 10更早的评论

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