I want these kinetic models to fit in MATLAB but I don't know how to proceed. I am in real worry. If anyone can solve/code, it would be of great help. Thanks a million!!

Question

Shubhangi Agrawal 2020-8-21

0
链接

此问题的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/582614-i-want-these-kinetic-models-to-fit-in-matlab-but-i-don-t-know-how-to-proceed-i-am-in-real-worry-i

评论： Shubhangi Agrawal 2021-4-15

采纳的回答： Alan Stevens

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

请先登录，再进行评论。

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

Answer 1

Alan Stevens 2020-8-23

0
链接

此回答的直接链接

https://ww2.mathworks.cn/matlabcentral/answers/582614-i-want-these-kinetic-models-to-fit-in-matlab-but-i-don-t-know-how-to-proceed-i-am-in-real-worry-i#answer_483782

编辑：Alan Stevens 2020-8-23

在 MATLAB Online 中打开

You need to define Haldane as follows (your text-book expression is misleading)

Haldane = @(S) 0.478*S / ((1.137 + S) .* (1 + (S/1.139)));

Then you get:

which matches the line on your text-book graph.

32 个评论
显示 30更早的评论隐藏 30更早的评论

Shubhangi Agrawal 2021-1-17

I tried the way you told me. But it's not working.

Here is the code I am attaching. Could you please try it once in your MATLAB and setting the y-axis as 0 2 4 6 8?

clc;

clear;

%model parameters from research paper

yxs=0.20;

mumax=0.34;

ks=4.21;

sm=23.41;

n=4.12;

%experimental extracted values from webplot digitalizer

%TA time of substrate consumption

%A substrate

%TB time of biomass growth

%B biomass

TA=[0.0;1.86;3.15;4.59;5.59;7.02;8.17;9.17;10.03;11.18;11.75;12.04;12.32;12.75;13.0;13.33;13.76;14.04;14.76;15.19];

A=[3.00;2.99;2.92;2.8;2.67;2.49;2.33;2.19;1.99;1.81;1.65;1.45;1.3;1.14;0.98;0.78;0.62;0.48;0.3;0.08;];

hold on % hold on to plot in one graph

plot(TA,A,'r s')

TB=[0;1;1.43;2.29;3.44;4.3;5.16;6.02;6.88;7.74;8.74;9.6;10.6;11.32;12.47;13.33;14.47;14.61;15.04;15.19];

B=[0.04;0.11;0.14;0.18;0.25;0.25;0.27;0.33;0.33;0.33;0.35;0.35;0.38;0.4;0.46;0.48;0.49;0.55;0.58;0.6];

hold on

plot(TB,B,'b o')

param = [yxs;mumax;ks;sm;n];

[t,X] = ode45(@Luong,linspace(0,96,2),[0.2 3],[],param);

hold on

plot(t,X)

set(gca, 'XTicklabel', [0 12 24 36 48 60 72 84 96])

title('Batch Process Simulation')

xlabel('Time(hr)')

ylabel('Concentation in g/l')

legend('S(Substrate Experimental)','X(Biomass Experimental)','X(Biomass predicted)','S(Substrate predicted)')

function dX=Luong(t,x,param)

x1 = x(1);

s = x(2);

yxs = param(1);

mumax = param(2);

ks = param(3);

sm = param(4);

n = param(5);

mu = (mumax* s /(ks + s) )*((1 -(s/sm))^n);

dx1 = mu * x1;

ds = (-(1/yxs)* dx1);

dX = [dx1 ds]'

end

Alan Stevens 2021-1-17

编辑：Alan Stevens 2021-1-17

在 MATLAB Online 中打开

Of course, one just needs to use the set(gca,...) command!

Try the following

