升温速率为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))