最新多元模型回归与分析精品课件.ppt

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

1、2一、实验数据分析一、实验数据分析q由实验数据回归模型，得到模型参数前，对数据自变量间的由实验数据回归模型，得到模型参数前，对数据自变量间的线性相关性进行检验，是发现回归模型应用的可靠性和准确线性相关性进行检验，是发现回归模型应用的可靠性和准确性受限制的有效方法。性受限制的有效方法。 因自变量间的线性相关性，使得无法区分它们对因变量的作用因自变量间的线性相关性，使得无法区分它们对因变量的作用; 回归模型参数时会遇到几乎是奇异的数据矩阵，这样的模型参数有很回归模型参数时会遇到几乎是奇异的数据矩阵，这样的模型参数有很大的不确定性大的不确定性(95的参数置信度范围宽的参数置信度范围宽)。 q例：例：

2、回归二氧化硫的催化氧化速率方程：回归二氧化硫的催化氧化速率方程： 装有载铂氧化铝催化剂颗粒的微分固定床反应器中，测定二装有载铂氧化铝催化剂颗粒的微分固定床反应器中，测定二氧化硫的催化氧化速率。总压为氧化硫的催化氧化速率。总压为790 mmHg时，记录流体相的时，记录流体相的组成分压，有下表所示的速率结果，通过这些数据求取二氧组成分压，有下表所示的速率结果，通过这些数据求取二氧化硫的催化氧化速率方程。化硫的催化氧化速率方程。 9Antoine方程拟合的残差方程拟合的残差 残差虽然很小，但其分残差虽然很小，但其分布不是随机的。布不是随机的。残差的分布同正态分布相残差的分布同正态分布相比，有较大的差

3、距。比，有较大的差距。 两方面的结果充分说明了拟合回归的两方面的结果充分说明了拟合回归的AntoineAntoine方程还不能充方程还不能充分反映蒸汽压与温度间的关系，造成残差间存在关联。分反映蒸汽压与温度间的关系，造成残差间存在关联。采用采用RiedelRiedel方程拟合得到的也是类似的结果。方程拟合得到的也是类似的结果。lnPResidual Values-0.003-0.002-0.0010.0000.0010.0020.00312345678Residual fromAntoine Eq fitExpectedNormalFrequency Distribution: Residua

4、lsNo of obs02468101214161820222426-0.003-0.003-0.002-0.002-0.001-5.000e-40.0005.000e-40.0010.0020.0020.0030.0030.00410改进改进Thek-Stiel方程方程的拟合结果方程方程的拟合结果拟合误差比拟合误差比AntoineAntoine方程方程小了近一个数量级，而且小了近一个数量级，而且残差分布是随机分布的。残差分布是随机分布的。误差分布基本符合正态误差分布基本符合正态分布。分布。改进改进Thek-StielThek-Stiel方程方程描述水饱和蒸汽压的合适模型。方程方程描述水饱和蒸

5、汽压的合适模型。 lnPvPrResidual Values-7e-4-5e-4-3e-4-1e-41e-43e-45e-4-11-10-9-8-7-6-5-4Eq.(8.2.6) fit Frequency Distribution: ResidualsNo of obs01020304050607080-7e-4-6e-4-5e-4-4e-4-3e-4-2e-4-1e-401e-42e-43e-44e-4Eq.(6) fitExpectedNormal11二、回归模型的选择（二、回归模型的选择（2）q前面说明了前面说明了模型参数较少模型参数较少时会出现拟合残差时会出现拟合残差的分布不是随机

6、的，而是呈现某种分布，相的分布不是随机的，而是呈现某种分布，相互关联。互关联。q在模型回归拟合数据的过程中，如在模型回归拟合数据的过程中，如模型参数模型参数过多过多会出现什么情况？如何判断回归拟合模会出现什么情况？如何判断回归拟合模型中有过多的参数呢？型中有过多的参数呢？12丙烷在氢型丝光沸石上的吸附平衡丙烷在氢型丝光沸石上的吸附平衡q例：例：选用不同吸附方程拟合丙烷在氢型丝光沸石体系选用不同吸附方程拟合丙烷在氢型丝光沸石体系303K的吸附平衡数的吸附平衡数据。据。q目标：目标：说明如何对模型拟合结果进行统计分析，确定模型拟合的好坏、模说明如何对模型拟合结果进行统计分析，确定模型拟合的好坏、模

