24小时热门版块排行榜

>论坛更新日志 (3640)
>文献求助 (338)
>虫友互识 (291)
>导师招生 (224)
>硕博家园 (128)
>招聘信息布告栏 (115)
>考博 (104)
>休闲灌水 (83)
>论文道贺祈福 (76)
>博后之家 (71)
>考研 (56)
>论文投稿 (54)
>教师之家 (51)
>基金申请 (43)
>公派出国 (43)
>找工作 (31)

返回列表

ZHM10

新虫 (小有名气)

应助: 0 (幼儿园)
金币: 3506.9
帖子: 240
在线: 274.7小时
虫号: 3796100
注册: 2015-04-08
专业: 应用有机化学

[求助] Matlab拟合动力学参数问题，求指导！已有1人参与

小弟现在要用matlab拟合出一个复杂反应体系的动力学参数，先用ode45求解微分方程组再用最小二乘拟合，参照文献编了一个程序，但运行老出一些问题，快急死了…
求大神们帮忙看看，有什么问题…还有个问题，就是我看文献上需要先给动力学参数一个初值，那这个初值该怎么选择啊？
下面是程序、动力学模型以及实验数据，先谢谢各位了！
clear all; %%清理内存
clc;    %%清屏
global C_A_EXP C_B_EXP C_C1_EXP C_C2_EXP C_C3_EXP C_D_EXP C_E_EXP %%各个物质浓度实验值
global C_A_CAL C_B_CAL C_C1_CAL C_C2_CAL C_C3_CAL C_D_CAL C_E_CAL %%各个物质浓度计算值
global R tspan C00  %%理想气体常数、积分区间、初值
R=8.314;
tspan={[0 5];[0 10];[0 20];[0 30];[0 45];[0 60];[0 90];[0 120];[0 180]};
%% 读取数据
display '读取50 下数据...';data_expx_50=xlsread('data_expx.xlsx','50','B5:H13');
display '读取60 下数据...';data_expx_60=xlsread('data_expx.xlsx','60','B5: H13');
display '读取70 下数据...';data_expx_70=xlsread('data_expx.xlsx','70','B5: H13');
display '读取80 下数据...';data_expx_80=xlsread('data_expx.xlsx','80','B5: H13');
display '读取初值...'
display '读取50 下初值...';data_expx_50_0=xlsread('data_expx.xlsx','50','B4:H4');
display '读取60 下初值...';data_expx_60_0=xlsread('data_expx.xlsx','60','B4:H4');
display '读取70 下初值...';data_expx_70_0=xlsread('data_expx.xlsx','70','B4:H4');
display '读取80 下初值...';data_expx_80_0=xlsread('data_expx.xlsx','80','B4:H4');
display '读取数据结束...'
%% 浓度实验值数据
C_A_EXP=[ data_expx_50(:,1)...data_expx_60(:,1)...data_expx_70(:,1)...
data_expx_80(:,1)];
C_B_EXP=[ data_expx_50(:,2)...data_expx_60(:,2)...data_expx_70(:,2)...
data_expx_80(:,2)];
C_C1_EXP=[ data_expx_50(:,3)...data_expx_60(:,3)...data_expx_70(:,3)...
data_expx_80(:,3)];
C_C2_EXP=[ data_expx_50(:,4)...data_expx_60(:,4)...data_expx_70(:,4)...
data_expx_80(:,4)];
C_C3_EXP=[ data_expx_50(:,5)...data_expx_60(:,5)...data_expx_70(:,5)...
data_expx_80(:,5)];
C_D_EXP=[ data_expx_50(:,6)...data_expx_60(:,6)...data_expx_70(:,6)...
data_expx_80(:,6)];
C_E_EXP=[ data_expx_50(:,7)...data_expx_60(:,7)...data_expx_70(:,7)...
data_expx_80(:,7)];
%% 初值
C_A_EXP_0=[data_expx_50_0(:,1)...data_expx_60_0(:,1)...data_expx_70_0(:,1)...
data_expx_80_0(:,1)];
C_B_EXP_0=[data_expx_50_0(:,2)...data_expx_60_0(:,2)...data_expx_70_0(:,2)...
data_expx_80_0(:,2)];
C_C1_EXP_0=[data_expx_50_0(:,3)...data_expx_60_0(:,3)...data_expx_70_0(:,3)...
data_expx_80_0(:,3)];
C_C2_EXP_0=[data_expx_50_0(:,4)...data_expx_60_0(:,4)...data_expx_70_0(:,4)...
data_expx_80_0(:,4)];
C_C3_EXP_0=[data_expx_50_0(:,5)...data_expx_60_0(:,5)...data_expx_70_0(:,5)...
data_expx_80_0(:,5)];
C_D_EXP_0=[data_expx_50_0(:,6)...data_expx_60_0(:,6)...data_expx_70_0(:,6)...
data_expx_80_0(:,6)];
C_E_EXP_0=[data_expx_50_0(:,7)...data_expx_60_0(:,7)...data_expx_70_0(:,7)...
data_expx_80_0(:,7)];
C00=[C_A_EXP_0 C_B_EXP_0 C_C1_EXP_0 C_C2_EXP_0 C_C3_EXP_0 C_D_EXP_0 C_E_EXP_0];
p0=[2.0e+5 1.0e+004 1.0...
2.0e+5 1.0e+004 1.0...
2.0e+5 1.0e+004 1.0...
2.0e+5 1.0e+004 1.0...
2.0e+5 1.0e+004 1.0...
2.0e+5 1.0e+004 1.0...
2.0e+5 1.0e+004 1.0...
2.0e+5 1.0e+004 1.0];
options=optimset('tolfun',0.1,'TolX',1e4,'MaxFunEvals',9e25,'display','iter');
[p,exitflag]=fminsearch(@minfun_,p0,options);
%% 显示结果
disp 'p1=';p(1)
disp 'p2=';p(2)
disp 'p3=';p(3)
disp 'p4=';p(4)
disp 'p5=';p(5)
disp 'p6=';p(6)
disp 'p7=';p(7)
disp 'p8=';p(8)
disp 'p9=';p(9)
disp 'p10=';p(10)
disp 'p11=';p(11)
disp 'p12=';p(12)
disp 'p13=';p(13)
disp 'p14=';p(14)
disp 'p15=';p(15)
disp 'p16=';p(16)
disp 'p17=';p(17)
disp 'p18=';p(18)
disp 'p19=';p(19)
disp 'p20=';p(20)
disp 'p21=';p(21)
disp 'p22=';p(22)
disp 'p23=';p(23)
disp 'p24=';p(24)

