function KineticsEst1_int
% 動(dòng)力學(xué)參數(shù)辨識(shí): 用積分法進(jìn)行反應(yīng)速率分析得到速率常數(shù)k和反應(yīng)級(jí)數(shù)n
% Analysis of kinetic rate data by using the integral method
%
%   Author: HUANG Huajiang
%   Copyright 2003 UNILAB Research Center,
%   East China University of Science and Technology, Shanghai, PRC
%   $Revision: 1.0 $  $Date: 2003/07/27 $
%
% 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

global CAm
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)

% 擬合效果圖(實(shí)驗(yàn)與擬合的比較)
[t4plot CA4plot] = ode45(@KineticsEqs,[tspan(1)  tspan(end)],CA0,[],beta);
plot(t,CAm,'bo',t4plot,CA4plot,'k-')
legend('Exp','Model')
xlabel('時(shí)間t, s')
ylabel('濃度C_A, mol/L')

% 殘差關(guān)于擬合值的殘差圖
[t CAc] = ode45(@KineticsEqs,tspan,CA0,[],beta);
figure
plot(CAc,resid,'*')
xlabel('濃度擬合值（mol/L）')
ylabel('殘差R (mol/L)')
refline(0,0)

% 參數(shù)辨識(shí)結(jié)果
fprintf('Estimated Parameters:\n')
fprintf('\tk = %.4f ± %.4f\n',beta(1),ci(1,2)-beta(1))
fprintf('\tn = %.2f ± %.2f\n',beta(2),ci(2,2)-beta(2))

% ------------------------------------------------------------------
function f = OptObjFunc(beta,tspan,CA0,CAm)
[t CAc] = ode45(@KineticsEqs,tspan,CA0,[],beta);
f = CAc - CAm;

% ------------------------------------------------------------------
function dCAdt = KineticsEqs(t,CA,beta)
dCAdt = -beta(1)*CA^beta(2);            % k = beta(1), n = beta(2)