|
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); |
|