function dCdt=myfun_(t,C,p,m)
global R
%%
%% 计算各个物质反应速率
%%
if m<=1
TC=323.15;
elseif m>1 && m<=2
TC=433.15;
elseif m>2 && m<=3
TC=443.15;
elseif m>3
TC=453.15;
end
R1=p(1).*exp(-p(2)./(R*TC)).*(C(1).^p(3));
R2=p(4).*exp(-p(5)./(R*TC)).*(C(2).^p(6));
R3=p(7).*exp(-p(8)./(R*TC)).*(C(2).^p(9));
R4=p(10).*exp(-p(11)./(R*TC)).*(C(2).^p(12));
R5=p(13).*exp(-p(14)./(R*TC)).*(C(3).^p(15));
R6=p(16).*exp(-p(17)./(R*TC)).*(C(4).^p(18));
R7=p(19).*exp(-p(20)./(R*TC)).*(C(5).^p(21));
R8=p(22).*exp(-p(23)./(R*TC)).*(C(6).^p(24));
dCdt=[-R1;R1-R2-R3-R4;R2-R5;R3-R6;R4-R7;R5+R6+R7-R8;R8];

function min=minfun_(p)
global R tspan C00
for i=1:1:9
for m=1:1:7
C0=[C00(m) C00(m+7) C00(m+14) C00(m+21)];
[t,C]=ode45(@myfun,tspan{i},C0,[],p,m);
C_A_CAL(i,m)=C(end,1);
C_B_CAL(i,m)=C(end,2);
C_C1_CAL(i,m)=C(end,3);
C_C2_CAL(i,m)=C(end,4);
C_C3_CAL(i,m)=C(end,5);
C_D_CAL(i,m)=C(end,6);
C_E_CAL(i,m)=C(end,7);
end
end
min = sum(sum((C_A_EXP-C_A_CAL).^2+(C_B_EXP-C_B_CAL).^2+…
(C_C1_EXP-C_C1_CAL).^2+(C_C2_ EXP-C_C2_CAL).^2+…
(C_C3_EXP-C_C3_CAL).^2+(C_D_EXP-C_D_CAL).^2+…
(C_E_ EXP-C_E_CAL).^2));

动力学模型1.jpg