function KineticsEst5
% 動力學ODE方程模型的參數估計
%
%   Author: HUANG Huajiang
%   Copyright 2003 UNILAB Research Center, 
%   East China University of Science and Technology, Shanghai, PRC
%   $Revision: 1.0 $  $Date: 2003/06/06 $
%
% [Ref]：Frerich Keil, et al. ed., Scientific computing in chemical 
%   engineering II,1999 （P.351）
%
%   The variables  y here are y(1)=x1, y(2)=x4, y(3)=x5,y(4)=x6 . 

clear all
clc

k0 = [0.5  0.5  0.5  0.5  0.5];         % 參數初值
lb = [0  0  0  0  0];                   % 參數下限
ub = [+inf  +inf  +inf  +inf  +inf];    % 參數上限
x0 = [0.1883  0.2507  0.0467  0.0899  0.1804  0.1394  0.1046];
KineticsData1;
yexp = ExpData(:,2:5);                  % yexp: 實驗數據[x1	x4	x5	x6]

% 使用函數fmincon()進行參數估計
[k,fval,flag] = fmincon(@ObjFunc4Fmincon,k0,[],[],[],[],lb,ub,[],[],x0,yexp);
fprintf('\n使用函數fmincon()估計得到的參數值為:\n')
fprintf('\tk1 = %.4f\n',k(1))
fprintf('\tk2 = %.4f\n',k(2))
fprintf('\tk3 = %.4f\n',k(3))
fprintf('\tk4 = %.4f\n',k(4))
fprintf('\tk5 = %.4f\n',k(5))
fprintf('  The sum of the squares is: %.1e\n\n',fval)
k_fmincon = k;

% 使用函數lsqnonlin()進行參數估計
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[],x0,yexp);       
ci = nlparci(k,residual,jacobian);
fprintf('\n\n使用函數lsqnonlin()估計得到的參數值為:\n')
Output

% 以函數fmincon()估計得到的結果為初值，使用函數lsqnonlin()進行參數估計
k0 = k_fmincon;
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[],x0,yexp);       
ci = nlparci(k,residual,jacobian);
fprintf('\n\n以fmincon()的結果為初值，使用函數lsqnonlin()估計得到的參數值為:\n')
Output


% ------------------------------------------------------------------
function f = ObjFunc4Fmincon(k,x0,yexp)
tspan = [0.00 : 0.01 : 0.20];
[t x] = ode45(@KineticEqs,tspan,x0,[],k);   
y(:,1) = x(:,1);
y(:,2:4) = x(:,4:6);
f = sum((y(:,1)-yexp(:,1)).^2) + sum((y(:,2)-yexp(:,2)).^2)   ...
    + sum((y(:,3)-yexp(:,3)).^2) + sum((y(:,4)-yexp(:,4)).^2);

% ------------------------------------------------------------------
function f = ObjFunc4LNL(k,x0,yexp)
tspan = [0.00 : 0.01 : 0.20]; 
[t x] = ode45(@KineticEqs,tspan,x0,[],k);   
y(:,1) = x(:,1);
y(:,2:4) = x(:,4:6);
f1 = y(:,1) - yexp(:,1); 
f2 = y(:,2) - yexp(:,2);
f3 = y(:,3) - yexp(:,3);
f4 = y(:,4) - yexp(:,4);
f = [f1; f2; f3; f4];

% ------------------------------------------------------------------
function dxdt = KineticEqs(t,x,k)
q = 8.75 + k(5);
dxdt =  ...
 [ ( k(5)-q*x(1)- k(1)*x(1)*x(2)-k(4)*x(1)*x(6)*sqrt(0.9) ) 
   ( 7.0-q*x(2) - k(1)*x(1)*x(2)-2*k(2)*x(2)*x(3) )
   ( 1.75 -q*x(3) - k(2)*x(2)*x(3) )
   ( -q*x(4) + 2*k(1)*x(1)*x(2)-k(3)*x(4)*x(5) )
   ( -q*x(5) + 3*k(2)*x(2)*x(3)-k(3)*x(4)*x(5) )
   ( -q*x(6) + 2*k(3)*x(4)*x(5)-k(4)*x(1)*x(6)*sqrt(0.9) )
   ( -q*x(7) + 2*k(4)*x(1)*x(6)*sqrt(0.9) )
];