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;