kineticsest4.m

matlab在化学工程中的应用实例程序,对大家研究学习有指导作用。

function KineticsEst4
% 一氧化氮催化还原反应的动力学参数估计
%
%   Author: HUANG Huajiang
%   Copyright 2003 UNILAB Research Center, 
%   East China University of Science and Technology, Shanghai, PRC
%   $Revision: 1.0 $  $Date: 2003/04/16 $
%
%  [Ref]：Englezos & Kalogerakis. Applied Parameter Estimation for Chemical
%   Engineers. Marcel Dekker, 2000（16.1.3）

clear all
clc

% Load experimental data
load ChemKineticsData
PH2 = ExpData(:,1);
PNO = ExpData(:,2);
r = ExpData(:,3);

% 用多变量线性回归方法估计动力学参数
R = sqrt(PH2.*PNO./r);
y = R;
X = [ones(size(y))  PH2  PNO];
[b,bint] = regress(y,X,0.05);    % 或b = X\y
denom = 1/b(1);
KH2 = b(2)*denom;
KNO = b(3)*denom;
k = denom^2/(KH2*KNO);

% 用lsqnonlin()--求解非线性最小二乘法(非线性数据拟合)问题
beta0 = [k  KH2  KNO];
lb = [0  0  0];
ub = [+inf  +inf  +inf]
[beta,resnorm,resid,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFun,beta0,lb,ub,[],PH2,PNO,r);
ci = nlparci(beta,resid,jacobian);

% 模型适定性判别
Ne = length(r);
Np = length(beta);
[rho2, F] = rho2_F(k, r, resnorm, Ne, Np);

% 残差关于拟合值的残差图
rc = RateEqs(beta,PH2,PNO);
plot(rc,resid,'*')
xlabel('反应速率拟合值, mol/(min g催化剂)')
ylabel('残差R, mol/(min g催化剂)')
refline(0,0)

% 参数辨识结果
fprintf('\n\nEstimated Parameters:\n')
fprintf('\tk = %.4f ± %.4f\n',beta(1),ci(1,2)-beta(1))
fprintf('\tKH2 = %.2f ± %.2f\n',beta(2),ci(2,2)-beta(2))
fprintf('\tKNO = %.2f ± %.2f\n',beta(3),ci(3,2)-beta(3))
fprintf('  决定性指标ρ^2: %.3f\n',rho2)
fprintf('  F比: %.3f\n\n',F)


% ------------------------------------------------------------------
function f = ObjFun(beta,PH2,PNO,r)
rc = RateEqs(beta,PH2,PNO);
f = r - rc;

% ------------------------------------------------------------------
function rc = RateEqs(beta,PH2,PNO)    % Rate equation 
rc = beta(1)*beta(2)*beta(3)*PH2.*PNO./(1+beta(3)*PNO+beta(2)*PH2).^2;