| 查看: 424 | 回复: 1 | |||
| 当前主题已经存档。 | |||
allenhero1228金虫 (小有名气)
|
[交流]
求助 matlab做模型拟合时的问题
|
||
|
我用matlab求四个动力学模型的参数，并输出理论和实验值及二者的相关系数，结果程序运行的结果是后三个模型算出的理论值和相关系数都是一模一样，而且这四个模型的理论值与实验值相差很大，不知道这是什么原因，望各位高手帮忙给看一下，谢谢啦 function KineticsEst % 动力学参数估计及模型辨识—判别四个模型中哪个最好 clear all clc global rA Pco Pco2 Ph2 Pm Ph2o Kf %553K Pco = [1.3487 1.3863 1.3253 1.3140 1.3072 1.3052 1.3674]; Ph2 = [1.6901 1.7165 1.6863 1.6728 1.6699 1.6587 1.7991]; Pco2 = [0.0599 0.0576 0.0637 0.0610 0.0581 0.0619 0.0675]; Pm =[0.0032 0.0032 0.0028 0.0029 0.0030 0.0029 0.0030]; Ph2o = [0.0463 0.0471 0.0491 0.0510 0.0516 0.0524 0.0525]; rA = [0.3118 0.3177 0.3301 0.3465 0.3510 0.3554 0.3560]; Kf = 0.0006; %---- IDmodel = 'a'; x0 = [1 1 1 1 1]; option1=optimset('MaxFunEvals',10000,'Maxiter',10000,'tolfun',1e-10); [x,resnorm,residual,exitflag,output,lambda,jacobian] = ... lsqnonlin(@ObjFunc,x0,zeros(size(x0)),inf*ones(size(x0)),option1,IDmodel); rAc=x(1)*Pco.*Ph2.*(1-Pm./(Kf*Pco.*Ph2.^2))./((1+x(2)*Pco+x(3)*Pco2).*(1+x(4)*Ph2+x(5)*Ph2o)) B=[rA;rAc] %rA为实验值，rAc为理论值 R=corrcoef(rA',rAc') %求模型的相关系数 exitflag %查看计算终值的条件 fprintf('\n Model (a):\n') fprintf(' k = %.5d Kco = %.3d Kco2 = %.3d Kh2 = %.3d Kh2o = %.3d\n',x(1),x(2),x(3),x(4),x(5)) fprintf(' The sum of the squares is: %.3f\n\n',resnorm) IDmodel = 'b'; x0 = [1 1 1 1 1 1]; option1=optimset('MaxFunEvals',30000,'Maxiter',10000,'tolfun',1e-6); [x,resnorm,residual,exitflag,output,lambda,jacobian] = ... lsqnonlin(@ObjFunc,x0,zeros(size(x0)),inf*ones(size(x0)),option1,IDmodel); rAc=x(1)*Pco.*Ph2.^2.*(1-Pm./(Kf*Pco.*Ph2.^2))./((1+x(2)*Pco+x(3)*Pco2+x(4)*Pm).*(1+x(5)*Ph2+x(6)*Ph2o)) B=[rA;rAc] R=corrcoef(rA',rAc') exitflag fprintf('\n Model (b):\n') fprintf(' k = %.3d Kco = %.3d Kco2 = %.3d Km = %.3d Kh2 = %.3d Kh2o = %.d\n',x(1),x(2),x(3),x(4),x(5),x(6)) fprintf(' The sum of the squares is: %.3f\n\n',resnorm) IDmodel = 'c'; x0 = [1 1 1 1 1]; option1=optimset('MaxFunEvals',10000,'Maxiter',10000,'tolfun',1e-10); [x,resnorm,residual,exitflag,output,lambda,jacobian] = ... lsqnonlin(@ObjFunc,x0,zeros(size(x0)),inf*ones(size(x0)),option1,IDmodel); B=[rA;rAc] R=corrcoef(rA',rAc') exitflag fprintf('\n Model (c):\n') fprintf(' k = %.3d Kco = %.3d Kco2 = %.3d Kh2 = %.3d Kh2o = %.3d\n',x(1),x(2),x(3),x(4),x(5)) fprintf(' The sum of the squares is: %.3f\n\n',resnorm) IDmodel = 'd'; x0 = [1 1 1 1]; option1=optimset('MaxFunEvals',15000,'Maxiter',10000,'tolfun',1e-15); [x,resnorm,residual,exitflag,output,lambda,jacobian] = ... lsqnonlin(@ObjFunc,x0,zeros(size(x0)),inf*ones(size(x0)),option1,IDmodel); B=[rA;rAc] R=corrcoef(rA',rAc') exitflag fprintf('\n Model (d):\n') fprintf(' k = %.3d Kco = %.3d Kco2 = %.3d Kh2=%.3d\n',x(1),x(2),x(3),x(4)) fprintf(' The sum of the squares is: %.3f\n\n',resnorm) % ------------------------------------------------------------------ function f = ObjFunc(x,IDmodel) global rA Pco Pco2 Ph2 Pm Ph2o Kf switch IDmodel case 'a' % Model (a), x(1)=k; x(2)=Kco; x(3)=Kco2; x(4)=Kh2; x(5)=Kh2o; f = rA - x(1)*Pco.*Ph2.*(1-Pm./(Kf*Pco.*Ph2.^2))./((1+x(2)*Pco+x(3)*Pco2).*(1+x(4)*Ph2+x(5)*Ph2o)); case 'b' % Model (b), x(1)=k; x(2)=Kco; x(3)=Kco2; x(4)=Km; x(5)=Kh2; x(6)=Kh2o; f = rA - x(1)*Pco.*Ph2.^2.*(1-Pm./(Kf*Pco.*Ph2.^2))./((1+x(2)*Pco+x(3)*Pco2+x(4)*Pm).*(1+x(5)*Ph2+x(6)*Ph2o)); case 'c' % Model (c), x(1)=k; x(2)=Kco; x(3)=Kco2;x(5)=Kh2; x(6)=Kh2o; f = rA - x(1)*Pco.*Ph2.*(1-Pm./(Kf*Pco.*Ph2.^2))./((1+x(2)*Pco+x(3)*Pco2).*(1+x(4).^0.5*Ph2.^0.5+x(5)*Ph2o)); case 'd' % Model (d), x(1)=k; x(2)=Kco; x(3)=Kco2; x(4)=Kh2; f = rA - x(1)*Pco.*Ph2.^2.*(1-Pm./(Kf*Pco.*Ph2.^2))./(1+x(2)*Pco+x(3)*Pco2+x(4)*Ph2).^3; otherwise disp('Unknown model.') end 计算结果如下： Optimization terminated: first-order optimality less than OPTIONS.TolFun, and no negative/zero curvature detected in trust region model. rAc = -1.1888e-012 -9.5004e-013 -7.3257e-013 -9.6625e-013 -1.143e-012 -1.06e-012 -4.0011e-013 B = 0.3118 0.3177 0.3301 0.3465 0.351 0.3554 0.356 -1.1888e-012 -9.5004e-013 -7.3257e-013 -9.6625e-013 -1.143e-012 -1.06e-012 -4.0011e-013 R = 1 0.3244 0.3244 1 exitflag = 1 Model (a): k = 5.19864e-010 Kco = 1.187e+001 Kco2 = 1.181e+001 Kh2 = 1.189e+001 Kh2o = 1.174e+001 The sum of the squares is: 0.803 Optimization terminated: first-order optimality less than OPTIONS.TolFun, and no negative/zero curvature detected in trust region model. B = 0.3118 0.3177 0.3301 0.3465 0.351 0.3554 0.356 -6.4994e-011 -5.2725e-011 -3.998e-011 -5.232e-011 -6.1787e-011 -5.6925e-011 -2.3264e-011 R = 1 0.35169 0.35169 1 exitflag = 1 Model (b): k = 3.277e-008 Kco = 1.680e+001 Kco2 = 1.676e+001 Km = 1.654e+001 Kh2 = 1.682e+001 Kh2o = 2e+001 The sum of the squares is: 0.803 Optimization terminated: first-order optimality less than OPTIONS.TolFun, and no negative/zero curvature detected in trust region model. B = 0.3118 0.3177 0.3301 0.3465 0.351 0.3554 0.356 -6.4994e-011 -5.2725e-011 -3.998e-011 -5.232e-011 -6.1787e-011 -5.6925e-011 -2.3264e-011 R = 1 0.35169 0.35169 1 exitflag = 1 Model (c): k = 8.178e-009 Kco = 2.015e+001 Kco2 = 2.017e+001 Kh2 = 2.007e+001 Kh2o = 2.011e+001 The sum of the squares is: 0.803 Optimization terminated: first-order optimality less than OPTIONS.TolFun, and no negative/zero curvature detected in trust region model. B = 0.3118 0.3177 0.3301 0.3465 0.351 0.3554 0.356 -6.4994e-011 -5.2725e-011 -3.998e-011 -5.232e-011 -6.1787e-011 -5.6925e-011 -2.3264e-011 R = 1 0.35169 0.35169 1 exitflag = 1 Model (d): k = 2.220e-014 Kco = 9.605e+000 Kco2 = 9.539e+000 Kh2=9.656e+000 The sum of the squares is: 0.803 [ Last edited by allenhero1228 on 2008-4-15 at 11:09 ] |
回复此楼» 猜你喜欢
295求调剂
已经有3人回复
-
一志愿C9的化学工程（085602） 340分，感觉校内调剂无望，求调剂
已经有9人回复
-
070300化学求调剂
已经有16人回复
-
363求调剂
已经有4人回复
-
283分求调剂
已经有3人回复
-
材料化工306分找合适调剂
已经有7人回复
-
一志愿北京化工085600 310分求调剂
已经有4人回复
-
292求调剂
已经有9人回复
-
调剂
已经有4人回复
-
一志愿9材料学硕297已过六级求调剂推荐
已经有5人回复