7、型参数的可靠性和准确性，从而进行拟合模型的选择。型参数的可靠性和准确性，从而进行拟合模型的选择。 P, kPaq, mmol/gP, kPaq, mmol/gP, kPaq, mmol/gP, kPaq, mmol/g0.100.09 1.080.4812.670.81115.891.140.140.12 1.470.5116.700.85140.071.170.220.18 1.510.5324.810.90158.901.190.330.24 2.270.5934.280.95176.761.200.410.30 3.220.6443.850.98193.371.220.490.31 4.

8、720.6954.621.02206.811.240.570.36 5.060.7065.791.040.770.41 7.390.7573.191.060.990.4410.260.7994.661.09表表8-4 303K8-4 303K时丙烷在氢型丝光沸石上的吸附平衡数据时丙烷在氢型丝光沸石上的吸附平衡数据 13具有代表性的、也是适用性较广的模型具有代表性的、也是适用性较广的模型 q1、Lanmuir (L)双参数方程：双参数方程： q2、Freundlich (F)双参数方程：双参数方程： nn aPaPm1(8.2.8) q3、BET双参数方程：双参数方程：q4、Langmuir-F

9、reundlich (LF)三参数方程：三参数方程：q5、三参数方程：、三参数方程： q6、Toth三参数方程：三参数方程： q7、扩展的、扩展的LF方程（五参数）：方程（五参数）： q8、(14)式的特殊形式（四参数）式的特殊形式（四参数） ： naPbnn cxxxcxm()()11xP P/0nn aPaPmbb1Pbc111exp() n nm/nn PbPmcc()/1nn aPdcPaPcPmbebe1nn aPbc bd PaPcPmbdbd()11(8.2.9) (8.2.10) (8.2.11) (8.2.12) (8.2.13) (8.2.14) (8.2.15) 14各模

10、型的计算结果各模型的计算结果 Eq.(8)Eq.(9)Eq.(10)Eq.(11)Eq.(12)Eq.(13)Eq.(14)Eq.(15)nm1.0840.051/0.9760.0250.4380.0534.62317.221.5350.2570.7580.0880.7690.068a0.5530.1310.4460.034/1.3820.107/1.3290.4961.4270.144b/0.200.018/0.4940.0620.6781.6580.5490.0760.9420.1070.9750.058c/812.1116.3/20.9652.650.3410.0591.9071.593

11、0.0170.003d/0.0240.0480.9120.046e/1.6380.586/s29.09010-27.72110-25.13410-23.79910-22.9933.15410-21.04810-28.68610-3R20.987900.991270.996140.997950.998580.998590.999860.99990表表8 85 5 吸附等温线关联的参数值、方差和回归系数吸附等温线关联的参数值、方差和回归系数 从表中可看出，方程从表中可看出，方程(8(814)14)拟合方差逐渐减少，回归系数更拟合方差逐渐减少，回归系数更接近接近1 1（方程（方程(12)(12)是通

12、过压力数据来拟合的，故拟合方差和其是通过压力数据来拟合的，故拟合方差和其它方程的结果不是在同一数量级上它方程的结果不是在同一数量级上) )。由方程。由方程(14)(14)的五参数形的五参数形式改进的方程式改进的方程(15)(15)式获得的结果最好，实验数据点几乎完全落式获得的结果最好，实验数据点几乎完全落在方程（在方程（1515）式的曲线上）式的曲线上( (见下图见下图) )。15方程方程(13)和方程和方程(15)的拟合结果的拟合结果 0.11.010.0100.0p, kPa0.00.51.0n, mmol/gExp. dataToth Eq.(13)Eq.(15)方程方程(15)(15)