%model parameters from research paper
yxs=0.20;
mumax=0.34;
ks=4.21;
sm=23.41;
n=4.12;
%experimental extracted values from webplot digitalizer
%TA time of substrate consumption
%A substrate
%TB time of biomass growth
%B biomass
TA=[0.0;1.86;3.15;4.59;5.59;7.02;8.17;9.17;10.03;11.18;11.75;12.04;12.32;12.75;13.0;13.33;13.76;14.04;14.76;15.19];
A=[3.00;2.99;2.92;2.8;2.67;2.49;2.33;2.19;1.99;1.81;1.65;1.45;1.3;1.14;0.98;0.78;0.62;0.48;0.3;0.08;];
hold on % hold on to plot in one graph
plot(TA,A,'r s')
TB=[0;1;1.43;2.29;3.44;4.3;5.16;6.02;6.88;7.74;8.74;9.6;10.6;11.32;12.47;13.33;14.47;14.61;15.04;15.19];
B=[0.04;0.11;0.14;0.18;0.25;0.25;0.27;0.33;0.33;0.33;0.35;0.35;0.38;0.4;0.46;0.48;0.49;0.55;0.58;0.6];
hold on
plot(TB,B,'b o')
param = [yxs;mumax;ks;sm;n];
[t,X] = ode45(@Luong,linspace(0,96,2),[0.2 3],[],param);
hold on
plot(t,X), grid
axis([0 100 0 8])
set(gca,'YTick',[0 2 4 6 8])
set(gca,'YTickLabel',[0 2 4 6 8])
set(gca,'XTick',[0 12 24 36 48 60 72 84 96])
set(gca, 'XTicklabel', [0 12 24 36 48 60 72 84 96])
title('Batch Process Simulation')
xlabel('Time(hr)')
ylabel('Concentation in g/l')
legend('S(Substrate Experimental)','X(Biomass Experimental)','X(Biomass predicted)','S(Substrate predicted)')
function dX=Luong(~,x,param)
x1 = x(1);
s = x(2);
yxs = param(1);
mumax = param(2);
ks = param(3);
sm = param(4);
n = param(5);
mu = (mumax* s /(ks + s) )*((1 -(s/sm))^n);
dx1 = mu * x1;
ds = (-(1/yxs)* dx1);
dX = [dx1 ds]';
end

This produces

Shubhangi Agrawal 2021-1-25

I am fitting these curves in my MATLAB.

with the following equations-

But, My curves are not exactly matching with the above figure. It's appearing like-

You can see here that the curves are not matching their y-axis.

Could you help me to find where the error lies.

I am attaching the code which I used to plot.

---------------------------------------------------------------------------------------------------------

% implement YanoandKoga1 Model as a Matlab anonymous function

YanoandKoga1 = @(S) 0.51*S./(1.413+S+(S.^2/1250)+(S.^3/4.561));

% calculate mu

mu = YanoandKoga1(S);

hold on

% plot result

plot(S,mu,'k --'), xlabel('S (g/l)'), ylabel('\mu (1/h)');

% implement YanoandKoga2 Model as a Matlab anonymous function

YanoandKoga2= @(S) 0.41*S./(5.23+S+(S.^3/2.64));

% calculate mu

mu = YanoandKoga2(S);

hold on

% plot result

plot(S,mu,'k-'), xlabel('S (g/l)'), ylabel('\mu (1/h)');

% implement AlagappanandCowan Model as a Matlab anonymous function

AlagappanandCowan= @(S)(0.82*S./(4.366+S+(S.^2/2.246)))-(0.004379*(S+2.985));

% calculate mu

mu = AlagappanandCowan(S);

hold on

% plot result

plot(S,mu,'k-o'), xlabel('S (g/l)'), ylabel('\mu (1/h)');

% implement WaymanandTseng Model as a Matlab anonymous function

WaymanandTseng= @(S)(0.3164*S./(0.7346+S))-((0.01618)*(S+0.2886));

% calculate mu

mu = WaymanandTseng(S);

hold on

% plot result

plot(S,mu,'k -.'), xlabel('S (g/l)'), ylabel('\mu (1/h)'), legend('YanoandKoga1','YanoandKoga2', 'AlagappanandCowan','WaymanandTseng');

axis([0 20 0 0.3])

