| 3 | 1/1 | 返回列表 |
| 查看: 783 | 回復(fù): 2 | ||
cyzeng金蟲 (小有名氣)
|
[求助]
縮聚反應(yīng)動(dòng)力學(xué)參數(shù)估計(jì) 置信區(qū)間太大 求程序改進(jìn)
|
|
縮聚反應(yīng)動(dòng)力學(xué)參數(shù)估計(jì) 置信區(qū)間太大 求程序改進(jìn) 程序 clc;clear; format long global K1 K2 K3 cin rin T R T=[483.15 493.15 503.15 513.15 523.15];%動(dòng)力學(xué)實(shí)驗(yàn)溫度點(diǎn),K R=8.3145;%J/(K*mol) lb=[1e2 3e4 1e2 3e4 1e2 3e4]; ub=[1e10 8e4 1e10 8e4 1e10 8e4]; %實(shí)驗(yàn)測(cè)定平衡常數(shù)隨溫度變化 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; % 動(dòng)力學(xué)數(shù)據(jù)-不同溫度下間歇反應(yīng)器內(nèi)各濃度隨時(shí)間變化 % 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]; %生成動(dòng)力學(xué)數(shù)據(jù)三維矩陣-轉(zhuǎn)換成行向量 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, ; %濃度插值及微分求導(dǎo)得反應(yīng)速率 for i=1:5 [cin(:,:,i) rin(:,:,i)]=Rate(time(:,:,i),cexp(:,:,i)); end cin; rin; %擬合EG和H2O反應(yīng)速率曲線 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進(jìn)行參數(shù)估計(jì) % options=optimset('MaxFunEvals',2400); [b,fval,flag]=fmincon(@objfuncfmincon,b0,[],[],[],[],lb,ub,[],[],c0,cin,ppp2,ppp5); fval b_fmincon=b; b0=b_fmincon %以fmincon估計(jì)的k0,使用lsqnonlin進(jìn)行參數(shù)估計(jì) 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)) %模型計(jì)算值與實(shí)驗(yàn)點(diǎn)比較 %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') %定義插值和反應(yīng)速率函數(shù) 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]; 結(jié)果 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 |
金蟲 (小有名氣)
|
本帖內(nèi)容被屏蔽 |
| 3 | 1/1 | 返回列表 |
| 最具人氣熱帖推薦 [查看全部] | 作者 | 回/看 | 最后發(fā)表 | |
|---|---|---|---|---|
|
[考研] 289求調(diào)劑 +5 | BrightLL 2026-03-29 | 5/250 |
|
|---|---|---|---|---|
|
[考研] 275求調(diào)劑 +13 | Micky11223 2026-03-25 | 18/900 |
|
|
[考研] 070305高分子化學(xué)與物理 304分求調(diào)劑 +12 | c297914 2026-03-28 | 12/600 |
|
|
[考研] 材料學(xué)碩333求調(diào)劑 +11 | 北道巷 2026-03-24 | 11/550 |
|
|
[考研] 調(diào)劑求院校招收 +6 | 鶴鯨鴿 2026-03-28 | 6/300 |
|
|
[考研] 0703化學(xué)調(diào)劑,求導(dǎo)師收 +9 | 天天好運(yùn)來(lái)上岸?/a> 2026-03-24 | 10/500 |
|
|
[考研] 11408軟件工程求調(diào)劑 +3 | Qiu學(xué)ing 2026-03-28 | 3/150 |
|
|
[考研] 312,生物學(xué)求調(diào)劑 +3 | 小譯同學(xué)abc 2026-03-28 | 3/150 |
|
|
[考研] 340求調(diào)劑 +5 | jhx777 2026-03-27 | 5/250 |
|
|
[考研] 0703化學(xué)求調(diào)劑,各位老師看看我。。 +5 | 祁祺祺 2026-03-25 | 5/250 |
|
|
[考研] 一志愿南師大0703化學(xué) 275求調(diào)劑 +4 | Ripcord上岸 2026-03-27 | 4/200 |
|
|
[考研] 298調(diào)劑 +3 | jiyingjie123 2026-03-27 | 3/150 |
|
|
[考研] 求調(diào)劑 一志愿 本科 北科大 化學(xué) 343 +6 | 13831862839 2026-03-24 | 7/350 |
|
|
[考研] 294分080500材料科學(xué)與工程求調(diào)劑 +4 | 柳溪邊 2026-03-26 | 4/200 |
|
|
[考研] 321求調(diào)劑 +6 | wasdssaa 2026-03-26 | 6/300 |
|
|
[考研] 機(jī)械學(xué)碩310分,數(shù)一英一,一志愿211本科雙非找調(diào)劑信息 +3 | @357 2026-03-25 | 3/150 |
|
|
[考研] 材料科學(xué)與工程 317求調(diào)劑 +4 | JKSOIID 2026-03-26 | 4/200 |
|
|
[考研] 求調(diào)劑 +6 | 研研,接電話 2026-03-24 | 7/350 |
|
|
[考研] 292求調(diào)劑 +4 | 鵝鵝鵝額額額額?/a> 2026-03-24 | 4/200 |
|
|
[考研] 求老師收我 +3 | zzh16938784 2026-03-23 | 3/150 |
|