13、式获得的结果最好，实验数据点几乎完式获得的结果最好，实验数据点几乎完全落在方程（全落在方程（1515）式的曲线上。）式的曲线上。16判断模型参数是否过少的依据判断模型参数是否过少的依据q通过对方程通过对方程(13)和五参数方程和五参数方程(15)的残差进行分析，的残差进行分析， 方程方程(13)因参数过少，吸附量的计算误差与实验吸附量之间存在着某种分布。因参数过少，吸附量的计算误差与实验吸附量之间存在着某种分布。 方程方程(15) 计算误差在零的两边是随机分布的，看不出规律性。因此，拟合计算计算误差在零的两边是随机分布的，看不出规律性。因此，拟合计算误差有无规律性的分布是判断模型参数是否过少的

14、依据。误差有无规律性的分布是判断模型参数是否过少的依据。因此，拟合计算误差有无规律性的分布是判断模型参数是因此，拟合计算误差有无规律性的分布是判断模型参数是否过少的依据。否过少的依据。 Predicted versus Residual Valuesn, mmol/gResidual Values-0.08-0.06-0.04-0.020.000.020.040.060.00.20.40.60.81.01.21.4Predicted versus Residual Valuesn, mmol/gResidual Values-0.020-0.015-0.010-0.0050.0000.0050

15、.0100.0150.0200.00.20.40.60.81.01.21.4方程方程(13)(13)的拟合误差的拟合误差 方程方程(15)(15)的拟合误差的拟合误差 17判断模型参数是否过多的依据判断模型参数是否过多的依据Eq.(14)Eq.(15)nm0.7580.0880.7690.068a1.3290.4961.4270.144b0.9420.1070.9750.058c1.9071.5930.0170.003d0.0240.0480.9120.046e1.6380.586/s21.04810-28.68610-3R20.999860.99990q在方程在方程(14)的计算结果中，有的

16、计算结果中，有些参数些参数95%的置信度较大，说的置信度较大，说明这些参数之间有联系，不是明这些参数之间有联系，不是独立的。独立的。q而对于方程而对于方程(14)的五参数形式，的五参数形式，即方程即方程(15)，其所有参数的，其所有参数的95%置信度都较小。事实上，置信度都较小。事实上，方程方程(15)就是据此分析对吸附就是据此分析对吸附平衡理论作进一步研究而获得平衡理论作进一步研究而获得的。的。因此，拟合参数因此，拟合参数95%95%的置信度是否较大是判断模型参数是否的置信度是否较大是判断模型参数是否过多的依据。过多的依据。18回归模型的选择总结回归模型的选择总结q模型参数较少模型参数较少时

17、会出现拟合残差的分布不是时会出现拟合残差的分布不是随机的，而是呈现某种分布，相互关联。残随机的，而是呈现某种分布，相互关联。残差的分布偏离正态分布较远。差的分布偏离正态分布较远。 q模型参数过多模型参数过多会出现某些参数会出现某些参数95%的置信度的置信度较大，说明这些参数之间有联系，不是独立较大，说明这些参数之间有联系，不是独立的。的。19习题习题q研究二硫化碳饱和蒸汽压的模型回归问题。（研究二硫化碳饱和蒸汽压的模型回归问题。（P266，Ex8.3）二硫化碳的基本性质：二硫化碳的基本性质： 临界温度为临界温度为273.05273.05 临界压力为临界压力为72.87 72.87 atmatm

18、。 温度，蒸汽压，mmHg温度，蒸汽压，mmHg-701.610198.0-603.520297.5-507.130432.7-4014.040616.7-3026.250995.6-2046.5601170.4-1078.8701558.00127.320Statistica的非线性估计的非线性估计q非线性估计方法非线性估计方法 User-Specified Regression, least square 可以计算95%置信区间21“least square”与与“Custom Loss ”比较比较q“least square”计算结果（采用计算结果（采用Levenberg-Marquar

