| 查看: 659 | 回复: 1 | |||
[求助]
关于动力学参数拟合的matlab程序的一些问题（懂动力学拟合的大神进）
|
| 各位大神，小弟正在编写matlab程序拟合动力学参数，在程序调试的过程中发现，当确定反应级数后，k0的初值无论怎么改变matlab算得的结果不变，这是怎么回事，不是说matlab程序是以最小二乘法计算的，属于局部最优，对于初值的选取很重要，因此一般初值变化结果会跟着变化的呀，可我在调试程序时，不论怎么改变k0的初值都不变，这是怎么回事呢，求大神告知，小弟感激不尽，金币送上！！！ |
回复此楼» 猜你喜欢
职称评审没过，求安慰
已经有51人回复
-
毕业后当辅导员了，天天各种学生超烦
已经有5人回复
-
26申博自荐
已经有3人回复
-
A期刊撤稿
已经有4人回复
-
垃圾破二本职称评审标准
已经有17人回复
-
投稿Elsevier的Neoplasia杂志，到最后选publishing options时页面空白，不能完成投稿
已经有22人回复
-
EST投稿状态问题
已经有7人回复
|
clear all clc format short global y_exp y_sim lspan = [0,0.0466]; % t=w/F(MEOH.0),W:催化剂的质量，F(MEOH.0)：甲醇的初始进料流率 k0 = [2000 100 20 1]; lb = [0 0 0 0]; ub = [+inf +inf +inf +inf ]; y0 = [0 0 0 0 ]; y_exp =[0.0288 0.0014 0.0005 0.0002; 0.0264 0.0015 0.0004 0.0004; 0.0159 0.0006 0.0002 0.0003]; %原数据 0.0109（更改） 0.0006 0.0002 0.0003 % 使用函数lsqnonlin()进行参数估计 [k,resnorm,residual,exitflag,output,lambda,jacobian] = lsqnonlin(@ObjFunc,k0,lb,ub,optimset('TolFun',1.0000e-20),lspan,y0,y_exp); ci = nlparci(k,residual,jacobian) %[k,resnorm,residual,exitflag,output,lambda,jacobian] = lsqnonlin(@ObjFunc4LNL,k0,lb,ub,optimset('TolFun',1.0000e-6),y0,y_exp); %ci = nlparci(k,residual,jacobian) fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n') fprintf('\tk1 = %.4f ± %.4f\n',k(1),ci(1,2)-k(1)) fprintf('\tk2 = %.4f ± %.4f\n',k(2),ci(2,2)-k(2)) fprintf('\tk3 = %.4f ± %.4f\n',k(3),ci(3,2)-k(3)) fprintf('\tk4 = %.4f ± %.4f\n',k(4),ci(4,2)-k(4)) fprintf('\tThe sum of the squares is: %.1e\n\n',resnorm) function f = ObjFunc(k,lspan,y0,y_exp) % 目标函数 [t,y_sim2] = ode45(@KineticsEqs2,lspan,y0,[],k) f5 = 1*(y_sim2(end,1)-y_exp(1,1)); f6= 1*(y_sim2(end,2)-y_exp(1,2)); f7= 1*(y_sim2(end,3)-y_exp(1,3)); f8= 1*(y_sim2(end,4)-y_exp(1,4)); [t,y_sim4] = ode45(@KineticsEqs4,lspan,y0,[],k) f13 = 1*(y_sim4(end,1)-y_exp(2,1)); f14 =1*(y_sim4(end,2)-y_exp(2,2)); f15 =1*(y_sim4(end,3)-y_exp(2,3)); f16 =1*(y_sim4(end,4)-y_exp(2,4)); [t,y_sim5] = ode45(@KineticsEqs5,lspan,y0,[],k) f17 = 1*(y_sim5(end,1)-y_exp(3,1)); f18 =1*(y_sim5(end,2)-y_exp(3,2)); f19 =1*(y_sim5(end,3)-y_exp(3,3)); f20 =1*(y_sim5(end,4)-y_exp(3,4)); f=[f5 f6 f7 f8; f13 f14 f15 f16; f17 f18 f19 f20] % ------------------------------------------------------------------ function dydl = KineticsEqs2(l,y,k) %P=[0.6 0.8 0.5 0.7 0.9];P为总压 %nT0=[0.154284 0.161404 0.161404 0.166152 0.332303 ]; % n:进料总摩尔数 %yO20=[0.0769 0.0882 0.1176 0.1429 0.0857 ]; %进料O2组成 %yMEOH0=[0.3077 0.2941 0.2941 0.2857 0.1429 ]; %进料甲醇（MEOH）组成 %A*Rho=0.05832 yO2=0.0882*(1+0.5*y(1)-0.5*y(3))-(0.5*y(1)+y(2)+0.5*y(3)); yMEOH=0.2941*(1+0.5*y(1)-0.5*y(3))-(3*y(1)+2*y(2)+y(3)+2*y(4)); % yO2=yO20*(1+0.5*y(1)-0.5*y(3))-(0.5*y(1)+y(2)+0.5*y(3)); % yMEOH=yMEOH0*(1+0.5*y(1)-0.5*y(3))-(3*y(1)+2*y(2)+y(3)+2*y(4)); dydl=... [0.05832*(2+y(1)-y(3))/4*0.161404*((2+y(1))*k(1)*0.8^2.5*yMEOH^2*yO2^0.5-y(1)*k(3)*0.8^2*yMEOH*yO2); 0.05832*(2+y(1)-y(3))/4*0.161404*(2*k(2)*0.8^2*yMEOH^1*yO2^1+y(2)*k(1)*0.8^2.5*yMEOH^2*yO2^0.5-y(2)*k(3)*0.8^2*yMEOH*yO2); 0.05832*(2+y(1)-y(3))/4*0.161404*((2-y(3))*k(3)*0.8^2*yMEOH*yO2+y(3)*k(1)*0.8^2.5*yMEOH^2*yO2^0.5); 0.05832*(2+y(1)-y(3))/4*0.161404*(2*k(4)*0.8^1*yMEOH+y(4)*k(1)*0.8^2.5*yMEOH^2*yO2^0.5-y(4)*k(3)*0.8^2*yMEOH*yO2); ]; % dydl=... % [A*Rho(2+y(1)-y(3))/4*nT0*((2+y(1))*k(1)*yMEOH^1*yO2^1-y(1)*k(3)*yMEOH*yO2); % A*Rho(2+y(1)-y(3))/4*nT0*(2*k(2)*yMEOH*yO2+y(2)*k(1)*yMEOH^1*yO2^1-y(2)*k(3)*yMEOH*yO2); % A*Rho(2+y(1)-y(3))/4*nT0*((2-y(3))*k(3)*yMEOH*yO2+y(3)*k(1)*yMEOH^1*yO2^1); % A*Rho(2+y(1)-y(3))/4*nT0*(2*k(4)*yMEOH+y(4)*k(1)*yMEOH^1*yO2^1-y(4)*k(3)*yMEOH*yO2); function dydl = KineticsEqs4(l,y,k) %P=[0.6 0.8 0.5 0.7 0.9];P为总压 %nT0=[0.154284 0.161404 0.161404 0.166152 0.332303 ]; % n:进料总摩尔数 %yO20=[0.0769 0.0882 0.1176 0.1429 0.0857 ]; %进料O2组成 %yMEOH0=[0.3077 0.2941 0.2941 0.2857 0.1429 ]; %进料甲醇（MEOH）组成 %A*Rho=0.05832 yO2=0.1429*(1+0.5*y(1)-0.5*y(3))-(0.5*y(1)+y(2)+0.5*y(3)); yMEOH=0.2857*(1+0.5*y(1)-0.5*y(3))-(3*y(1)+2*y(2)+y(3)+2*y(4)); % yO2=yO20*(1+0.5*y(1)-0.5*y(3))-(0.5*y(1)+y(2)+0.5*y(3)); % yMEOH=yMEOH0*(1+0.5*y(1)-0.5*y(3))-(3*y(1)+2*y(2)+y(3)+2*y(4)); dydl=... [0.05832*(2+y(1)-y(3))/4*0.166152*((2+y(1))*k(1)*0.7^2.5*yMEOH^2*yO2^0.5-y(1)*k(3)*0.7^2*yMEOH*yO2); 0.05832*(2+y(1)-y(3))/4*0.166152*(2*k(2)*0.7^2*yMEOH^1*yO2^1+y(2)*k(1)*0.7^2.5*yMEOH^2*yO2^0.5-y(2)*k(3)*0.7^2*yMEOH*yO2); 0.05832*(2+y(1)-y(3))/4*0.166152*((2-y(3))*k(3)*0.7^2*yMEOH*yO2+y(3)*k(1)*0.7^2.5*yMEOH^2*yO2^0.5); 0.05832*(2+y(1)-y(3))/4*0.166152*(2*k(4)*0.7^1*yMEOH+y(4)*k(1)*0.7^2.5*yMEOH^2*yO2^0.5-y(4)*k(3)*0.7^2*yMEOH*yO2); ]; % dydl=... % [A*Rho(2+y(1)-y(3))/4*nT0*((2+y(1))*k(1)*P^2*yMEOH^1*yO2^1-y(1)*k(3)*P^2*yMEOH*yO2); % A*Rho(2+y(1)-y(3))/4*nT0*(2*k(2)*P^2*yMEOH*yO2+y(2)*k(1)*P^2*yMEOH^1*yO2^1-y(2)*k(3)*P^2*yMEOH*yO2); % A*Rho(2+y(1)-y(3))/4*nT0*((2-y(3))*k(3)*P^2*yMEOH*yO2+y(3)*k(1)*P^2*yMEOH^1*yO2^1); % A*Rho(2+y(1)-y(3))/4*nT0*(2*k(4)*P^2*yMEOH+y(4)*k(1)*P^2*yMEOH^1*yO2^1-y(4)*k(3)*P^2*yMEOH*yO2); function dydl = KineticsEqs5(l,y,k) %P=[0.6 0.8 0.5 0.7 0.9];P为总压 %nT0=[0.154284 0.161404 0.161404 0.166152 0.332303 ]; % n:进料总摩尔数 %yO20=[0.0769 0.0882 0.1176 0.1429 0.0857 ]; %进料O2组成 %yMEOH0=[0.3077 0.2941 0.2941 0.2857 0.1429 ]; %进料甲醇（MEOH）组成 %A*Rho=0.05832 yO2=0.0857*(1+0.5*y(1)-0.5*y(3))-(0.5*y(1)+y(2)+0.5*y(3)); yMEOH=0.1429*(1+0.5*y(1)-0.5*y(3))-(3*y(1)+2*y(2)+y(3)+2*y(4)); % yO2=yO20*(1+0.5*y(1)-0.5*y(3))-(0.5*y(1)+y(2)+0.5*y(3)); % yMEOH=yMEOH0*(1+0.5*y(1)-0.5*y(3))-(3*y(1)+2*y(2)+y(3)+2*y(4)); dydl=... [0.05832*(2+y(1)-y(3))/4*0.332303*((2+y(1))*k(1)*0.9^2.5*yMEOH^2*yO2^0.5-y(1)*k(3)*0.9^2*yMEOH*yO2); 0.05832*(2+y(1)-y(3))/4*0.332303*(2*k(2)*0.9^2*yMEOH^1*yO2^1+y(2)*k(1)*0.9^2.5*yMEOH^2*yO2^0.5-y(2)*k(3)*0.9^2*yMEOH*yO2); 0.05832*(2+y(1)-y(3))/4*0.332303*((2-y(3))*k(3)*0.9^2*yMEOH*yO2+y(3)*k(1)*0.9^2.5*yMEOH^2*yO2^0.5); 0.05832*(2+y(1)-y(3))/4*0.332303*(2*k(4)*0.9^1*yMEOH+y(4)*k(1)*0.9^2.5*yMEOH^2*yO2^0.5-y(4)*k(3)*0.9^2*yMEOH*yO2); ]; % dydl=... % [A*Rho(2+y(1)-y(3))/4*nT0*((2+y(1))*k(1)*P^2*yMEOH^1*yO2^1-y(1)*k(3)*P^2*yMEOH*yO2); % A*Rho(2+y(1)-y(3))/4*nT0*(2*k(2)*P^2*yMEOH*yO2+y(2)*k(1)*P^2*yMEOH^1*yO2^1-y(2)*k(3)*P^2*yMEOH*yO2); % A*Rho(2+y(1)-y(3))/4*nT0*((2-y(3))*k(3)*P^2*yMEOH*yO2+y(3)*k(1)*P^2*yMEOH^1*yO2^1); % A*Rho(2+y(1)-y(3))/4*nT0*(2*k(4)*P^2*yMEOH+y(4)*k(1)*P^2*yMEOH^1*yO2^1-y(4)*k(3)*P^2*yMEOH*yO2); |
2楼2015-10-26 09:10:18