set(gca,'YTick',[0 0.05 .1 0.15 0.2 0.25 0.3 ])

set(gca,'YTickLabel',[0 0.05 .1 0.15 0.2 0.25 0.3 ])

set(gca,'XTick',[0 2 4 6 8 10 12 14 16 18 20 ])

set(gca, 'XTicklabel', [0 2 4 6 8 10 12 14 16 18 20 ])

title('\mu vs S')

Shubhangi Agrawal 2021-3-7

Hi Alan,

I am plotting these two graphs in one graph in MATLAB.

I am getting my output as-

I want X and Y axis to be exactly similar as my textbook graphs i.e 0 6 12 18 24 30 (X-axis) and 0 2 4 6 8(Y-axis) but not able to do so.

The code I wrote is--

------------------------------------------------------------------------------------------------------------

yxs=0.228;% biomass yield

mumax=0.34;%maximum specific growth rate

Sin=1;%input substate concentration

ks= 4.21;

Sm=23.41 ;

n = 4.12;

param = [yxs mumax Sin ks Sm n];

Opt = odeset('Events', @myEventsFcn);

[t1,X1]= ode45(@(t,y)fed_batch(t,y,param),linspace(0,36,2),[0.2 3 ],Opt);% matlab ode45 function

[t2,X2]= ode45(@(t,y)fed_batch(t,y,param),linspace(t1(end),64,2),[X1(end,1) 6.5],Opt);% matlab ode45 function

[t3,X3]= ode45(@(t,y)fed_batch(t,y,param),linspace(t2(end),84,2),[X2(end,1) 6.5],Opt);% matlab ode45 function

%[t4,X4]= ode45(@(t,y)fed_batch(t,y,param),linspace(t3(end),1000,2),[X3(end,1) 100 X3(end,3:4)],Opt);% matlab ode45 function

%[t5,X5]= ode45(@(t,y)fed_batch(t,y,param),linspace(t4(end),1000,2),[X4(end,1) 100 X4(end,3:4)],Opt);% matlab ode45 function

%[t6,X6]= ode45(@(t,y)fed_batch(t,y,param),linspace(t5(end),1000,2),[X5(end,1) 100 X5(end,3:4)],Opt);% matlab ode45 function

X = [X1;X2;X3];

t = [t1;t2;t3];

axis([0 40 0 8])

set(gca,'YTick',[0 2 4 6 8 ])

set(gca,'YTickLabel',[0 2 4 6 8 ])

set(gca,'XTick',[0 6 12 18 24 36 ])

set(gca, 'XTicklabel', [0 6 12 18 24 36])

plot(t,X)

title('Concentration versus Time')

xlabel('Time(h)')

ylabel('Concentration (g/L)')

legend('Biomass','Substrate','Location','northwest')

grid on

%TA time of substrate consumption

%A substrate

%TB time of biomass growth

%B biomass

function dX= fed_batch(t,y, param)

x1 = y(1);

s = y(2);

yxs = param(1);

mumax=param(2);

Sin=param(3);

ks = param(4);

sm = param(5);

n = param(6);

mu = (mumax* s /(ks + s) )*((1 -(s/sm))^n);% luong model to estimate specific growth rate

%D = mu;%assuming at steady state, dilution rate equal to mu

%dx1 = (( mu *x1) - (x1* (f/v)));% biomass model equation

dx1 = (mu.*x1);

ds = - (dx1*(1/ yxs));%substatre model equation

dX = [dx1 ds]';

end

function [value, isterminal, direction] = myEventsFcn(t, X,param)

value = X(2) <= 0;

isterminal = 1;

direction = [];

end

----------------------------------------------------------------------------------------------------------------------------------------

I request you to please help me out in the problem. I am really feeling helpless.

Thanks in advance.

Alan Stevens 2021-3-7

在 MATLAB Online 中打开

The following gives something like the text-book values (though not exactly the same), but, in order to achieve this you need to use the Luong constants.

