Kinetic Models (Monod Equations)

Question

Learning 2022-10-9

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1820635-kinetic-models-monod-equations

评论： Star Strider 2022-10-12

采纳的回答： Star Strider

Hello,

I saw the equations below in litrature. This is a process model for cultivation/fermentation process. Question is, being new to Matlab, how do I code this and fit it to experiemtal data?
Can these equations be implemented as is in Simulink and the unknown parameters estimated based on the test/experiment data? like can Simulink solve these equations? I am not sure how to implement this and would any guiance you can provide.
Example video from MathWorks is How to Estimate Model Parameters from Test Data with Simulink Video - MATLAB & Simulink (mathworks.com).

Thanks!

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

Ghazwan 2022-10-9

Hello! Yes, you can use both Matlab and Simulink. With Matlab, you can use ode45 and fminsearch to fit the equations to the experimental and find the parameters. Basically, the optimization tool functions.

请先登录，再进行评论。

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

Answer 1

Star Strider 2022-10-9

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/1820635-kinetic-models-monod-equations#answer_1069960

See: Monod kinetics and curve fitting

A somewhat improved version of that is in: Parameter Estimation for a System of Differential Equations

.

17 个评论
显示 15更早的评论隐藏 15更早的评论

Star Strider 2022-10-9

As always, my pleasure!

It might be necessary for you to re-code your differential equations if they differ from those in my original code (that simply illustrates a way to solve the problem).

The code in my first reference deals specifically with Monod kinetics. My code in the second reference demonstrates an improved version of the code and displays the fittted parameter responses with the data. If you want to make predictions for times beyond the original tme vector, run the same code with the estimated parameters with a new time vector that extends the original time vector. That would require saving the estimated parameters and using them with the existing code, with a new time vector.

If you output the result of the differential equation integration to a structure instead of the way I have the results, you can use the deval function with a new time vector to plot predictions beyone the existing time vector. That would require a slight re-write of my existing code to display the result with the existing time vector and data, however it is certainly an option.

Star Strider 2022-10-10

As always, my pleasure!

The ‘[1;1;1;1]’ are simply starting parameter estimates. They can be anything reasonable, although they should have a lower bound of zeros to prevent them from becoming negative (since kinetic parameters are always positive by defintion in most models). My code most recently estimates the initial concentrations as well as the kinetic constants, since that produces a better fit to the data. They should also be bounded as positive.

I am not certain what you mean by ‘validating it with different experimental data’. It is likely best to estimate the parameters for each set of data, although it is certainly possible to use a given set of parameters with different data. To use the same parameters, use the new data with a given set of parameters and see how well they estimate the new data. I would need a specific example (and data) to expand on this in order to understand what you want to do.

It is certainly possible to use a different time vector with a given set of estimated parameters, and the code does that in order to produce a smooth curve for the plots at the end. That time vector does not have to end when the original time vector ends, and can easily be longer. This extrapolates beyond the region-of-fit, and while that is not recommended (because that amounts to ‘guessing’ what the model will do beyond the actual data), it is certianly possible.

The original parameter estimation uses only the data at the time instants specified, and those can be produced at the end of the parameter estimation process. Interpolating data at new times within the original interval is permitted (as opposed to extrapolating beyond the original time interval that can be done but is not recommended), and the code does that to draw the plots at the end.

I can provide more specific comments with specific data. I will then have a better idea of what the specific issues are.

.

Star Strider 2022-10-10

在 MATLAB Online 中打开

If you already know the parameters, there is no need to do the estimation, so comment out the lsqcurvefit call. Just load them as the ‘theta’ vector (by readint a file containing them or simply pasting the vector into the code) to pass them to the function, similar to what I did after the lsqcurvefit call.

That might be something like this (using your own ‘kinetics’ code) —

function C=kinetics(theta,t)
    c0=theta(end-3:end);
    [T,Cv]=ode45(@DifEq,t,c0);
    %
    function dC=DifEq(t,c)
        dcdt=zeros(4,1);
        dcdt(1)=-theta(1).*c(1)-theta(2).*c(1);
        dcdt(2)= theta(1).*c(1)+theta(4).*c(3)-theta(3).*c(2)-theta(5).*c(2);
        dcdt(3)= theta(2).*c(1)+theta(3).*c(2)-theta(4).*c(3)+theta(6).*c(4);
        dcdt(4)= theta(5).*c(2)-theta(6).*c(4);
        dC=dcdt;
    end
    C=Cv;
