04_多元模型回归与分析.ppt

《04_多元模型回归与分析.ppt》由会员分享，可在线阅读，更多相关《04_多元模型回归与分析.ppt（17页珍藏版）》请在淘文阁 - 分享文档赚钱的网站上搜索。

1、1使用使用MATLAB进行参数估计示例进行参数估计示例q例例1：等温积分反应器的参数估计：等温积分反应器的参数估计在一等温积分反应器的动力学实验中，发生如下反应在一等温积分反应器的动力学实验中，发生如下反应 ：已知：当已知：当t=0时，时，CA0=1，CB0=0，CC0=0。实验数据如下表：。实验数据如下表：ti051015205080100200CAi,exp1.0000.8400.6790.5030.4250.1370.0300.0060.000CBi,exp0.0000.1550.3270.4260.4770.5710.5020.3850.145q模型：模型：因为因为CA+CB+CC=1

2、，所以描述反应器出口反应物浓度变化的模，所以描述反应器出口反应物浓度变化的模型方程只需要两个：型方程只需要两个：初始条件为：初始条件为：t=0，CA0=1，CB0=0，CC0=0。2使用使用MATLAB的的dsolve求积分求积分 s=dsolve(Dy1=-k1*y1,Dy2=k1*y1-k2*y2,y1(0)=1,y2(0)=0)s=y1:1x1 sym y2:1x1 sym s.y1 s.y2 simplify(s.y2)ans=exp(-k1*t)ans=-(-k1+k2)/(k1-k2)2*k1*exp(-k2*t)-1/(k1-k2)*k1*exp(-k1*t)ans=k1*(ex

3、p(-k2*t)-exp(-k1*t)/(k1-k2)3对积分式进行参数估计的源程序对积分式进行参数估计的源程序 function seqcurvefit11clear all;load seqdata;beta0=0.005 0.001;lb=0 0;ub=inf inf;beta,resnorm,residual,exitflag,output,lambda,jacobian=.lsqnonlin(seqfun,beta0,lb,ub,t,c);ci=nlparci(beta,residual,jacobian);function y=seqfun(beta,t,c)ca=exp(-bet

4、a(1)*t);cb=beta(1)/(beta(2)-beta(1)*(exp(-beta(1)*t)-exp(-beta(2)*t);y=ca-c(:,1)cb-c(:,2);4输出计算结果的源程序输出计算结果的源程序%print resultfprintf(n Estimated Parameters by Lsqnonlin():n)fprintf(t k1=%.4f%.4fn,beta(1),ci(1,2)-beta(1)fprintf(t k2=%.4f%.4fn,beta(2),ci(2,2)-beta(2)fprintf(The sum of the squares is:%.

5、1enn,sum(residual.2)%plot of fit resultstc=linspace(0,max(t),200);y_row,y_col=size(c);zeroc=zeros(200,y_col);yc=seqfun(beta,tc,zeroc);plot(t,c(:,1),ro,tc,yc(:,1),r-);hold onplot(t,c(:,2),b+,tc,yc(:,2),b-);xlabel(Time);ylabel(Concentration);hold off用函数求拟合值。用函数求拟合值。注意转置，与函数定义注意转置，与函数定义的矩阵维数一致！的矩阵维数一致！

6、resnorm5参数拟合结果参数拟合结果 seqcurvefit11Optimization terminated:relative function value changing by less than OPTIONS.TolFun.Estimated Parameters by Lsqnonlin():k1=0.0412 0.0018 k2=0.0121 0.0008 The sum of the residual squares is:2.7e-0036方法方法2：对微分式进行参数估计：对微分式进行参数估计function seqcurvefit21clear all;load seq

7、data;y_row,y_col=size(c);beta0=0.005 0.001;c0=1 0;lb=0 0;ub=inf inf;beta,resnorm,residual,exitflag,output,lambda,jacobian=.lsqnonlin(seqfun,beta0,lb,ub,t,c,y_col,c0);ci=nlparci(beta,residual,jacobian);7拟合模型和目标函数的定义拟合模型和目标函数的定义function y=seqfun(beta,t,c,y_col,c0)%Objective functiontspan=0 max(t);tt y

8、y=ode45(modeleqs,tspan,c0,beta);for col=1:y_col yc(:,col)=spline(tt,yy(:,col),t);endy=c(:,1)-yc(:,1);c(:,2)-yc(:,2);function dydt=modeleqs(t,y,beta)%Model equationdydt=-beta(1)*y(1);beta(1)*y(1)-beta(2)*y(2);beta,resnorm,residual,exitflag,output,lambda,jacobian=.lsqnonlin(seqfun,beta0,lb,ub,t,c,y_co

9、l,c0);8输出计算结果的源程序输出计算结果的源程序%print resultfprintf(n Estimated Parameters by Lsqnonlin():n)fprintf(t k1=%.4f%.4fn,beta(1),ci(1,2)-beta(1)fprintf(t k2=%.4f%.4fn,beta(2),ci(2,2)-beta(2)fprintf(The sum of the squares is:%.1enn,sum(residual.2)%plot of fit resultstspan=0 max(t);tt yc=ode45(modeleqs,tspan,c0

10、,beta);tc=linspace(0,max(t),200);yca=spline(tt,yc(:,1),tc);plot(t,c(:,1),ro,tc,yca,r-);hold onycb=spline(tt,yc(:,2),tc);plot(t,c(:,2),b+,tc,ycb,b-);xlabel(Time);ylabel(Concentration);hold off 首先解微分方程求拟合首先解微分方程求拟合值。然后用样条插值逼值。然后用样条插值逼近。近。9参数拟合结果参数拟合结果 seqcurvefit21Optimization terminated:relative func

11、tion value changing by less than OPTIONS.TolFun.Estimated Parameters by Lsqnonlin():k1=0.0412 0.0018 k2=0.0121 0.0008 The sum of the residual squares is:5.0e-003 10使用使用MATLAB进行参数估计示例进行参数估计示例q例例2：青霉素发酵过程动力学参数：青霉素发酵过程动力学参数估计：估计：在间歇发酵罐中研究青霉素发酵在间歇发酵罐中研究青霉素发酵过程动力学，微生物过程动力学，微生物Penicillium chrysogenum在一定的控

12、制条件在一定的控制条件下生长，细胞生长速率可以用逻下生长，细胞生长速率可以用逻辑模型描述：辑模型描述：青霉素的生产速率模型为：青霉素的生产速率模型为：q初始条件为：初始条件为：t=0，y1=0.29，y2=0。TimehoursCellconcentrationdry weightPenicillinconcentrationunits/ml00.180100.120220.480.0089341.460.0642461.560.2266581.730.4373701.990.6943822.621.2459942.881.43151063.432.04021183.371.92781303.

13、922.18481423.962.42041543.582.46151663.582.2831783.342.70781903.472.654211使用使用MATLAB的的dsolve求积分求积分 s=dsolve(Dy1=k1*y1*(1-y1/k1),Dy2=k3*y1-k4*y2,y1(0)=0.29,y2(0)=0)s=y1:1x1 sym y2:1x1 sym s.y1 s.y2 simplify(s.y2)ans=k1/(1+1/29*exp(-k1*t)*(100*k1-29)ans=Int(k3*k1/(1+1/29*exp(-k1*_z1)*(100*k1-29)*exp(k

14、4*_z1),_z1=0.t)*exp(-k4*t)ans=29*k3*k1*Int(exp(k4*_z1)/(29+100*exp(-k1*_z1)*k1-29*exp(-k1*_z1),_z1=0.t)*exp(-k4*t)没有确定的积分表达式没有确定的积分表达式12参数估计的源程序参数估计的源程序function PenicilliumEstclear all;load PenicilliumData%Load experimental data y_row,y_col=size(y);%Nonlinear least square estimate using lsqnonlin()b

15、eta0=0.1 4.0 0.02 0.02;y0=0.29 0.0;lb=0 0.01 0 0;ub=inf inf inf inf;beta,resnorm,residual,exitflag,output,lambda,jacobian=.lsqnonlin(Func,beta0,lb,ub,x,y,y_col,y0);ci=nlparci(beta,residual,jacobian);13拟合模型和目标函数的定义拟合模型和目标函数的定义%=function f=Func(beta,x,y,y_col,y0)%Define objective functiontspan=0 max(x

16、);tt yy=ode45(ModelEqs,tspan,y0,beta);for col=1:y_col yc(:,col)=spline(tt,yy(:,col),x);endf1=y(:,1)-yc(:,1);f2=y(:,2)-yc(:,2);f=f1;f2;%=function dydt=ModelEqs(t,y,beta)%Model equationsdydt=beta(1)*y(1)*(1-y(1)/beta(2);beta(3)*y(1)-beta(4)*y(2);14输出计算结果的源程序输出计算结果的源程序%resultfprintf(n Estimated Paramet

17、ers by Lsqnonlin():n)fprintf(t k1=%.4f%.4fn,beta(1),ci(1,2)-beta(1)fprintf(t k2=%.4f%.4fn,beta(2),ci(2,2)-beta(2)fprintf(t k3=%.4f%.4fn,beta(3),ci(3,2)-beta(3)fprintf(t k4=%.4f%.4fn,beta(4),ci(4,2)-beta(4)fprintf(The sum of the residual squares is:%.1enn,sum(residual.2)15输出计算结果的源程序输出计算结果的源程序%plotnam

18、ex=Time,hours;namey(1,1:length(Cell concentration,%dry weight)=.Cell concentration,%dry weight;namey(2,1:length(Pencillin concentration,units/ml)=.Pencillin concentration,units/ml;tspan=0 max(x);tt yy=ode45(ModelEqs,tspan,y0,beta);xc=linspace(0,max(x),200);for col=1:y_col yc(:,col)=spline(tt,yy(:,co

19、l),xc);figure(col)plot(x,y(:,col),o,xc,yc(:,col)xlabel(namex)ylabel(namey(col,:)end 16参数拟合结果参数拟合结果 penicilliumestOptimization terminated:relative function value changing by less than OPTIONS.TolFun.Estimated Parameters by Lsqnonlin():k1=0.0423 0.0040 k2=3.6850 0.1949 k3=0.0186 0.0081 k4=0.0235 0.0142 The sum of the residual squares is:1.4e+000 17参数拟合结果参数拟合结果