kineticsest2.m

来自「《实用化工计算机模拟：MATLAB在化学工程中的应用 》这本书光盘里的程序~」· M 代码 · 共 52 行

function  KineticsEst2
% 动力学参数估计
%
%   Author: HUANG Huajiang
%   Copyright 2003 UNILAB Research Center, 
%   East China University of Science and Technology, Shanghai, PRC
%   $Revision: 1.0 $  $Date: 2003/03/23 $
%
%   [Ref] H. Scott Fogler, Elements of Chemical Reaction Engineering, Prentice Hall PTR, 
%   Upper Saddle River, New Jersey, 1999

clear all
clc

rA = [1.04 3.13 5.21 3.82 4.19 2.391 3.867 2.199 0.75];  
Pe = [1 1 1 3 5 0.5 0.5 0.5 0.5];  
Ph = [1 3 5 3 3 3 5 3 1];

% 线性拟合
RA = Pe.*Ph./rA;
y = RA';
X = [ones(size(y))  Pe'];
b = X\y;
k = 1/b(1);
KE = b(2)*k;

% 非线性拟合
beta0 = [k KE]
[beta,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc,beta0,[],[],[],rA,Pe,Ph)           
ci = nlparci(beta,residual,jacobian)

% 残差关于拟合值的残差图
rAc = RateA(beta,Pe,Ph)
plot(rAc,residual,'*')
xlabel('反应速率（-rA）拟合值, mol/(kg cat. s)')
ylabel('残差R, mol/(kg cat. s)')
refline(0,0)

% 参数辨识结果
fprintf('Estimated Parameters:\n')
fprintf('\tk = %.3f ± %.3f\n',beta(1),ci(1,2)-beta(1))
fprintf('\tKE = %.3f ± %.3f\n',beta(2),ci(2,2)-beta(2))
fprintf('\tThe sum of the squares is: %.3f',resnorm)

% ------------------------------------------------------------------
function f = ObjFunc(beta,rA,Pe,Ph)
f = rA - RateA(beta,Pe,Ph);

% ------------------------------------------------------------------
function rA = RateA(beta,Pe,Ph)
rA = beta(1)*Pe.*Ph./(1+beta(2)*Pe);