end
% LOAD THE 'theta' VECTOR AND THE 'time' VECTOR (OR AT LEAST A TWO-ELEMENT VERSION OF THE 'time' VECTOR) HERE
fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(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')
end

I have not specifically tested it, however this should work.

.

Learning 2022-10-11

Sure. Thanks for your patience for this. your code is below:

function Igor_Moura

% 2016 12 03

% NOTES:

%

% 1. The ‘theta’ (parameter) argument has to be first in your

% ‘kinetics’ funciton,

% 2. You need to return ALL the values from ‘DifEq’ since you are fitting

% all the values

function C=kinetics(theta,t)

c0=[1;0;0;0];

[T,Cv]=ode45(@DifEq,t,c0);

%

function dC=DifEq(t,c)

dcdt=zeros(4,1);

dcdt(1)=-theta(1).*c(1)-theta(2).*c(1);

dcdt(2)= theta(1).*c(1)+theta(4).*c(3)-theta(3).*c(2)-theta(5).*c(2);

dcdt(3)= theta(2).*c(1)+theta(3).*c(2)-theta(4).*c(3)+theta(6).*c(4);

dcdt(4)= theta(5).*c(2)-theta(6).*c(4);

dC=dcdt;

end

C=Cv;

end

t=[0.1

0.2

0.4

0.6

0.8

1

1.5

2

3

4

5

6];

c=[0.902 0.06997 0.02463 0.00218

0.8072 0.1353 0.0482 0.008192

0.6757 0.2123 0.0864 0.0289

0.5569 0.2789 0.1063 0.06233

0.4297 0.3292 0.1476 0.09756

0.3774 0.3457 0.1485 0.1255

0.2149 0.3486 0.1821 0.2526

0.141 0.3254 0.194 0.3401

0.04921 0.2445 0.1742 0.5277

0.0178 0.1728 0.1732 0.6323

0.006431 0.1091 0.1137 0.7702

0.002595 0.08301 0.08224 0.835];

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, '\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

now if I remove write it again an remove t and c data from the body of the code, it looks like this below. I then write t and c to the workspace and I try calling it in the code, it doesn't work. One thing also is that, c0 = [1;0;0;0], is this the initial conitions for C(1), C(2), C(3) and c(4)?

c=[0.902 0.06997 0.02463 0.00218

0.8072 0.1353 0.0482 0.008192

0.6757 0.2123 0.0864 0.0289

0.5569 0.2789 0.1063 0.06233

0.4297 0.3292 0.1476 0.09756

0.3774 0.3457 0.1485 0.1255

0.2149 0.3486 0.1821 0.2526

0.141 0.3254 0.194 0.3401

0.04921 0.2445 0.1742 0.5277

0.0178 0.1728 0.1732 0.6323

0.006431 0.1091 0.1137 0.7702

0.002595 0.08301 0.08224 0.835];

t=[0.1

0.2

0.4

0.6

0.8

1

1.5

2

3

4

5

6];

function Igor_Moura

% 2016 12 03

% NOTES:

%

% 1. The ‘theta’ (parameter) argument has to be first in your

% ‘kinetics’ funciton,

% 2. You need to return ALL the values from ‘DifEq’ since you are fitting

% all the values

function C=kinetics(theta,t)

c0=[1;0;0;0];

[T,Cv]=ode45(@DifEq,t,c0);

%

function dC=DifEq(t,c)

dcdt=zeros(4,1);

dcdt(1)=-theta(1).*c(1)-theta(2).*c(1);

dcdt(2)= theta(1).*c(1)+theta(4).*c(3)-theta(3).*c(2)-theta(5).*c(2);

dcdt(3)= theta(2).*c(1)+theta(3).*c(2)-theta(4).*c(3)+theta(6).*c(4);

dcdt(4)= theta(5).*c(2)-theta(6).*c(4);

dC=dcdt;

end

C=Cv;

end

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, '\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

Star Strider 2022-10-12

在 MATLAB Online 中打开

I changed the code slightly to reflect updates to my original code, and added a lower bound of zeros to constrain all the parameters to be positive.

This example seems to use the original data, so pasting your data into it where I put ‘c’ and ‘t’ should work.

Update ‘theta0’ to reflect your differential equation system so the length is (number of parameters + number of differential equations).

This version of the code also estimates the initial conditions of the differrential equations as parameters, as well as the kinetic constants, since that generally provides a better fit to the data. The initial conditions are the last N parameters of theta0’ and ‘theta’ where N is the number of differential equations in the system. The rest of the code adapts to the data, with the plot legend matching the columns of ‘c’.

In the original version of the code, the functions and the calling code all had to be within a separate function (or different function files, however that was inconvenient to post here so I put them all in one function and then called that function to execute the code). Several versions/releases ago MATLAB allowed the functions to be in the script, providing they are all at the end of the script, so this version of my original code follows that later convention.

Try this —

c=[0.902 0.06997 0.02463 0.00218

0.8072 0.1353 0.0482 0.008192

0.6757 0.2123 0.0864 0.0289

0.5569 0.2789 0.1063 0.06233

0.4297 0.3292 0.1476 0.09756

0.3774 0.3457 0.1485 0.1255

0.2149 0.3486 0.1821 0.2526

0.141 0.3254 0.194 0.3401

0.04921 0.2445 0.1742 0.5277

0.0178 0.1728 0.1732 0.6323

0.006431 0.1091 0.1137 0.7702

0.002595 0.08301 0.08224 0.835];

t=[0.1

0.2

0.4

0.6

0.8

1

1.5

2

3

4

5

6];

theta0=rand(10,1);

[theta,Rsdnrm,Rsd,ExFlg,OptmInfo,Lmda,Jmat]=lsqcurvefit(@kinetics,theta0,t,c,zeros(size(theta0)));

Local minimum found. Optimization completed because the size of the gradient is less than the value of the optimality 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.76482 Theta( 2) = 0.23492 Theta( 3) = 0.20877 Theta( 4) = 0.49179 Theta( 5) = 0.62211 Theta( 6) = 0.00000 Theta( 7) = 0.90287 Theta( 8) = 0.07146 Theta( 9) = 0.02840 Theta(10) = 0.00000

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(7:10);

[T,Cv]=ode45(@DifEq,t,c0);

function dC=DifEq(t,c)

dcdt=zeros(4,1);

dcdt(1)=-theta(1).*c(1)-theta(2).*c(1);

dcdt(2)= theta(1).*c(1)+theta(4).*c(3)-theta(3).*c(2)-theta(5).*c(2);

dcdt(3)= theta(2).*c(1)+theta(3).*c(2)-theta(4).*c(3)+theta(6).*c(4);

dcdt(4)= theta(5).*c(2)-theta(6).*c(4);

dC=dcdt;

end

C=Cv;

end

Put your ‘c’ and ‘t’ in place of these, and provide your own version of ‘kinetics’ and it should work.

.

Learning 2022-10-12

编辑：Learning 2022-10-12

This is very helpful. You’ve done great and I am only grateful. I cannot get Cfit, theta vector and Tv to be output into the workspace. Logic is after finding the rate constants, I want to save these Cfit, tv and theta vector to the workspace so I can plot new data on it without using the Isqcurvefit (will comment out this call in the code).

But thank you brother. I appreciate it.

Star Strider 2022-10-12

My pleasure!

I do not understand the problems you are having.

请先登录，再进行评论。

Kinetic Models (Monod Equations)

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

采纳的回答

17 个评论
显示 15更早的评论隐藏 15更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

Kinetic Models (Monod Equations)

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

采纳的回答

17 个评论 显示 15更早的评论隐藏 15更早的评论

更多回答（0 个）

另请参阅

类别

标签

产品

版本

Community Treasure Hunt

WeChat

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

17 个评论
显示 15更早的评论隐藏 15更早的评论