24小时热门版块排行榜

返回列表

kidingkids

新虫 (初入文坛)

应助: 0 (幼儿园)
金币: 2.5
散金: 10
帖子: 5
在线: 4.1小时
虫号: 1292152
注册: 2011-05-11
性别: GG
专业: 化学反应工程

[求助] matlab最小二乘插值拟合，跪求大神！！

function benan5
clear all;
load cc  %load experimental data

[y_row,y_col]=size(y);
[x_row,x_col]=size(x);

global b

beta0=[2 60 1 45 1.5 90 5 45 1 100 1 0 1 1 1 1 1 1 0]
y0=[0 0 0];
lb=[0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0];
ub=[+inf +inf +inf +inf +inf +inf +inf +inf +inf +inf +inf +inf +inf +inf +inf +inf +inf];

[beta,resnorm,residual,exitflag,output,lambda,jacobian]=...
lsqnonlin(@func,beta0,lb,ub,[],x,y0,y)
ci=nlparci(beta,residual,jacobian);

namex='W/Fa0 , g.h/mol';

for i=1:1:3
T1=180+273;
T2=200+273;
T3=220+273;
T4=240+273;
x1=x(1:8,i);
x2=x(9:16,i);
x3=x(17:24,i);
x4=x(25:32,i);
b=3*i-2;
x5=linspace(0,500,1000);
tspan1=[0 max(x1)];
tspan2=[0 max(x2)];
tspan3=[0 max(x3)];
tspan4=[0 max(x4)];

y1=y(1:8,b);
y2=y(9:16,b);
y3=y(17:24,b);
y4=y(25:32,b);
y5=y(1:8,b+1);
y6=y(9:16,b+1);
y7=y(17:24,b+1);
y8=y(25:32,b+1);
y9=y(1:8,b+2);
y10=y(9:16,b+2);
y11=y(17:24,b+2);
y12=y(25:32,b+2);

[t1 yy1]=ode23(@modeleqs,tspan1,y0,[],beta,T1);
[t2 yy2]=ode23(@modeleqs,tspan2,y0,[],beta,T2);
[t3 yy3]=ode23(@modeleqs,tspan3,y0,[],beta,T3);
[t4 yy4]=ode23(@modeleqs,tspan4,y0,[],beta,T4);

yc1=spline(t1,yy1(:,1),x1);
yc2=spline(t2,yy2(:,1),x2);
yc3=spline(t3,yy3(:,1),x3);
yc4=spline(t4,yy4(:,1),x4);

yc5=spline(t1,yy1(:,2),x1);
yc6=spline(t2,yy2(:,2),x2);
yc7=spline(t3,yy3(:,2),x3);
yc8=spline(t4,yy4(:,2),x4);

yc9=spline(t1,yy1(:,3),x5);
yc10=spline(t2,yy2(:,3),x5);
yc11=spline(t3,yy3(:,3),x5);
yc12=spline(t4,yy4(:,3),x5);

figure(6*i-5)
plot(x1,yc1,'r-',x1,y1,'o');hold on;
plot(x2,yc2,'b-',x2,y2,'s');hold on;
plot(x3,yc3,'m-',x3,y3,'v');hold on;
plot(x4,yc4,'g-',x4,y4,'x');hold off;
xlabel(namex)
ylabel('苯胺转化率')

figure(6*i-4)
plot(x1,yc5,'r-',x1,y5,'o');hold on;
plot(x2,yc6,'b-',x2,y6,'s');hold on;
plot(x3,yc7,'m-',x3,y7,'v');hold on;
plot(x4,yc8,'g-',x4,y8,'x');hold off;
xlabel(namex)
ylabel('二环己胺收率')

