24小时热门版块排行榜

返回列表

cyzeng

金虫 (小有名气)

应助: 4 (幼儿园)
金币: 1515.9
散金: 4
红花: 2
帖子: 73
在线: 33.5小时
虫号: 704956
注册: 2009-02-20
性别: MM
专业: 无机纳米化学

[求助] 缩聚反应动力学参数估计置信区间太大求程序改进

缩聚反应动力学参数估计置信区间太大求程序改进
程序
clc;clear;
format long
global K1 K2 K3 cin rin T R
T=[483.15 493.15 503.15 513.15 523.15];%动力学实验温度点，K
R=8.3145;%J/(K*mol)
lb=[1e2 3e4  1e2 3e4 1e2 3e4];
ub=[1e10  8e4  1e10  8e4 1e10  8e4];

%实验测定平衡常数随温度变化
deltS1=43.2237;
deltH1=1.58e+4;
deltS2=50.4599;
deltH2=2.01e+4;
K1=exp(deltS1/R-deltH1./(R*T));
K2=exp(deltS2/R-deltH2./(R*T));
K3=K2./K1;
% 动力学数据-不同温度下间歇反应器内各浓度随时间变化
%          time/h c1    c2    c3    c4    c5
ExpData_210=[0    0.372 0.076 1.630 4.439 0.0249
         1    0.146 0.087 1.339 4.718 0.0127
         2    0.059 0.073 1.270 4.820 0.0100
         3    0.025 0.070 1.245 4.856 0.0052];
ExpData_220=[0    0.400 0.079 1.683 4.387 0.0306;
         1    0.125 0.054 1.016 4.990 0.0106;
         2    0.042 0.033 0.948 5.100 0.0067;
         3    0.015 0.028 0.886 5.162 0.0033];
ExpData_230=[0    0.290 0.070 1.620 4.488 0.0239
         0.5 0.102 0.061 1.071 4.954 0.0099
         1    0.040 0.039 1.033 5.035 0.0054
         2    0.008 0.036 1.011 5.385 0.0030];
ExpData_240=[0    0.285 0.065 1.386 4.656 0.0229
         0.5 0.088 0.033 0.823 5.165 0.0086
         1    0.035 0.028 0.781 5.224 0.0043
         1.5 0.016 0.025 0.750 5.258 0.0020];
ExpData_250=[0    0.228 0.050 1.090 4.902 0.0236
         0.5 0.054 0.018 0.635 5.329 0.0086
         0.75 0.038 0.016 0.608 5.357 0.0064
         1    0.024 0.014 0.589 5.378 0.0041];
