function BioKineticsEst
% 生化反應動力學參數估計
%
%   Author: HUANG Huajiang
%   Copyright 2003 UNILAB Research Center, 
%   East China University of Science and Technology, Shanghai, PRC
%   $Revision: 1.0 $  $Date: 2003/07/10 $
%
%   [Ref]: J. R. Backhurst and J. H. Harker, Chemical Engineering Vol. 5, 
%   2nd ed., Butterworth-heinemann, 1997, 世界圖書出版公司1999年重印

clear all
clc
global Y kd mum Ks D X0
% Experimental data：
%           F (l/h)  S (mg/l) X (mg/l)
ExpData = [ 0.7900   36.9000  487.0000
            1.0300   49.1000  490.0000
            1.3100   64.4000  489.0000
            1.7800   93.4000  482.0000
            2.3900  138.8000  466.0000
            2.6800  164.2000  465.0000  ];
V = 9.8;    % liter
S0 = 1200;  % limiting substrate concentration
X0 = 0;
F = ExpData(:,1);
S = ExpData(:,2);
X = ExpData(:,3);
D = F/V

% 用多變量線性回歸方法估計動力學參數
% 對方程（5）進行線性回歸
y = X./D./(S0 - S);
XX = [ones(size(y)) 1./S];
[a,aint] = regress(y,XX,0.05); 
Ks = a(2)/a(1)
Y2mum = a(1);
% 對方程（6）進行線性回歸
y = (S0 - S)./X;
XX = [ones(size(y)) 1./D];
[b,bint] = regress(y,XX,0.05); 
Y = 1/b(1)
kd = Y*b(2)
mum = Y/a(1)
fprintf('\nEstimated Parameters by Regress():\n')
fprintf('\tY = %.3f\n',Y)
fprintf('\tkd = %.3f\n',kd)
fprintf('\tmum = %.2f\n',mum)
fprintf('\tKs = %.0f\n',Ks)

% 用lsqnonlin()--求解非線性最小二乘法(非線性數據擬合)問題
lb = [0  0  0  0];
beta0 = [Y, kd, mum, Ks];
[beta,resnorm,resid,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFun,beta0,lb,[],[],S,X,D,S0,X0);
ci = nlparci(beta,resid,jacobian);

% 參數辨識結果
fprintf('\nEstimated Parameters by Lsqnonlin():\n')
fprintf('\tY = %.3f ± %.3f\n',beta(1),ci(1,2)-beta(1))
fprintf('\tkd = %.3f ± %.3f\n',beta(2),ci(2,2)-beta(2))
fprintf('\tmum = %.2f ± %.2f\n',beta(3),ci(3,2)-beta(3))
fprintf('\tKs = %.0f ± %.0f\n',beta(4),ci(4,2)-beta(4))
fprintf('  The sum of the squares is: %.1e\n\n',sum(resid.^2))

% 繪圖比較實驗值與模擬值
for i=1:length(D)
    x = fsolve(@Eqs,[S(1),X(1)],[],beta(1),beta(2),beta(3),beta(4),D(i),S0,X0);
    Scal(i) = x(1);
    Xcal(i) = x(2);
end
plot(F,S,'o',F,Scal,'*');
xlabel('Feed Flowrate (l/h)')
ylabel('Substrate Conc.(mg/l)')
figure
plot(F,X,'o',F,Xcal,'*');
xlabel('Feed Flowrate (l/h)')
ylabel('Cell Density (mg/l)')

% 殘差關于擬合值的殘差圖
plot(Scal,S'-Scal,'*')
xlabel('底物濃度擬合值, mg/l')
ylabel('底物濃度殘差R, mg/l')
refline(0,0)
figure
plot(Xcal,X'-Xcal,'*')
xlabel('干細胞密度擬合值, mg/l')
ylabel('干細胞密度殘差R, mg/l')
refline(0,0)

% ------------------------------------------------------------------
function f = ObjFun(beta,S,X,D,S0,X0)               % 目標函數
Y = beta(1);
kd = beta(2);
mum = beta(3);
Ks = beta(4);
for i=1:length(D)
    x = fsolve(@Eqs,[S(1),X(1)],[],Y,kd,mum,Ks,D(i),S0,X0);
    Scal(i) = x(1);
    Xcal(i) = x(2);
end
f = [S - Scal'; X - Xcal'];

% ------------------------------------------------------------------
function f = Eqs(x,Y,kd,mum,Ks,D,S0,X0)             % 模型方程
S = x(1);
X = x(2);
tmp = mum*S*X/(Ks + S);
f1 = D*(X - X0) - tmp + kd*X;
f2 = D*Y*(S0 - S) - tmp;
f = [f1; f2];