figure(6*i-3)
plot(x5,yc9,'r-',x1,y9,'o');hold on;
plot(x5,yc10,'b-',x2,y10,'s');hold on;
plot(x5,yc11,'m-',x3,y11,'v');hold on;
plot(x5,yc12,'g-',x4,y12,'x');hold off;
xlabel(namex)
ylabel('N-环己基苯胺收率')
%{
figure(6*i-2)
b=3*i-2
d1=residual(96*(i-1)+1:96*(i-1)+8,1);hold on;
d2=residual(96*(i-1)+9:96*(i-1)+16,1);hold on;
d3=residual(96*(i-1)+17:96*(i-1)+24,1);hold on;
d4=residual(96*(i-1)+25:96*(i-1)+32,1);hold off;
plot(y1,d1,'ko',y2,d2,'kv',y3,d3,'k*',y4,d4,'ks');
namex1='苯胺转化率';
namey1(1,1:length('residual values'))='Residual values';
xlabel(namex1)
ylabel(namey1)

figure(6*i-1)
d5=residual(96*(i-1)+33:96*(i-1)+40,1);hold on;
d6=residual(96*(i-1)+32+9:96*(i-1)+32+16,1);hold on;
d7=residual(96*(i-1)+32+17:96*(i-1)+32+24,1);hold on;
d8=residual(96*(i-1)+32+25:96*(i-1)+32+32,1);hold off;
plot(y5,d5,'ko',y6,d6,'kv',y7,d7,'k*',y8,d8,'ks');
namex2='二环己胺收率';
namey2(1,1:length('residual values'))='Residual values';
xlabel(namex2)
ylabel(namey2)

figure(6*i)
d9=residual(96*(i-1)+64+1:96*(i-1)+64+8,1);hold on;
d10=residual(96*(i-1)+64+9:96*(i-1)+64+16,1);hold on;
d11=residual(96*(i-1)+64+17:96*(i-1)+64+24,1);hold on;
d12=residual(96*(i-1)+64+25:96*(i-1)+64+32,1);hold off;
plot(y9,d9,'ko',y10,d10,'kv',y11,d11,'k*',y12,d12,'ks');
namex3='N-环己基苯胺收率';
namey3(1,1:length('residual values'))='Residual values';
xlabel(namex3)
ylabel(namey3)
%}
end

%output of parameters needed
fprintf('\n estimated parameters by lsqnonlin():\n')
fprintf('\t k1=%10.6f±%10.4f\n',beta(1),ci(1,2)-beta(1))
fprintf('\t k2=%10.6f±%10.4f\n',beta(2),ci(2,2)-beta(2))
fprintf('\t k3=%10.6f±%10.4f\n',beta(3),ci(3,2)-beta(3))
fprintf('\t k4=%10.6f±%10.4f\n',beta(4),ci(4,2)-beta(4))
fprintf('\t k5=%10.6f±%10.4f\n',beta(5),ci(5,2)-beta(5))
fprintf('\t k6=%10.4f±%10.4f\n',beta(6),ci(6,2)-beta(6))
fprintf('\t k7=%10.4f±%10.4f\n',beta(7),ci(7,2)-beta(7))
fprintf('\t k8=%10.4f±%10.4f\n',beta(8),ci(8,2)-beta(8))
fprintf('\t k9=%10.4f±%10.4f\n',beta(9),ci(9,2)-beta(9))
fprintf('\t k10=%10.4f±%10.4f\n',beta(10),ci(10,2)-beta(10))
fprintf('\t k11=%10.6f±%10.4f\n',beta(11),ci(11,2)-beta(11))
fprintf('\t k12=%10.6f±%10.4f\n',beta(12),ci(12,2)-beta(12))
fprintf('\t k13=%10.6f±%10.4f\n',beta(13),ci(13,2)-beta(13))
fprintf('\t k14=%10.6f±%10.4f\n',beta(14),ci(14,2)-beta(14))
fprintf('\t k15=%10.6f±%10.4f\n',beta(15),ci(15,2)-beta(15))
fprintf('\t k16=%10.4f±%10.4f\n',beta(16),ci(16,2)-beta(16))
fprintf('\t k17=%10.4f±%10.4f\n',beta(17),ci(17,2)-beta(17))
fprintf('\t k18=%10.4f±%10.4f\n',beta(18),ci(18,2)-beta(18))
fprintf('\t k19=%10.4f±%10.4f\n',beta(19),ci(19,2)-beta(19))
fprintf(' the sum of the residual squares is:%.4e\n\n',sum((residual).^2))

%=======================================================
function f=func(beta,x,y0,y)
%define objective function
global b
for i=1:1:3
T1=180+273;
T2=200+273;
T3=220+273;
T4=240+273;
x1=x(1:8,i);
x2=x(9:16,i);
x3=x(17:24,i);
x4=x(25:32,i);
b=3*i-2;

tspan1=[0 max(x1)];
tspan2=[0 max(x2)];
tspan3=[0 max(x3)];
tspan4=[0 max(x4)];