19、dt方法）方法）Model is: q=(nm*a*p*b*(b+c*(b+d)*p*d)/(1+a*p*b*(1+c*p*d) (propane-mor.sta)Dep. Var. : qLevel of confidence: 95.0% ( alpha=0.050)EstimateStandarderrort-valuedf = 28p-levelLo. ConfLimitUp. ConfLimitnmabcd0.7637600.032019 23.85354 0.000000 0.698173 0.8293471.4359130.069004 20.80900 0.000000 1.2

20、94564 1.5772620.9822270.028027 35.04555 0.000000 0.924816 1.0396380.0163330.0016779.73910 0.000000 0.012898 0.0197680.9168660.022152 41.38886 0.000000 0.871489 0.962243Model: q=(nm*a*p*b*(b+c*(b+d)*p*d)/(1+a*p*b*(1+c*p*d). (propane-mor.sta)Dep. var: q Loss: (OBS-PRED)*2Final loss: .002032166 R=.9997

21、5 Variance explained: 99.950%N=33nmabcdEstimate0.763768 1.435910 0.982209 0.016336 0.916840q“Custom Loss ”计算结果（采用计算结果（采用Quasi-Newton法）法）Matrix ill conditioned; cannot compute standard errors.EConf22“Custom Loss ”的方差分析与迭代步骤的方差分析与迭代步骤 “Custom “Custom Loss ”Loss ”无迭无迭代历史纪录，代历史纪录，协方差分析结协方差分析结果已出现病态。果已出现

22、病态。Model is: q=(nm*a*p*b*(b+c*(b+d)*p*d)/(1+a*p*b*(1+c*p*d) (propane-mor.sta)Dep. Var. : qEffect1Sum of Sqares2DF3Mean Squares4F-value5p-valueRegressionResidualTotalCorrected TotalRegression vs.Corrected Total21.505475.000004.30109459262.190.0000000.0020328.000000.00007321.5075033.000004.0675032.000

23、0021.505475.000004.30109433.840.000000Model is: q=(nm*a*p*b*(b+c*(b+d)*p*d)/(1+a*p*b*(1+c*p*d) (propane-mor.sta)Dep. Var. : qLossFunctionnmabcd123456789104.629683 0.100000 0.100000 0.100000 0.100000 0.1000004.294027 0.184348 0.196978 0.173970 0.528173 0.5570214.227658 0.191476 0.221943 0.188990 0.58

24、7318 0.5716294.102290 0.204036 0.266635 0.217677 0.697367 0.5994433.874137 0.224709 0.339318 0.270956 0.885848 0.6511153.429648 0.259087 0.453795 0.373779 1.198671 0.7516122.556049 0.312350 0.616018 0.561417 1.656081 0.9378231.036472 0.392841 0.756563 0.881424 1.932017 1.2618210.899994 0.404754 0.76

25、4595 0.929852 1.734310 1.3103860.755383 0.425557 0.992230 0.962118 0.684109 1.321683Covariance matrix cannot be computed. 23Statistica非线性估计的残差分析非线性估计的残差分析残差分布情况残差分布情况残差对预测值作图残差对预测值作图24对方程对方程(15)的残差分析的残差分析Predicted versus Residual Values123456789101112131415161718192021222324252627282930313233-0.20.0

26、0.20.40.60.81.01.21.4Predicted Values-0.020-0.015-0.010-0.0050.0000.0050.0100.0150.020Residual ValuesFrequency Distribution:-0.020-0.015-0.010-0.0050.0000.0050.0100.0150.020Residual Values01234567891011No of obsHistogram of Histogram of residualsresidualsResidual vs. Residual vs. PredictedPredicted2

27、5误差正态分布图误差正态分布图Frequency Distribution: Residuals Expected Normal-0.0030-0.0025-0.0020-0.0015-0.0010-0.00050.00000.00050.00100.00150.00200.00250.00300.003502468101214161820222426No of obsFrequency Distribution:-0.0030-0.0025-0.0020-0.0015-0.0010-0.00050.00000.00050.00100.00150.00200.00250.00300.00350

28、.0040Residual Values0510152025303540No of obs“least square”计算的误差正态分布图计算的误差正态分布图“Custom Loss ”计算的误差正态分布图计算的误差正态分布图qAntoine 方程拟合结果方程拟合结果26非线性函数的管理非线性函数的管理27三、三、MATLAB的拟合函数的拟合函数q多项式拟合函数多项式拟合函数polyfitq非线性最小二乘法非线性最小二乘法 lsqnonlin()非线性最小二乘（优化问题）非线性最小二乘（优化问题） lsqcurvefit()非线性最小二乘曲线拟合非线性最小二乘曲线拟合 nlinfit()前两种

29、的简化版本前两种的简化版本 nlparci()计算参数的置信区间计算参数的置信区间 nlpredci()计算预测值的置信区间计算预测值的置信区间 nlintool()nlinfit()、nlparci()、nlpredci()的集成图形用的集成图形用户界面拟合函数户界面拟合函数q注意：不同的拟合函数命令，其优化目标函数定义注意：不同的拟合函数命令，其优化目标函数定义以及调用形式不同，注意区分！以及调用形式不同，注意区分！28（一）（一）Polyfit的使用的使用qp = polyfit(x0, y0, n) 其中其中x0和和y0分别为观察节点和观察值向量；分别为观察节点和观察值向量；n表示插值

30、多项式的次数；表示插值多项式的次数；输出值输出值p表示插值多项式的系数。表示插值多项式的系数。q例：某实验中测得一组数据，其值如下：例：某实验中测得一组数据，其值如下： x l 2 3 4 5 y 1.3 1.8 2.2 2.9 3.5 已知已知x和和y成线性关系，即成线性关系，即y=kx+b，求系数，求系数k和和b x=1 2 3 4 5;y=1.3 1.8 2.2 2.9 3.5;p=polyfit(x,y,1)y1=polyval(p,x);plot(x,y1)hold on;plot(x,y,b*);p = 0.55 0.6911.522.533.544.5511.522.533.5也

31、就是，也就是，k为为0.55，b为为0.6929（二）（二）lsqcurvefit的使用的使用q方程（目标函数）方程（目标函数） Find coefficients beta that best fit the equationq调用形式调用形式 beta= lsqcurvefit(fun, beta0,xdata,ydata) beta= lsqcurvefit(fun, beta0,xdata,ydata,lb,ub) beta = lsqcurvefit(fun, beta0,xdata,ydata,lb,ub,options) xdata和和ydata分别为观察节点和观察值向量；分别为

32、观察节点和观察值向量；fun自定义的非线性拟合模型；自定义的非线性拟合模型；beta0拟合参数的初始值；拟合参数的初始值；beta拟合模型中的参数；拟合模型中的参数；lb,ub拟合参数初值的边界值拟合参数初值的边界值, lb = beta 1 % two output arguments J = . % Jacobian of the function evaluated at betaend options = optimset(Jacobian,on)qFor examplefunction F = myfun(beta,xdata)F = beta(1)*exp(beta(2)*xdat

33、a);33lsqcurvefit应用示例应用示例% Assume you determined xdata and ydata experimentallyxdata = 0.9 1.5 13.8 19.8 24.1 28.2 35.2 60.3 74.6 81.3;ydata = 455.2 428.6 124.1 67.3 43.2 28.1 13.1 -0.4 -1.3 -1.5;beta0 = 100; -1 % Starting guessbeta,resnorm = lsqcurvefit(myfun, beta0,xdata,ydata)function F = myfun(be

34、ta,xdata)F = beta (1)*exp(beta(2)*xdata);beta = 498.8309 -0.1013 resnorm = 9.504921b xybe34（三）（三）lsqnonlin的使用的使用q方程（目标函数）方程（目标函数） Find coefficients beta that best fit the equationq调用形式调用形式 beta = lsqnonlin(fun, beta 0) beta = lsqnonlin(fun, beta 0,lb,ub) beta = lsqnonlin(fun, beta 0,lb,ub,options) f

35、un自定义的优化函数；自定义的优化函数；beta 0优化参数的初始值；优化参数的初始值；beta拟合模型中的参数；拟合模型中的参数；lb,ub拟合参数初值的边界值拟合参数初值的边界值, lb = beta 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*

36、(exp(-k2*t)-exp(-k1*t)/(k1-k2)tkAeC1)(21121tktkBeekkkC42对积分式进行参数估计的源程序对积分式进行参数估计的源程序 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); fu

37、nction y = seqfun(beta,t,c)ca=exp(-beta(1)*t);cb=beta(1)/(beta(2)-beta(1)*(exp(-beta(1)*t)-exp(-beta(2)*t);y=ca-c(:,1) cb-c(:,2); 43输出计算结果的源程序输出计算结果的源程序 % 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

38、(2,2)-beta(2)fprintf( The sum of the squares is: %.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 o

39、ff用函数求拟合值。用函数求拟合值。注意转置，与函数定义注意转置，与函数定义的矩阵维数一致！的矩阵维数一致！44参数拟合结果参数拟合结果 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-003020406080100

40、12014016018020000.10.20.30.40.50.60.70.80.91TimeConcentration45方法方法2：对微分式进行参数估计：对微分式进行参数估计function seqcurvefit21clear all;load seqdata;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);

41、ci = nlparci(beta,residual,jacobian); 46拟合模型和目标函数的定义拟合模型和目标函数的定义function y = seqfun(beta,t,c,y_col,c0) % Objective functiontspan = 0 max(t);tt yy = 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,be

42、ta) % 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_col,c0);47输出计算结果的源程序输出计算结果的源程序 % print resultfprintf(n Estimated Parameters by Lsqnonlin():n)fprintf(t k1 = %.4f %.4fn,beta(1),ci(1,2)-bet

43、a(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,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),

44、tc);plot(t,c(:,2),b+,tc,ycb,b-);xlabel(Time);ylabel(Concentration);hold off 首先解微分方程求拟合首先解微分方程求拟合值。然后用样条插值逼值。然后用样条插值逼近。近。48参数拟合结果参数拟合结果 seqcurvefit21Optimization terminated: relative function value changing by less than OPTIONS.TolFun. Estimated Parameters by Lsqnonlin(): k1 = 0.0412 0.0018 k2 = 0.01

45、21 0.0008 The sum of the residual squares is: 5.0e-003 02040608010012014016018020000.10.20.30.40.50.60.70.80.91TimeConcentration49使用使用MATLAB进行参数估计示例进行参数估计示例q例例2：青霉素发酵过程动力学参数：青霉素发酵过程动力学参数估计：估计： 在间歇发酵罐中研究青霉素发酵在间歇发酵罐中研究青霉素发酵过程动力学，微生物过程动力学，微生物Penicillium chrysogenum在一定的控制条件在一定的控制条件下生长，细胞生长速率可以用逻下生长，细胞生长

46、速率可以用逻辑模型描述：辑模型描述： 青霉素的生产速率模型为：青霉素的生产速率模型为： q初始条件为：初始条件为：t=0， y1=0.29，y2=0 。211111kyykdtdy24132ykykdtdyTimehoursCellconcentrationdry weightPenicillinconcentrationunits/ml00.180100.120220.480.0089341.460.0642461.560.2266581.730.4373701.990.6943822.621.2459942.881.43151063.432.04021183.371.92781303.92

47、2.18481423.962.42041543.582.46151663.582.2831783.342.70781903.472.654250使用使用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

48、-29)*exp(k4*_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) 没有确定的积分表达式没有确定的积分表达式51参数估计的源程序参数估计的源程序function PenicilliumEstclear all;load PenicilliumData % Load experimental data y_row,y_col=size(y); % Nonlinear least square e

49、stimate using lsqnonlin()beta0 = 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); 52拟合模型和目标函数的定义拟合模型和目标函数的定义% =function f = Func(beta,x,y,y_co

50、l,y0) % Define objective functiontspan = 0 max(x);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)