bdt模型 matlab,matlab程序求助-连串反应动力学参数拟合

升温速率为20℃/min

function g=g(x)

g=(x-159.3100)/20    %通过温度数据求出对应时间

%A数据为随意补充，通过A和C可求B

%建立数据m文件：

%时间t            温度T               C                      B                      A

%t   y1=x1   y2=x2   y3=x3     y4=x4

Kinetics =...

[   0.0205      159.7200        0.0016  0.00015   1.0000

0.0965       161.2400    0.0002    0.0002    0.9996

0.1725      162.7600    0.0008    0.0003    0.9989

0.2480      164.2700    0.0021    0.0005    0.9974

0.3240       165.7900    0.0046    0.0020    0.9934

0.3995        167.3000    0.0084    0.0030    0.9886

0.4755        168.8200    0.0134    0.0040    0.9826

0.5515       170.3400    0.0197    0.0050    0.9753

0.6270        171.8500    0.0275    0.0060    0.9665

0.7030        173.3700    0.0373    0.0070    0.9557

0.7790         174.8900    0.0498    0.0080    0.9422

0.8545         176.4000    0.0638    0.0090    0.9272

0.9305          177.9200    0.0799    0.0100    0.9101

1.0065          179.4400    0.0981    0.0200    0.8819

1.0820          180.9500    0.1191    0.0300    0.8509

1.1580         182.4700    0.1426    0.0310    0.8264

1.2335          183.9800    0.1688    0.0340    0.7972

1.3095          185.5000    0.1974    0.0350    0.7676

1.3855          187.0200    0.2284    0.0360    0.7356

1.4610          188.5300    0.2619    0.0370    0.7011

1.5370          190.0500    0.2978    0.0400    0.6622

1.6130          191.5700    0.3360    0.0600    0.6040

1.6885          193.0800    0.3764    0.0640    0.5596

1.7645          194.6000    0.4189    0.0670    0.5141

1.8405           196.1200    0.4633    0.0690    0.4677

1.9160           197.6300    0.5097    0.0700    0.4203

1.9920           199.1500    0.5576    0.1000    0.3424

2.0675          200.6600    0.6061    0.0150    0.3789

2.1435           202.1800    0.6540    0.0200    0.3260

2.2195           203.7000    0.7008    0.0230    0.2762

2.2950          205.2100    0.7458    0.0500    0.2042

2.3710            206.7300    0.7882    0.1000    0.1118

2.4470             208.2500    0.8270    0.0900    0.0830

2.5225              209.7600    0.8624    0.0800    0.0576

2.5985              211.2800    0.8931    0.0700    0.0369

2.6740               212.7900    0.9192    0.0600    0.0208

2.7500               214.3100    0.9404    0.0500    0.0096

2.8260               215.8300    0.9590    0.0400    0.0010

2.9015               217.3400    0.9744    0.0200    0.0056

2.9775               218.8600    0.9864    0.0100    0.0036

3.0535               220.3800    0.9950    0.0050         0

3.1290              221.8900    0.9988         0    0.0012

]

%      m文件  ：

function f = f(k,tspan,x0,yexp)           % 目标函数

[t Xsim] = ode45(@g,tspan,x0,[],k);

ysim(:,1) = Xsim(2:end,1);

ysim(:,2) = Xsim(2:end,2);

ysim(:,3) = Xsim(2:end,3);

ysim(:,4) = Xsim(2:end,4);

f = [ysim(:,1)-yexp(:,1); ysim(:,2)-yexp(:,2);  ysim(:,3)-yexp(:,3);ysim(:,4)-yexp(:,4)];

function dCdt =g(t,C,k)              % ODE模型方程

dCTdt =20;

dCCdt=k(1)*exp(-k(2)./(8.314*C(1))).*C(4).^k(3)-k(1)*exp(-k(2)./(8.314*C(1))).*C(4).^k(3)+k(4)*exp(-k(5)./(8.314*C(1))).*C(4).^k(6).*(1+k(7).*C(2));

dCBdt=k(1)*exp(-k(2)./(8.314*C(1))).*C(4).^k(3)-k(4)*exp(-k(5)./(8.314*C(1))).*C(4).^k(6).*(1+k(7).*C(2));

dCAdt=-k(1)*exp(-k(2)./(8.314*C(1))).*C(4).^k(3);

dCdt = [dCTdt; dCCdt; dCBdt;dCAdt];

%       command    ：

clear

tspan =[0 0.0205 0.0965 0.1725 0.2480 0.3240 0.3995 0.4755 0.5515 0.6270 0.7030 0.7790 0.8545 0.9305 1.0065 1.0820 1.1580 1.2335 1.3095 1.3855 1.4610 1.5370 1.6130 1.6885 1.7645 1.8405 1.9160 1.9920 2.0675 2.1435 2.2195 2.2950 2.3710 2.4470 2.5225 2.5985 2.6740 2.7500 2.8260 2.9015 2.9775 3.0535 3.1290];

x0=[159.7200 0.0016 0.00015 1.0000];

k0=[1 1 1 1 1 1 1];

lb=[0 0 0 0 0 0 0];               % 简化把参数的范围均假设为0:10

ub = [10 10 10 10 10];

qqData1;

yexp = qq(:,2:5);

[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...

lsqnonlin(@f,k0,lb,ub,[],tspan,x0,yexp);

ci = nlparci(k,residual,jacobian);

fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')

fprintf('\tk1 = %.4f ± %.4f\n',k(1),ci(1,2)-k(1))

fprintf('\tk2 = %.4f ± %.4f\n',k(2),ci(2,2)-k(2))

fprintf('\tk3 = %.4f ± %.4f\n',k(3),ci(3,2)-k(3))

fprintf('\tk4 = %.4f ± %.4f\n',k(4),ci(4,2)-k(4))

fprintf('\tk5 = %.4f ± %.4f\n',k(5),ci(5,2)-k(5))

fprintf('\tk6 = %.4f ± %.4f\n',k(6),ci(6,2)-k(6))

fprintf('\tk7= %.4f ± %.4f\n',k(7),ci(7,2)-k(7))