yxs=0.228;% biomass yield
% Luong values
mumax = 0.2407;
ks = 0.5995;
Sm = 20.17; 
n = 2.962;
param = [yxs mumax Sin ks Sm n];
Opt = odeset('Events', @myEventsFcn);
[t1,X1]= ode45(@(t,y)fed_batch(t,y,param),linspace(0,13),[0.2 3 ],Opt);% matlab ode45 function
[t2,X2]= ode45(@(t,y)fed_batch(t,y,param),linspace(t1(end),23),[X1(end,1) 6.5],Opt);% matlab ode45 function
[t3,X3]= ode45(@(t,y)fed_batch(t,y,param),linspace(t2(end),30),[X2(end,1) 6.5],Opt);% matlab ode45 function
X = [X1;X2;X3];
t = [t1;t2;t3];
plot(t,X)  %%%%%% Put plot command before axis, YTick etc commands
axis([0 30 0 8])
set(gca,'YTick',[0 2 4 6 8 ])
set(gca,'YTickLabel',[0 2 4 6 8  ])
set(gca,'XTick',[0 6 12 18 24 30 ])
set(gca, 'XTicklabel', [0 6 12 18 24 30])
title('Concentration versus Time')
xlabel('Time(h)')
ylabel('Concentration (g/L)')
legend('Biomass','Substrate','Location','northwest')
grid on
function dX= fed_batch(~,y, param)
x1 = y(1); 
s = y(2); 
yxs = param(1);
mumax=param(2);
%Sin=param(3);
ks = param(4);
sm = param(5);
n = param(6);
mu = (mumax* s /(ks + s) )*((1 -(s/sm))^n);% luong model to estimate specific growth rate
dx1 = (mu.*x1);
ds =  - (dx1*(1/ yxs));%substatre model equation
dX = [dx1 ds]';
end
function [value, isterminal, direction] = myEventsFcn(~, X,~)
value = X(2) <= 0;
isterminal = 1;
direction = [];
end

This results in

Shubhangi Agrawal 2021-3-9

Also, when I am trying for another graph i.e.-

I am getting the output as-

Could you please tell me why I am not getting initial 4-5 curves (till 36 h) exactly as my textbook graph (only for caffeine).

The code I used is--

-----------------------------------------------------------------------------------------------------------------------------------

yxs=0.228;% biomass yield

% Luong values

mumax = 0.2407;

ks = 0.5995;

Sm = 20.17;

n = 2.962;

Sin=3;

%TA time of substrate consumption

%A substrate

%TB time of biomass growth

%B biomass

param = [yxs mumax Sin ks Sm n];

Opt = odeset('Events', @myEventsFcn);

[t1,X1]= ode45(@(t,y)fed_batch(t,y,param),linspace(0,15),[0.2 3 ],Opt);% matlab ode45 function

[t2,X2]= ode45(@(t,y)fed_batch(t,y,param),linspace(t1(end),21),[X1(end,1) 6.5],Opt);% matlab ode45 function

[t3,X3]= ode45(@(t,y)fed_batch(t,y,param),linspace(t2(end),26),[X2(end,1) 6.5],Opt);% matlab ode45 function

[t4,X4]= ode45(@(t,y)fed_batch(t,y,param),linspace(t3(end),30),[X3(end,1) 6.5],Opt);% matlab ode45 function

[t5,X5]= ode45(@(t,y)fed_batch(t,y,param),linspace(t4(end),33),[X4(end,1) 6.5],Opt);% matlab ode45 function

[t6,X6]= ode45(@(t,y)fed_batch(t,y,param),linspace(t5(end),36),[X5(end,1) 6.5],Opt);% matlab ode45 function

[t7,X7]= ode45(@(t,y)fed_batch(t,y,param),linspace(t6(end),39),[X6(end,1) 6.5],Opt);% matlab ode45 function

[t8,X8]= ode45(@(t,y)fed_batch(t,y,param),linspace(t7(end),42),[X7(end,1) 6.5],Opt);% matlab ode45 function

