function KineticsEst6
% 動力學ODE方程模型的參數估計
% 
%   Author: HUANG Huajiang
%   Copyright 2003 UNILAB Research Center, 
%   East China University of Science and Technology, Shanghai, PRC
%   $Revision: 1.0 $  $Date: 2003/05/08 $
%
%  [Ref]: Belohlav Z., Zamostny P., Kluson P. and Volf J. Application of 
%   Random-Search Algorithm for Regression Analysis of Catalytic 
%   Hydrogenations. The Canadian Journal of Chemical Engineering, V.75, 
%   1997:738-739.

clear all
clc
tspan = [0 15  30  45  60  75  90  120  180  240  320  360  380  400];
x0 = [1; 0; 0];
k0 = [0.02  0.005  0.025  1.5  1.25];   
lb = [0  0  0  0  0];
ub = [1  1  1  10  10];

KineticsData2;
yexp = Kinetics(:,2:4);

% 使用函數lsqnonlin()進行參數估計
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc,k0,lb,ub,[],tspan,x0,yexp);       
ci = nlparci(k,residual,jacobian);
fprintf('\n\n使用函數lsqnonlin()估計得到的參數值為:\n')
fprintf('\tk1 = %.4f ± %.4f\n',k(1),ci(1,2)-k(1))
fprintf('\tk2 = %.4f ± %.4f\n',k(2),ci(2,2)-k(2))
fprintf('\tk3 = %.4f ± %.4f\n',k(3),ci(3,2)-k(3))
fprintf('\tk4 = %.4f ± %.4f\n',k(4),ci(4,2)-k(4))
fprintf('\tk5 = %.4f ± %.4f\n',k(5),ci(5,2)-k(5))
fprintf('  The sum of the squares is: %.1e\n\n',resnorm)

% ------------------------------------------------------------------
function f = ObjFunc(k,tspan,x0,yexp)           % 目標函數
[t Xsim] = ode45(@KineticsEqs,tspan,x0,[],k);   
ysim(:,1) = Xsim(2:end,1);
ysim(:,2) = Xsim(2:end,2);
ysim(:,3) = Xsim(2:end,3);
f = [ysim(:,1)-yexp(:,1); ysim(:,2)-yexp(:,2);  ysim(:,3)-yexp(:,3)];

% ------------------------------------------------------------------
function dCdt = KineticsEqs(t,C,k)              % ODE模型方程
denom = k(4)*C(1)+C(2)+k(5)*C(3);               % k(4) = KA, k(5) = KC
theA = k(4)*C(1) / denom;
theB = C(2) / denom;
r1 = k(1)*theA;
r2 = k(2)*theB;
r3 = k(3)*theB;

dCAdt = - r1 + r2;
dCBdt = r1 - r2 - r3;
dCCdt = r3;

dCdt = [dCAdt; dCBdt; dCCdt];