| 查看: 395 | 回复: 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 ] |
回复此楼» 猜你喜欢
有没有人能给点建议
已经有5人回复
-
假如你的研究生提出不合理要求
已经有12人回复
-
实验室接单子
已经有7人回复
-
全日制（定向）博士
已经有5人回复
-
萌生出自己或许不适合搞科研的想法，现在跑or等等看？
已经有4人回复
-
Materials Today Chemistry审稿周期
已经有4人回复
-
参与限项
已经有3人回复
-
对氯苯硼酸纯化
已经有3人回复
-
所感
已经有4人回复
-
要不要辞职读博？
已经有7人回复