%生成动力学数据三维矩阵-转换成行向量
ExpData=ones(6,4,5);
ExpData(:,:,1)=ExpData_210';
ExpData(:,:,2)=ExpData_220';
ExpData(:,:,3)=ExpData_230';
ExpData(:,:,4)=ExpData_240';
ExpData(:,:,5)=ExpData_250';
time=ExpData(1,:,

;
cexp=ExpData(2:6,:,

;
c0=cexp(:,1,

;
%浓度插值及微分求导得反应速率
for i=1:5
[cin(:,:,i) rin(:,:,i)]=Rate(time(:,:,i),cexp(:,:,i));
end
cin;
rin;
%拟合EG和H2O反应速率曲线
for i=1:5
ti2=linspace(time(:,1,i),time(:,4,i),61);
ppp2(:,:,i)=polyfit(ti2,rin(2,:,i),4);
ppp5(:,:,i)=polyfit(ti2,rin(5,:,i),4);
end
ppp2;ppp5;

b0=[2000.4 30547.1 16289.2 42486.9 13232.6 48038.2];

%使用fmincon进行参数估计
% options=optimset('MaxFunEvals',2400);
[b,fval,flag]=fmincon(@objfuncfmincon,b0,[],[],[],[],lb,ub,[],[],c0,cin,ppp2,ppp5);
fval
b_fmincon=b;
b0=b_fmincon
%以fmincon估计的k0,使用lsqnonlin进行参数估计
options=optimset('MaxFunEvals',2400);
[b,resnorm,resid,exitflag,output,lambda,jacobian]=lsqnonlin(@objfunclsqnonlin,b0,lb,ub,options,c0,cin,ppp2,ppp5);
resnorm
b
bi=nlparci(b,resid,jacobian);
fprintf('\tk10=%.1f±%.1f\n',b(1),bi(1,2)-b(1))
fprintf('\tEa1=%.1f±%.1f\n',b(2),bi(2,2)-b(2))
fprintf('\tk30=%.1f±%.1f\n',b(3),bi(3,2)-b(3))
fprintf('\tEa3=%.1f±%.1f\n',b(4),bi(4,2)-b(4))
fprintf('\tk50=%.1f±%.1f\n',b(5),bi(5,2)-b(5))
fprintf('\tEa5=%.1f±%.1f\n',b(6),bi(6,2)-b(6))
%模型计算值与实验点比较
%210度
ti3=linspace(time(:,1,1),time(:,4,1),61);
[t,c]=ode45(@masslsq1,ti3,c0(:,:,1),[],b,ppp2,ppp5);
figure
plot(time(:,:,1),cexp(1,:,1),'ks',time(:,:,1),cexp(2,:,1),'k<',time(:,:,1),cexp(3,:,1),'k>',time(:,:,1),cexp(4,:,1),'k*',time(:,:,1),cexp(5,:,1),'ko',...
ti3,c(:,1)','k-',ti3,c(:,2)','k-',ti3,c(:,3)','k-',ti3,c(:,4)','k-',ti3,c(:,5)','k-')
legend('c_1(-COOH)','c_2(EG)','c_3(-OH)','c_4(-COOC-)','c_5(H_2O)')
legend('boxoff')
xlabel('t,h'),ylabel('c,mol/kg'),text(1,4,'T=483.15K')
%220度
ti3=linspace(time(:,1,2),time(:,4,2),61);
[t,c]=ode45(@masslsq2,ti3,c0(:,:,2),[],b,ppp2,ppp5);
figure
plot(time(:,:,2),cexp(1,:,2),'ks',time(:,:,2),cexp(2,:,2),'k<',time(:,:,2),cexp(3,:,2),'k>',time(:,:,2),cexp(4,:,2),'k*',time(:,:,2),cexp(5,:,2),'ko',...
ti3,c(:,1)','k-',ti3,c(:,2)','k-',ti3,c(:,3)','k-',ti3,c(:,4)','k-',ti3,c(:,5)','k-')
legend('c_1(-COOH)','c_2(EG)','c_3(-OH)','c_4(-COOC-)','c_5(H_2O)')
legend('boxoff')
xlabel('t,h'),ylabel('c,mol/kg'),text(1,4,'T=493.15K')
%230度
ti3=linspace(time(:,1,3),time(:,4,3),61);
[t,c]=ode45(@masslsq3,ti3,c0(:,:,3),[],b,ppp2,ppp5);
figure
plot(time(:,:,3),cexp(1,:,3),'ks',time(:,:,3),cexp(2,:,3),'k<',time(:,:,3),cexp(3,:,3),'k>',time(:,:,3),cexp(4,:,3),'k*',time(:,:,3),cexp(5,:,3),'ko',...
ti3,c(:,1)','k-',ti3,c(:,2)','k-',ti3,c(:,3)','k-',ti3,c(:,4)','k-',ti3,c(:,5)','k-')
legend('c_1(-COOH)','c_2(EG)','c_3(-OH)','c_4(-COOC-)','c_5(H_2O)')
legend('boxoff')
xlabel('t,h'),ylabel('c,mol/kg'),text(1,4,'T=503.15K')
%240度
ti3=linspace(time(:,1,4),time(:,4,4),61);
[t,c]=ode45(@masslsq4,ti3,c0(:,:,4),[],b,ppp2,ppp5);
figure
plot(time(:,:,4),cexp(1,:,4),'ks',time(:,:,4),cexp(2,:,4),'k<',time(:,:,4),cexp(3,:,4),'k>',time(:,:,4),cexp(4,:,4),'k*',time(:,:,4),cexp(5,:,4),'ko',...
ti3,c(:,1)','k-',ti3,c(:,2)','k-',ti3,c(:,3)','k-',ti3,c(:,4)','k-',ti3,c(:,5)','k-')
legend('c_1(-COOH)','c_2(EG)','c_3(-OH)','c_4(-COOC-)','c_5(H_2O)')
legend('boxoff')
xlabel('t,h'),ylabel('c,mol/kg'),text(1,4,'T=513.15K')
%250度
ti3=linspace(time(:,1,5),time(:,4,5),61);
[t,c]=ode45(@masslsq5,ti3,c0(:,:,5),[],b,ppp2,ppp5);
figure
plot(time(:,:,5),cexp(1,:,5),'ks',time(:,:,5),cexp(2,:,5),'k<',time(:,:,5),cexp(3,:,5),'k>',time(:,:,5),cexp(4,:,5),'k*',time(:,:,5),cexp(5,:,5),'ko',...
ti3,c(:,1)','k-',ti3,c(:,2)','k-',ti3,c(:,3)','k-',ti3,c(:,4)','k-',ti3,c(:,5)','k-')
legend('c_1(-COOH)','c_2(EG)','c_3(-OH)','c_4(-COOC-)','c_5(H_2O)')
legend('boxoff')
xlabel('t,h'),ylabel('c,mol/kg'),text(1,4,'T=523.15K')

%定义插值和反应速率函数
function [cin rin]=Rate(t,c)
% knots=1; K=4;
% sp=spap2(knots,K,t,c);
% pp=fnder(sp);
% ti=linspace(t(1),t(end),61);
% cin=fnval(sp,ti);
% rin=fnval(pp,ti);
sp=csaps(t,c,0.99);
pp=fnder(sp);
ti=linspace(t(1),t(end),61);
cin=fnval(sp,ti);
rin=fnval(pp,ti);

%定义objfuncfmincon
function ffmi=objfuncfmincon(b,c0,cin,ppp2,ppp5)
tspan=linspace(0,3,61);
[t,c]=ode45(@masslsq1,tspan,c0(:,:,1),[],b,ppp2,ppp5);
f1=sum((c(:,1)-cin(1,:,1)').^2)+sum((c(:,3)-cin(3,:,1)').^2)+sum((c(:,4)-cin(4,:,1)').^2);
[t,c]=ode45(@masslsq2,tspan,c0(:,:,2),[],b,ppp2,ppp5);
f2=sum((c(:,1)-cin(1,:,2)').^2)+sum((c(:,3)-cin(3,:,2)').^2)+sum((c(:,4)-cin(4,:,2)').^2);
tspan2=linspace(0,2,61);
[t,c]=ode45(@masslsq3,tspan2,c0(:,:,3),[],b,ppp2,ppp5);
f3=sum((c(:,1)-cin(1,:,3)').^2)+sum((c(:,3)-cin(3,:,3)').^2)+sum((c(:,4)-cin(4,:,3)').^2);
tspan3=linspace(0,1.5,61);
[t,c]=ode45(@masslsq4,tspan3,c0(:,:,4),[],b,ppp2,ppp5);
f4=sum((c(:,1)-cin(1,:,4)').^2)+sum((c(:,3)-cin(3,:,4)').^2)+sum((c(:,4)-cin(4,:,4)').^2);
tspan4=linspace(0,1,61);
[t,c]=ode45(@masslsq5,tspan4,c0(:,:,5),[],b,ppp2,ppp5);
f5=sum((c(:,1)-cin(1,:,5)').^2)+sum((c(:,3)-cin(3,:,5)').^2)+sum((c(:,4)-cin(4,:,5)').^2);
ffmi=f1+f2+f3+f4+f5;

%定义objfunclsqnonlin
function flsq=objfunclsqnonlin(b,c0,cin,ppp2,ppp5)
tspan=linspace(0,3,61);
[t,c]=ode45(@masslsq1,tspan,c0(:,:,1),[],b,ppp2,ppp5);
f1=c(:,1)-cin(1,:,1)';
f2=c(:,3)-cin(3,:,1)';
f3=c(:,4)-cin(4,:,1)';
[t,c]=ode45(@masslsq2,tspan,c0(:,:,2),[],b,ppp2,ppp5);
f4=c(:,1)-cin(1,:,2)';
f5=c(:,3)-cin(3,:,2)';
f6=c(:,4)-cin(4,:,2)';
tspan2=linspace(0,2,61);
[t,c]=ode45(@masslsq3,tspan2,c0(:,:,3),[],b,ppp2,ppp5);
f7=c(:,1)-cin(1,:,3)';
f8=c(:,3)-cin(3,:,3)';
f9=c(:,4)-cin(4,:,3)';
tspan3=linspace(0,1.5,61);
[t,c]=ode45(@masslsq4,tspan3,c0(:,:,4),[],b,ppp2,ppp5);
f10=c(:,1)-cin(1,:,4)';
f11=c(:,3)-cin(3,:,4)';
f12=c(:,4)-cin(4,:,4)';
tspan4=linspace(0,1,61);
[t,c]=ode45(@masslsq5,tspan4,c0(:,:,5),[],b,ppp2,ppp5);
f13=c(:,1)-cin(1,:,5)';
f14=c(:,3)-cin(3,:,5)';
f15=c(:,4)-cin(4,:,5)';
flsq=[f1;f2;f3;f4;f5;f6;f7;f8;f9;f10;f11;f12;f13;f14;f15];

%定义masslsq1-210度
function dcdt1=masslsq1(t,c,b,ppp2,ppp5)
global K1 K2 K3 T R
k1=b(1)*exp(-b(2)/(R*T(1)));
k2=k1/K1(1);
k3=b(3)*exp(-b(4)/(R*T(1)));
k4=k3/K2(1);
k5=b(5)*exp(-b(6)/(R*T(1)));
k6=k5/K3(1);
dc1dt=-2*k1*c(1)*c(2)+k2*c(3)*c(5)-k3*c(1)*c(3)+2*k4*c(4)*c(5);
dc2dt=ppp2(:,5,1)+ppp2(:,4,1)*t+ppp2(:,3,1)*t.^2+ppp2(:,2,1)*t.^3+ppp2(:,1,1)*t.^4;
dc3dt=2*k1*c(1)*c(2)-k2*c(3)*c(5)-k3*c(1)*c(3)+2*k4*c(4)*c(5)-2*k5*c(3)*c(3)+8*k6*c(2)*c(4);
dc4dt=k3*c(1)*c(3)-2*k4*c(4)*c(5)+k5*c(3)*c(3)-4*k6*c(2)*c(4);
dc5dt=ppp5(:,5,1)+ppp5(:,4,1)*t+ppp5(:,3,1)*t.^2+ppp5(:,2,1)*t.^3+ppp5(:,1,1)*t.^4;
dcdt1=[dc1dt;dc2dt;dc3dt;dc4dt;dc5dt];

%定义masslsq2-220度
function dcdt2=masslsq2(t,c,b,ppp2,ppp5)
global K1 K2 K3 T R
k1=b(1)*exp(-b(2)/(R*T(2)));
k2=k1/K1(2);
k3=b(3)*exp(-b(4)/(R*T(2)));
k4=k3/K2(2);
k5=b(5)*exp(-b(6)/(R*T(2)));
k6=k5/K3(2);
dc1dt=-2*k1*c(1)*c(2)+k2*c(3)*c(5)-k3*c(1)*c(3)+2*k4*c(4)*c(5);
dc2dt=ppp2(:,5,2)+ppp2(:,4,2)*t+ppp2(:,3,2)*t.^2+ppp2(:,2,2)*t.^3+ppp2(:,1,2)*t.^4;
dc3dt=2*k1*c(1)*c(2)-k2*c(3)*c(5)-k3*c(1)*c(3)+2*k4*c(4)*c(5)-2*k5*c(3)*c(3)+8*k6*c(2)*c(4);
dc4dt=k3*c(1)*c(3)-2*k4*c(4)*c(5)+k5*c(3)*c(3)-4*k6*c(2)*c(4);
dc5dt=ppp5(:,5,2)+ppp5(:,4,2)*t+ppp5(:,3,2)*t.^2+ppp5(:,2,2)*t.^3+ppp5(:,1,2)*t.^4;
dcdt2=[dc1dt;dc2dt;dc3dt;dc4dt;dc5dt];

%定义masslsq3-230度
function dcdt3=masslsq3(t,c,b,ppp2,ppp5)
global K1 K2 K3 T R
k1=b(1)*exp(-b(2)/(R*T(3)));
k2=k1/K1(3);
k3=b(3)*exp(-b(4)/(R*T(3)));
k4=k3/K2(3);
k5=b(5)*exp(-b(6)/(R*T(3)));
k6=k5/K3(3);
dc1dt=-2*k1*c(1)*c(2)+k2*c(3)*c(5)-k3*c(1)*c(3)+2*k4*c(4)*c(5);
dc2dt=ppp2(:,5,3)+ppp2(:,4,3)*t+ppp2(:,3,3)*t.^2+ppp2(:,2,3)*t.^3+ppp2(:,1,3)*t.^4;
dc3dt=2*k1*c(1)*c(2)-k2*c(3)*c(5)-k3*c(1)*c(3)+2*k4*c(4)*c(5)-2*k5*c(3)*c(3)+8*k6*c(2)*c(4);
dc4dt=k3*c(1)*c(3)-2*k4*c(4)*c(5)+k5*c(3)*c(3)-4*k6*c(2)*c(4);
dc5dt=ppp5(:,5,3)+ppp5(:,4,3)*t+ppp5(:,3,3)*t.^2+ppp5(:,2,3)*t.^3+ppp5(:,1,3)*t.^4;
dcdt3=[dc1dt;dc2dt;dc3dt;dc4dt;dc5dt];

%定义masslsq4-240度
function dcdt4=masslsq4(t,c,b,ppp2,ppp5)
global K1 K2 K3 T R
k1=b(1)*exp(-b(2)/(R*T(4)));
k2=k1/K1(4);
k3=b(3)*exp(-b(4)/(R*T(4)));
k4=k3/K2(4);
k5=b(5)*exp(-b(6)/(R*T(4)));
k6=k5/K3(4);
dc1dt=-2*k1*c(1)*c(2)+k2*c(3)*c(5)-k3*c(1)*c(3)+2*k4*c(4)*c(5);
dc2dt=ppp2(:,5,4)+ppp2(:,4,4)*t+ppp2(:,3,4)*t.^2+ppp2(:,2,4)*t.^3+ppp2(:,1,4)*t.^4;
dc3dt=2*k1*c(1)*c(2)-k2*c(3)*c(5)-k3*c(1)*c(3)+2*k4*c(4)*c(5)-2*k5*c(3)*c(3)+8*k6*c(2)*c(4);
dc4dt=k3*c(1)*c(3)-2*k4*c(4)*c(5)+k5*c(3)*c(3)-4*k6*c(2)*c(4);
dc5dt=ppp5(:,5,4)+ppp5(:,4,4)*t+ppp5(:,3,4)*t.^2+ppp5(:,2,4)*t.^3+ppp5(:,1,4)*t.^4;
dcdt4=[dc1dt;dc2dt;dc3dt;dc4dt;dc5dt];

%定义masslsq5-250度
function dcdt5=masslsq5(t,c,b,ppp2,ppp5)
global K1 K2 K3 T R
k1=b(1)*exp(-b(2)/(R*T(5)));
k2=k1/K1(5);
k3=b(3)*exp(-b(4)/(R*T(5)));
k4=k3/K2(5);
k5=b(5)*exp(-b(6)/(R*T(5)));
k6=k5/K3(5);
dc1dt=-2*k1*c(1)*c(2)+k2*c(3)*c(5)-k3*c(1)*c(3)+2*k4*c(4)*c(5);
dc2dt=ppp2(:,5,5)+ppp2(:,4,5)*t+ppp2(:,3,5)*t.^2+ppp2(:,2,5)*t.^3+ppp2(:,1,5)*t.^4;
dc3dt=2*k1*c(1)*c(2)-k2*c(3)*c(5)-k3*c(1)*c(3)+2*k4*c(4)*c(5)-2*k5*c(3)*c(3)+8*k6*c(2)*c(4);
dc4dt=k3*c(1)*c(3)-2*k4*c(4)*c(5)+k5*c(3)*c(3)-4*k6*c(2)*c(4);
dc5dt=ppp5(:,5,5)+ppp5(:,4,5)*t+ppp5(:,3,5)*t.^2+ppp5(:,2,5)*t.^3+ppp5(:,1,5)*t.^4;
dcdt5=[dc1dt;dc2dt;dc3dt;dc4dt;dc5dt];

结果
Optimization terminated: directional derivative predicts change in
  objective value less than options.TolFun and maximum constraint violation
  is less than options.TolCon.
No active inequalities.

fval =

9.975490368849357

b0 =

  1.0e+004 *

0.021390412595752 3.140646199334575 1.805623620390919 3.561747395742712 1.562837065436626 4.047673875392536

Optimization terminated: relative function value
changing by less than OPTIONS.TolFun.

resnorm =

8.289356826509787

b =

  1.0e+007 *

0.000824325058983 0.003320473746231 0.141719950795143 0.005636884012449 3.314837171169622 0.006838741499027

k10=8243.3±936367.3
Ea1=33204.7±461448.7
k30=1417199.5±20061412.4
Ea3=56368.8±58203.6
k50=33148371.7±584453034.8
Ea5=68387.4±72746.8

回复此楼