[t9,X9]= ode45(@(t,y)fed_batch(t,y,param),linspace(t8(end),45),[X8(end,1) 6.5],Opt);% matlab ode45 function

[t10,X10]= ode45(@(t,y)fed_batch(t,y,param),linspace(t9(end),48),[X9(end,1) 6.5],Opt);% matlab ode45 function

[t11,X11]= ode45(@(t,y)fed_batch(t,y,param),linspace(t10(end),51),[X10(end,1) 6.5],Opt);% matlab ode45 function

[t12,X12]= ode45(@(t,y)fed_batch(t,y,param),linspace(t11(end),53),[X11(end,1) 6.5],Opt);% matlab ode45 function

[t13,X13]= ode45(@(t,y)fed_batch(t,y,param),linspace(t12(end),54),[X12(end,1) 6.5],Opt);% matlab ode45 function

[t14,X14]= ode45(@(t,y)fed_batch(t,y,param),linspace(t13(end),57),[X13(end,1) 6.5],Opt);% matlab ode45 function

[t15,X15]= ode45(@(t,y)fed_batch(t,y,param),linspace(t14(end),59),[X14(end,1) 6.5],Opt);% matlab ode45 function

[t16,X16]= ode45(@(t,y)fed_batch(t,y,param),linspace(t15(end),60),[X15(end,1) 6.5],Opt);% matlab ode45 function

[t17,X17]= ode45(@(t,y)fed_batch(t,y,param),linspace(t16(end),62),[X16(end,1) 6.5],Opt);% matlab ode45 function

[t18,X18]= ode45(@(t,y)fed_batch(t,y,param),linspace(t17(end),64),[X17(end,1) 6.5],Opt);% matlab ode45 function

[t19,X19]= ode45(@(t,y)fed_batch(t,y,param),linspace(t18(end),66),[X18(end,1) 6.5],Opt);% matlab ode45 function

[t20,X20]= ode45(@(t,y)fed_batch(t,y,param),linspace(t19(end),68),[X19(end,1) 6.5],Opt);% matlab ode45 function

[t21,X21]= ode45(@(t,y)fed_batch(t,y,param),linspace(t20(end),70),[X20(end,1) 6.5],Opt);% matlab ode45 function

X = [X1;X2;X3;X4;X5;X6;X7;X8;X9;X10;X11;X12;X13;X14;X15;X16;X17;X18;X19;X20;X21];

t = [t1;t2;t3;t4;t5;t6;t7;t8;t9;t10;t11;t12;t13;t14;t15;t16;t17;t18;t19;t20;t21];

plot(t,X) %%%%%% Put plot command before axis, YTick etc commands

axis([0 72 0 10])

set(gca,'YTick',[0 2.5 5.0 7.5 10 ])

set(gca,'YTickLabel',[0 2.5 5.0 7.5 10 ])

set(gca,'XTick',[0 12 24 36 48 60 72 ])

set(gca, 'XTicklabel', [0 12 24 36 48 60 72])

title('Concentration versus Time')

xlabel('Time(h)')

ylabel('Concentration (g/L)')

legend('Biomass', 'Caffeine','Location','northwest')

function dX= fed_batch(~,y, param)

x1 = y(1);

s = y(2);

yxs = param(1);

mumax=param(2);

%Sin=param(3);

ks = param(4);

sm = param(5);

n = param(6);

mu = (mumax* s /(ks + s) )*((1 -(s/sm))^n);% luong model to estimate specific growth rate

dx1 = (mu.*x1);

ds = - (dx1*(1/ yxs));%substatre model equation

dX = [dx1 ds]';

end

function [value, isterminal, direction] = myEventsFcn(~, X,~)

value = X(2) <= 0;

isterminal = 1;

direction = [];

end

Steven Lord 2021-4-14

As an FYI, you could use the xticks and xticklabels functions (and the obvious equivalent for the Y axes) as long as you're using release R2016b or later. You could put those calls in the same place as your other labeling calls (xlabel and ylabel.)

Shubhangi Agrawal 2021-4-15