姓名：方健学号：652081701073 专业：化学工程班级：化工1503
问题：在一个等温间歇反应器中进行复杂反应，描述该反应系统的七个微分方程共有5个参数k1，k2，k3，k4和k5，初试状态x（0）=[0.1883 0.2507 0.0467 0.0899
T
0.1804 0.1394 0.1046],试根据下表的数据，利用MATLAB最优化方法函数估计出k参数，表中的数据存于文件KineticsData1.m中 % 动力学数据:
% t y1 y2 y3 y4 ExpData= ...
[ 0 0.1883 0.0899 0.1804 0.1394 0.0100 0.2047 0.0866 0.1729 0.1297 0.0200 0.2181 0.0856 0.1680 0.1205 0.0300 0.2291 0.0863 0.1647 0.1123 0.0400 0.2382 0.0878 0.1623 0.1053 0.0500 0.2459 0.0899 0.1604 0.0995 0.0600 0.2523 0.0921 0.1588 0.0948 0.0700 0.2576 0.0945 0.1574 0.0911 0.0800 0.2622 0.0968 0.1561 0.0882 0.0900 0.2660 0.0989 0.1548 0.0859 0.1000 0.2692 0.1010 0.1537 0.0842 0.1100 0.2719 0.1028 0.1525 0.0830 0.1200 0.2742 0.1045 0.1515 0.0821 0.1300 0.2761 0.1060 0.1505 0.0814 0.1400 0.2777 0.1074 0.1495 0.0810 0.1500 0.2790 0.1086 0.1487 0.0807 0.1600 0.2801 0.1096 0.1479 0.0805 0.1700 0.2811 0.1106 0.1471 0.0803 0.1800 0.2819 0.1114 0.1465 0.0803 0.1900 0.2825 0.1121 0.1458 0.0803 0.2000 0.2830 0.1127 0.1453 0.0803 ]
x = lsqnonlin(fun,x0)
x = lsqnonlin(fun,x0,lb,ub)
x = lsqnonlin(fun,x0,lb,ub,options) x = lsqnonlin(problem)
[x,resnorm,residual,exitflag,output,lambda,jacobian] = lsqnonlin(...)
x = lsqcurvefit(fun,x0,xdata,ydata)
x = lsqcurvefit(fun,x0,xdata,ydata,lb,ub)
x = lsqcurvefit(fun,x0,xdata,ydata,lb,ub,options) x = lsqcurvefit(problem)
[x,resnorm,residual,exitflag,output,lambda,jacobian] = lsqcurvefit(...)
根据问题首先进行方程参数的估计编写KineticsEst5,m文件如下：
function KineticsEst5
% 动力学ODE方程模型的参数估计 %
% The variables y here are y(1)=x1, y(2)=x4, y(3)=x5,y(4)=x6 . clear all clc
k0 = [0.5 0.5 0.5 0.5 0.5]; % 参数初值 lb = [0 0 0 0 0]; % 参数下限 ub = [+inf +inf +inf +inf +inf]; % 参数上限
x0 = [0.1883 0.2507 0.0467 0.0899 0.1804 0.1394 0.1046]; KineticsData1;
yexp = ExpData(:,2:5); % yexp: 实验数据[x1 x4 x5 x6]
% 使用函数fmincon()进行参数估计
[k,fval,flag] = fmincon(@ObjFunc4Fmincon,k0,[],[],[],[],lb,ub,[],[],x0,yexp); fprintf('\\n使用函数fmincon()估计得到的参数值为:\\n') fprintf('\\tk1 = %.4f\\n',k(1)) fprintf('\\tk2 = %.4f\\n',k(2)) fprintf('\\tk3 = %.4f\\n',k(3)) fprintf('\\tk4 = %.4f\\n',k(4)) fprintf('\\tk5 = %.4f\\n',k(5))
fprintf(' The sum of the squares is: %.1e\\n\\n',fval) k_fmincon = k;
% 使用函数lsqnonlin()进行参数估计
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ... lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[],x0,yexp); k
ci = nlparci(k,residual,jacobian)
fprintf('\\n\\n使用函数lsqnonlin()估计得到的参数值为:\\n') output
% 以函数fmincon()估计得到的结果为初值，使用函数lsqnonlin()进行参数估计 k0 = k_fmincon;
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ... lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[],x0,yexp); k
ci = nlparci(k,residual,jacobian)
fprintf('\\n\\n以fmincon()的结果为初值，使用函数lsqnonlin()估计得到的参数值为:\\n') output
% ------------------------------------------------------------------ function f = ObjFunc4Fmincon(k,x0,yexp)
tspan = [0.00 : 0.01 : 0.20];
[t x] = ode45(@KineticEqs,tspan,x0,[],k); y(:,1) = x(:,1); y(:,2:4) = x(:,4:6);
f = sum((y(:,1)-yexp(:,1)).^2) + sum((y(:,2)-yexp(:,2)).^2) ... + sum((y(:,3)-yexp(:,3)).^2) + sum((y(:,4)-yexp(:,4)).^2);
% ------------------------------------------------------------------ function f = ObjFunc4LNL(k,x0,yexp) tspan = [0.00 : 0.01 : 0.20];
[t x] = ode45(@KineticEqs,tspan,x0,[],k); y(:,1) = x(:,1); y(:,2:4) = x(:,4:6); f1 = y(:,1) - yexp(:,1); f2 = y(:,2) - yexp(:,2); f3 = y(:,3) - yexp(:,3); f4 = y(:,4) - yexp(:,4); f = [f1; f2; f3; f4];
% ------------------------------------------------------------------ functiondxdt = KineticEqs(t,x,k) q = 8.75 + k(5); dxdt = ...
[ ( k(5)-q*x(1)- k(1)*x(1)*x(2)-k(4)*x(1)*x(6)*sqrt(0.9) ) ( 7.0-q*x(2) - k(1)*x(1)*x(2)-2*k(2)*x(2)*x(3) ) ( 1.75 -q*x(3) - k(2)*x(2)*x(3) )
( -q*x(4) + 2*k(1)*x(1)*x(2)-k(3)*x(4)*x(5) ) ( -q*x(5) + 3*k(2)*x(2)*x(3)-k(3)*x(4)*x(5) )
( -q*x(6) + 2*k(3)*x(4)*x(5)-k(4)*x(1)*x(6)*sqrt(0.9) ) ( -q*x(7) + 2*k(4)*x(1)*x(6)*sqrt(0.9) ) ];
在MATLAB中运行结果如下： ExpData =
0 0.1883 0.0899 0.1804 0.0100 0.2047 0.0866 0.1729 0.0200 0.2181 0.0856 0.1680 0.0300 0.2291 0.0863 0.1647 0.0400 0.2382 0.0878 0.1623 0.0500 0.2459 0.0899 0.1604 0.0600 0.2523 0.0921 0.1588 0.0700 0.2576 0.0945 0.1574 0.0800 0.2622 0.0968 0.1561 0.0900 0.2660 0.0989 0.1548 0.1394 0.1297 0.1205 0.1123 0.1053 0.0995 0.0948 0.0911 0.0882 0.0859
0.1000 0.2692 0.1010 0.1537 0.0842 0.1100 0.2719 0.1028 0.1525 0.0830 0.1200 0.2742 0.1045 0.1515 0.0821 0.1300 0.2761 0.1060 0.1505 0.0814 0.1400 0.2777 0.1074 0.1495 0.0810 0.1500 0.2790 0.1086 0.1487 0.0807 0.1600 0.2801 0.1096 0.1479 0.0805 0.1700 0.2811 0.1106 0.1471 0.0803 0.1800 0.2819 0.1114 0.1465 0.0803 0.1900 0.2825 0.1121 0.1458 0.0803
0.2000 0.2830 0.1127 0.1453 0.0803 使用函数fmincon()估计得到的参数值为: k1 = 17.5835 k2 = 72.7641 k3 = 51.1175 k4 = 22.9402 k5 = 5.9949
The sum of the squares is: 6.4e-07 k =
17.6085 73.0645 51.3262 23.0250 6.0012 ci =
17.5954 17.6215 72.9804 73.1485 51.2783 51.3742 22.9682 23.0817 5.9989 6.0036
使用函数lsqnonlin()估计得到的参数值为: output =
firstorderopt: 2.6911e-07 iterations: 9 funcCount: 60 cgiterations: 0
algorithm: 'trust-region-reflective' message: [1x425 char] k =
17.6063 73.0504 51.3184 23.0156 6.0009 ci =
17.5933 17.6194 72.9662 73.1345 51.2704 51.3663 22.9588 23.0724 5.9985 6.0032
以fmincon()的结果为初值，使用函数lsqnonlin()估计得到的参数值为:
output =
firstorderopt: 8.3868e-07 iterations: 1 funcCount: 12 cgiterations: 0
algorithm: 'trust-region-reflective' message: [1x425 char] 在微分法的情况下： function KineticsEst1_int
% 动力学参数辨识: 用积分法进行反应速率分析得到速率常数k和反应级数n
% Analysis of kinetic rate data by using the integral method %
% Reaction of the type -- rate = kCA^order % order - reaction order % rate -- reaction rate vector
% CA -- concentration vector for reactant A % T -- vector of reaction time % N -- number of data points % k- reacion rate constant clear all clc
globalCAm
t = [0 20 40 60 120 180 300]; CAm = [10 8 6 5 3 2 1]';
% 非线性拟合
beta0 = [0.0053 1.39];
tspan = [0 20 40 60 120 180 300]; CA0 = 10;
[beta,resnorm,resid,exitflag,output,lambda,jacobian] = ... lsqnonlin(@OptObjFunc,beta0,[],[],[],tspan,CA0,CAm) ci = nlparci(beta,resid,jacobian)
% 拟合效果图(实验与拟合的比较)
[t4plot CA4plot] = ode45(@KineticsEqs,[tspan(1) tspan(end)],CA0,[],beta); plot(t,CAm,'bo',t4plot,CA4plot,'k-') legend('Exp','Model') xlabel('时间t, s')
ylabel('浓度C_A, mol/L')
% 残差关于拟合值的残差图
[tCAc] = ode45(@KineticsEqs,tspan,CA0,[],beta); figure