生化反应动力学建模报告模板.docx
《生化反应动力学建模报告模板.docx》由会员分享,可在线阅读,更多相关《生化反应动力学建模报告模板.docx(15页珍藏版)》请在冰点文库上搜索。
生化反应动力学建模报告模板
生化反应动力学建模报告
题目:
在某分批发酵过程中测得不同时刻的细胞浓度及底物浓度如下表所示,假设其底物消耗动力学符合
(Yx/s<0.3,m<4h-1),试分别用Tessier方程和Contois方程建立其细胞生长动力学模型,并判断哪个模型更能描述该发酵过程,最终建立细胞生长及底物消耗动力学模型。
2
时间/h
3
10
17
24
31
38
45
52
59
66
细胞浓度(g/L)
0.050
0.076
0.075
0.17
0.26
0.39
0.57
0.83
1.57
2.00
底物浓度(g/L)
35.00
34.70
34.24
33.76
27.81
26.32
24.03
20.57
15.51
8.59
1、模型的求解:
1.1根据题目的表中底数据得到曲线图如下
1.2底物抑制建模
根据底物抑制可建立方程为:
将上式化简得到:
里面有4个参数
,
,
,
,考虑用龙格-库塔法解出这四个参数值,利用MATLAB软件求解参数。
求得:
Um=1.003538876503289e-001
Ks=8.232639071826206e-002
Yx/s=2.293097184299229e-001
m=8.674713815945129e-001
4.352235520691609e-003
对以上参数根据初始点条件:
时间t=5的时候求解其微分方程,用4阶的龙格-库塔法,在MATLAB用ode45可求:
[t,f]=ode45(@fun1,t,f0);
其中fun1为带有微分方程的M文件。
计算误差:
test=[Cx;Cs]';
b1=f-test;
b2=b1./test;
avererror=sum(sum(abs(b2)))/11
求得误差:
平均误差=
4.352235520691609e-003
注:
程序见附页。
1.3底物抑制方程建模
根据Contois方程
可建立以下方程
化简得到
同理求得:
对以上参数可以用ode45求得:
同理,可以求得:
平均误差=3.809e-003
注:
程序见附页。
2.模型的对比
由以上的两个误差可以知道:
Tessier方程所拟合后的参数模型的准确度比较高。
问题与讨论:
(以Tessier方程为例)
不同的上限(m、Yx/s、Ks、μm)初值对模型的影响,将赋值的结果给制成表,如下所示:
项目
k=[1e61e60.34]
k=[1e21e20.34]
k=[10100.34]
k=[110.34]
m
0.7921597
0.7921597
0.7921597
0.7600367
Yx/s
0.2937506
0.2937506
0.2937506
0.2536834
Ks
0.3760055
0.3760055
0.3760055
0.3793838
μm
0.06009261
0.6009261
0.6009261
0.6009230
误差
3.625e-003
3.625e-003
3.625e-003
3.646e-003
从表中的数据可以看出,在上限值对结果没有影响,但一旦接近于最优解时,m、Yx/s、Ks、μm及误差的精度就降低了。
附录:
(1)Tessier方程
主程序:
%Tessier主程序
clc;
clear;
k=[1e61e60.34];
result1(k);
求解参数程序:
functiony=result1(k)
%全局参数
Um=k
(1);
Ks=k
(2);
Yxs=k(3);
m=k(4);
%
formatlonge;
%初值的设置
k0=[1e-11e-11e-21e-2];
%范围的设置
Lb=[1e-41e-41e-41e-4];%下限
Ub=[UmKsYxsm];%上限
%优化参数设置
options=optimset('largescale','off','display','final','tolx',1e-8,'tolfun',1e-8,'MaxFunEvals',6000);
%求解四个参数
[p,fval,exitflag,output]=fmincon(@objfun1,k0,[],[],[],[],Lb,Ub,[],options);
y=p;
disp('fval=');
disp(fval);
disp('exitflag=');
disp(exitflag);
disp('output=');
disp(output);
disp('revelenterror=');
error=error1(p);
disp(error);
优化所调用的函数objfun1:
functiondata=objfun1(k)
%定义全局参数
globalUm;
globalKs;
globalYxs;
globalm;
Um=k
(1);
Ks=k
(2);
Yxs=k(3);
m=k(4);
%时间
t=[491419242934394449];
%细胞浓度实验值
Cx=[0.030.0480.0750.120.190.290.450.681.011.41];
%底物浓度实验值
Cs=[30.0029.7629.3728.7627.8126.3224.0320.5715.518.59];
%Cx,Cs初值
f0=[0.03030.00];
%用ODE45求解微分方程得出Cx,Cs计算值
[t,f]=ode45(@fun1,t,f0);
%计算值和实验值的残差平方和
tp=[Cx;Cs]';
a1=f-tp;
a2=a1./tp;
data=sum(sum(a2.^2));
微分方程的函数fun1:
functiondf=fun1(t,f)
globalUm;
globalKs;
globalYxs;
globalm;
%Cx=>f
(1),Cs=>f
(2)
%dCx/dt=>df
(1),dCs/dt=>df
(2),
%微分方程
df=zeros(2,1);
df
(1)=Um*(1-exp(-Ks*f
(2)))*f
(1);
df
(2)=(-1)*df
(1)/Yxs-m*f
(1);
求解误差以及拟合曲线的函数error1
functionavererror=error1(k)
globalUm;
globalKs;
globalYxs;
globalm;
Um=k
(1);
Ks=k
(2);
Yxs=k(3);
m=k(4);
%时间
t=[491419242934394449];
%细胞浓度实验值
Cx=[0.030.0480.0750.120.190.290.450.681.011.41];
%底物浓度实验值
Cs=[30.0029.7629.3728.7627.8126.3224.0320.5715.518.59];
%Cx,Cs初值
f0=[0.03030.00];
%用ODE45求解微分方程得出Cx,Cs计算值
[t,f]=ode45(@fun1,t,f0);
%平均误差
test=[Cx;Cs]';
b1=f-test;
b2=b1./test;
avererror=sum(sum(abs(b2)))/22;
disp('平均误差=');
disp(avererror);
disp('CX=');
disp(f(:
1));
disp('CS=');
disp(f(:
2));
disp('Um=');
disp(Um);
disp('Ks=');
disp(Ks);
disp('Yx/s=');
disp(Yxs);
disp('m=');
disp(m);
figure
(1);
plot(t,Cs,t,Cx,t,f,'o');
figure
(2);
plot(t,Cs,t,Cx);
(2)Contois方程
主程序:
%Tessier主程序
clc;
clear;
k=[1e61e60.34];
result2(k);
求解参数程序:
functiony=result2(k)
%全局参数
Um=k
(1);
Ks=k
(2);
Yxs=k(3);
m=k(4);
%
formatlonge;
%初值的设置
k0=[1e-11e21e-21e-2];
%范围的设置
Lb=[1e-41e-41e-41e-4];%下限
Ub=[UmKsYxsm];%上限
%优化参数设置
options=optimset('largescale','off','display','final','tolx',1e-8,'tolfun',1e-8,'MaxFunEvals',6000);
%求解四个参数
[p,fval,exitflag,output]=fmincon(@objfun2,k0,[],[],[],[],Lb,Ub,[],options);
y=p;
disp('fval=');
disp(fval);
disp('exitflag=');
disp(exitflag);
disp('output=');
disp(output);
disp('revelenterror=');
error=error2(p);
disp(error);
优化所调用的函数objfun2:
functiondata=objfun2(k)
%定义全局参数
globalUm;
globalKs;
globalYxs;
globalm;
Um=k
(1);
Ks=k
(2);
Yxs=k(3);
m=k(4);
%时间
t=[491419242934394449];
%细胞浓度实验值
Cx=[0.030.0480.0750.120.190.290.450.681.011.41];
%底物浓度实验值
Cs=[30.0029.7629.3728.7627.8126.3224.0320.5715.518.59];
%Cx,Cs初值
f0=[0.03030.00];
%用ODE45求解微分方程得出Cx,Cs计算值
[t,f]=ode45(@fun2,t,f0);
%计算值和实验值的残差平方和
tp=[Cx;Cs]';
a1=f-tp;
a2=a1./tp;
data=sum(sum(a2.^2));
微分方程的函数fun2:
functiondf=fun2(t,f)
globalUm;
globalKs;
globalYxs;
globalm;
%Cx=>f
(1),Cs=>f
(2)
%dCx/dt=>df
(1),dCs/dt=>df
(2),
%微分方程
df=zeros(2,1);
df
(1)=Um*f
(2)/(Ks*f
(1)+f
(2))*f
(1);
df
(2)=(-1)*df
(1)/Yxs-m*f
(1);
求解误差以及拟合曲线的函数error2
functionavererror=error2(k)
globalUm;
globalKs;
globalYxs;
globalm;
Um=k
(1);
Ks=k
(2);
Yxs=k(3);
m=k(4);
%时间
t=[491419242934394449];
%细胞浓度实验值
Cx=[0.030.0480.0750.120.190.290.450.681.011.41];
%底物浓度实验值
Cs=[30.0029.7629.3728.7627.8126.3224.0320.5715.518.59];
%Cx,Cs初值
f0=[0.03030.00];
%用ODE45求解微分方程得出Cx,Cs计算值
[t,f]=ode45(@fun2,t,f0);
%平均误差
test=[Cx;Cs]';
b1=f-test;
b2=b1./test;
avererror=sum(sum(abs(b2)))/22;
disp('平均误差=');
disp(avererror);
disp('CX=');
disp(f(:
1));
disp('CS=');
disp(f(:
2));
disp('Um=');
disp(Um);
disp('Ks=');
disp(Ks);
disp('Yx/s=');
disp(Yxs);
disp('m=');
disp(m);
figure
(1);
plot(t,Cs,t,Cx,t,f,'o');
figure
(2);
plot(t,Cs,t,Cx);