cupczp
2012-03-16, 15:07
function Fourlump
%由实验室固定床反应器数据估计三集总催化裂化动力学模型参数
format long e
k0=[1.0 0.5 1.0];%动力学参数初值
[k,resnorm,residual,exitflag]=lsqnonlin(@ObjFun,k0)
zspan=[0 1];
function Delta=ObjFun(k)
%参数估计目标函数
global xx0 M NH V21 V23 V24 theta A
xx0==[31.423 42.861 6.403 19.315];
M=[93.99229 83.27657 102.8053 116.2428];
NH=[0.2271 0.2271 0.2271 0.0852 0.0852 0.0852 0.1135 0.1135 0.1419 0.1419 0.1419 0.1419 0.1703 0.1703 0.1703];
V21=[0.8837 1.0897 0.9285 0.8446 0.9119 0.8925 0.9612 0.8992 0.8601 0.9473 0.9115 0.9747 0.8795 0.9365 0.9162];
V23=[0.8853 0.9914 0.8874 0.8228 0.8501 0.8609 0.9220 0.8669 0.8354 0.8654 0.8592 0.9637 0.8455 0.8638 0.9014];
V24=[0.7570 0.8668 0.7843 0.7332 0.7551 0.7666 0.7936 0.7678 0.7411 0.7644 0.7638 0.8576 0.7503 0.7658 0.8033];
theta=[166.8411 250.2616 286.0133 155.4816 248.7706 266.5399 169.4721 593.1523 80.9890 377.9487 566.9230 566.9229 180.9717 361.9435 434.3321];
A=[0.1562 0.1562 0.1562 0.0728 0.0728 0.0728 0.0925 0.0925 0.1106 0.1106 0.1106 0.1106 0.1270 0.1270 0.1270];
%--------------------------------------试验产率
%烷烃产率实验值
Pexp=[0.002553 0.003354 0.002758 0.003226 0.003495 0.003105 0.003571 0.003073 0.002857 0.003621 0.003228 0.003121 0.002850 0.003406 0.003142];
%烯烃产率实验值
Oexp=[0.002189 0.001745 0.002128 0.003156 0.002759 0.002974 0.001938 0.002763 0.002761 0.001888 0.002415 0.001553 0.002348 0.001851 0.001809];
%环烷烃产率实验值
Nexp=[0.000629 0.000684 0.000585 0.000672 0.000707 0.000678 0.000791 0.000675 0.000630 0.000679 0.000629 0.000734 0.000636 0.000648 0.000677];
%芳烃产率实验值
Aexp=[0.001467 0.001973 0.001336 0.001585 0.001658 0.001677 0.001780 0.001594 0.001569 0.001637 0.001522 0.001877 0.001601 0.001570 0.001687];
%--------------------------------------
for i=1:length(theta)
f0(1)=xx0(1)/((1+2*NH(i))*M(1));%计算进料浓度：mole oil/g gas
f0(2)=xx0(2)/((1+2*NH(i))*M(2));
f0(3)=xx0(3)/((1+2*NH(i))*M(3));
f0(4)=xx0(4)/((1+2*NH(i))*M(4));
z0=0;z1=1+eps;
[z,f]=ode45(@FixBed,[z0 z1],f0,[],k,theta(i),V21(i),V23(i),V24(i),A(i),NH(i));%估计产率
fout=f(end,:);
Pout(i)=fout(1);
Oout(i)=fout(2);
Nout(i)=fout(3);
Aout(i)=fout(4);
end
DeltaG=Pexp-Pout;
DeltaC=Oexp-Oout;
DeltaG=Nexp-Nout;
DeltaC=Aexp-Aout;
Delta=[DeltaP;DeltaO;DeltaN;DeltaA];
%--------------------------------------
function df=FixBed(z,f,k)
%固定床反应器物料衡算微分方程
global xx0 M NH V21 V23 V24 theta A
Mwave=1.0/(f(1)+f(2)+f(3)+f(4)+A);
df(1)=theta*Mwave*k(1)*V21*f(2);
df(2)=-theta*Mwave*(k(1)+k(2)+k(3))*f(2);
df(3)=theta*Mwave*k(2)*V23*f(2);
df(4)=theta*Mwave*k(3)*V24*f(2);
df=df';
运行后：??? Error using ==> optim\private\lsqncommon
User supplied function failed with the following error:
Error using ==> eq
Matrix dimensions must agree.
Error in ==> lsqnonlin at 147
[x,Resnorm,FVAL,EXITFLAG,OUTPUT,LAMBDA,JACOB] = ...
Error in ==> Fourlump11111 at 5
[k,resnorm,residual,exitflag]=lsqnonlin(@ObjFun,k0)
%由实验室固定床反应器数据估计三集总催化裂化动力学模型参数
format long e
k0=[1.0 0.5 1.0];%动力学参数初值
[k,resnorm,residual,exitflag]=lsqnonlin(@ObjFun,k0)
zspan=[0 1];
function Delta=ObjFun(k)
%参数估计目标函数
global xx0 M NH V21 V23 V24 theta A
xx0==[31.423 42.861 6.403 19.315];
M=[93.99229 83.27657 102.8053 116.2428];
NH=[0.2271 0.2271 0.2271 0.0852 0.0852 0.0852 0.1135 0.1135 0.1419 0.1419 0.1419 0.1419 0.1703 0.1703 0.1703];
V21=[0.8837 1.0897 0.9285 0.8446 0.9119 0.8925 0.9612 0.8992 0.8601 0.9473 0.9115 0.9747 0.8795 0.9365 0.9162];
V23=[0.8853 0.9914 0.8874 0.8228 0.8501 0.8609 0.9220 0.8669 0.8354 0.8654 0.8592 0.9637 0.8455 0.8638 0.9014];
V24=[0.7570 0.8668 0.7843 0.7332 0.7551 0.7666 0.7936 0.7678 0.7411 0.7644 0.7638 0.8576 0.7503 0.7658 0.8033];
theta=[166.8411 250.2616 286.0133 155.4816 248.7706 266.5399 169.4721 593.1523 80.9890 377.9487 566.9230 566.9229 180.9717 361.9435 434.3321];
A=[0.1562 0.1562 0.1562 0.0728 0.0728 0.0728 0.0925 0.0925 0.1106 0.1106 0.1106 0.1106 0.1270 0.1270 0.1270];
%--------------------------------------试验产率
%烷烃产率实验值
Pexp=[0.002553 0.003354 0.002758 0.003226 0.003495 0.003105 0.003571 0.003073 0.002857 0.003621 0.003228 0.003121 0.002850 0.003406 0.003142];
%烯烃产率实验值
Oexp=[0.002189 0.001745 0.002128 0.003156 0.002759 0.002974 0.001938 0.002763 0.002761 0.001888 0.002415 0.001553 0.002348 0.001851 0.001809];
%环烷烃产率实验值
Nexp=[0.000629 0.000684 0.000585 0.000672 0.000707 0.000678 0.000791 0.000675 0.000630 0.000679 0.000629 0.000734 0.000636 0.000648 0.000677];
%芳烃产率实验值
Aexp=[0.001467 0.001973 0.001336 0.001585 0.001658 0.001677 0.001780 0.001594 0.001569 0.001637 0.001522 0.001877 0.001601 0.001570 0.001687];
%--------------------------------------
for i=1:length(theta)
f0(1)=xx0(1)/((1+2*NH(i))*M(1));%计算进料浓度：mole oil/g gas
f0(2)=xx0(2)/((1+2*NH(i))*M(2));
f0(3)=xx0(3)/((1+2*NH(i))*M(3));
f0(4)=xx0(4)/((1+2*NH(i))*M(4));
z0=0;z1=1+eps;
[z,f]=ode45(@FixBed,[z0 z1],f0,[],k,theta(i),V21(i),V23(i),V24(i),A(i),NH(i));%估计产率
fout=f(end,:);
Pout(i)=fout(1);
Oout(i)=fout(2);
Nout(i)=fout(3);
Aout(i)=fout(4);
end
DeltaG=Pexp-Pout;
DeltaC=Oexp-Oout;
DeltaG=Nexp-Nout;
DeltaC=Aexp-Aout;
Delta=[DeltaP;DeltaO;DeltaN;DeltaA];
%--------------------------------------
function df=FixBed(z,f,k)
%固定床反应器物料衡算微分方程
global xx0 M NH V21 V23 V24 theta A
Mwave=1.0/(f(1)+f(2)+f(3)+f(4)+A);
df(1)=theta*Mwave*k(1)*V21*f(2);
df(2)=-theta*Mwave*(k(1)+k(2)+k(3))*f(2);
df(3)=theta*Mwave*k(2)*V23*f(2);
df(4)=theta*Mwave*k(3)*V24*f(2);
df=df';
运行后：??? Error using ==> optim\private\lsqncommon
User supplied function failed with the following error:
Error using ==> eq
Matrix dimensions must agree.
Error in ==> lsqnonlin at 147
[x,Resnorm,FVAL,EXITFLAG,OUTPUT,LAMBDA,JACOB] = ...
Error in ==> Fourlump11111 at 5
[k,resnorm,residual,exitflag]=lsqnonlin(@ObjFun,k0)