[t1 yy1]=ode23(@modeleqs,tspan1,y0,[],beta,T1);
[t2 yy2]=ode23(@modeleqs,tspan2,y0,[],beta,T2);
[t3 yy3]=ode23(@modeleqs,tspan3,y0,[],beta,T3);
[t4 yy4]=ode23(@modeleqs,tspan4,y0,[],beta,T4);

y1=y(1:8,b);
y2=y(9:16,b);
y3=y(17:24,b);
y4=y(25:32,b);

y5=y(1:8,b+1);
y6=y(9:16,b+1);
y7=y(17:24,b+1);
y8=y(25:32,b+1);

y9=y(1:8,b+2);
y10=y(9:16,b+2);
y11=y(17:24,b+2);
y12=y(25:32,b+2);

j=1;
yc1=spline(t1,yy1(:,j),x1);
yc2=spline(t2,yy2(:,j),x2);
yc3=spline(t3,yy3(:,j),x3);
yc4=spline(t4,yy4(:,j),x4);
%苯胺转化率拟合
j=2;
yc5=spline(t1,yy1(:,j),x1);
yc6=spline(t2,yy2(:,j),x2);
yc7=spline(t3,yy3(:,j),x3);
yc8=spline(t4,yy4(:,j),x4);
   %二环己胺收率拟合
j=3;
yc9=spline(t1,yy1(:,j),x1);
yc10=spline(t2,yy2(:,j),x2);
yc11=spline(t3,yy3(:,j),x3);
yc12=spline(t4,yy4(:,j),x4);
%N-环己基苯胺收率拟合

f11=y1-yc1;
f12=y2-yc2;
f13=y3-yc3;
f14=y4-yc4;
f21=y5-yc5;
f22=y6-yc6;
f23=y7-yc7;
f24=y8-yc8;
f31=y9-yc9;
f32=y10-yc10;
f33=y11-yc11;
f34=y12-yc12;

g(i,

=[f11' f12' f13' f14' f21' f22' f23' f24' f31' f32' f33' f34'];
end
f=[g(1,

';g(2,

';g(3,

'];
%===================================================================================
function dydt=modeleqs(t,yt,beta,T)
global b
P=0.3;
N=5*(b+2)/3+5;%N为氢胺比
X1=yt(1);
Y4=yt(2);
Y6=yt(3);
X1
Y4
Y6
P1=P.*(1-X1)./(N+1-3.*X1+3.*Y6);
P2=P.*(N-3.*X1+3.*Y6)./(N+1-3.*X1+3.*Y6);
P3=P.*(X1-Y4-2.*Y6)./(N+1-3.*X1+3.*Y6);
P4=P.*Y4./(2.*N+2-6.*X1+6.*Y6);
P5=P.*(Y4+2.*Y6)./(2.*N+2-6.*X1+6.*Y6);
P6=P.*Y6./(N+1-3.*X1+3.*Y6);
P1
P3
P4
P5
P6
K(1)=beta(1).*10.^6.*exp(-beta(2).*10.^3./(8.314.*T)); %k(1)、k(2)第一步反应的指前因子、表观活化能
K(2)=beta(3).*10.^7.*exp(-beta(4).*10.^3./(8.314.*T)); %k(3)、k(4)第二步反应的指前因子、表观活化能
K(3)=beta(5).*10.^10.*exp(-beta(6).*10.^3./(8.314.*T)); %k(5)、k(6)第三步反应的指前因子、表观活化能
K(4)=beta(7).*10.^11.*exp(-beta(8).*10.^3./(8.314.*T)); %k(7)、k(8)第四步反应的指前因子、表观活化能
K(5)=beta(9).*10.^12.*exp(-beta(10).*10.^3./(8.314.*T)); %k(9)、k(10)第五步反应的指前因子、表观活化能
r1=K(1).*P1.^beta(11).*P2.^beta(12);
r2=K(2).*P3.^beta(13)-K(3).*P4.^beta(14).*P5.^beta(15);
r3=K(4).*P1.^beta(16).*P3.^beta(17);
r4=K(5).*P6.^beta(18).*P2.^beta(19);
R1=r1+r3
R4=r2+r4
R6=r3-r4
dydt=[r1+r3;r2+r4;r3-r4];
求大神检查一下程序是否有问题。我现在拟合不出来结果，急！实验数据在附件里

回复此楼