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.00160.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 &nbp; 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))