function FermKineticsEst
% 青霉素發酵過程動力學（fermentation kinetics）
%
%   Author: HUANG Huajiang
%   Copyright 2003 UNILAB Research Center, 
%   East China University of Science and Technology, Shanghai, PRC
%   $Revision: 1.0 $  $Date: 2003/04/09 $
%
%   [Ref]：Alkis Constantinides & Navid Mostoufi. Numerical Methods for 
%   Chemical Engineers with MATLAB Applications. Prentice Hall, 1999 

clear all; clc
load FermKineticsData           % Load experimental data

% Sorting data
for i = 1:length(n);
    [x(1:n(i),i),index] = sort(x(1:n(i),i));
    y(1:n(i),i) = y(index,i);
end
w = Weights(x,y,n)

% 用lsqnonlin()進行非線性最小二乘參數估計(非線性數據擬合)
beta0 = [0.1  4.0  0.02  0.02];
y0 = [0.29  0.0];
lb = [0  0.01  0  0];
ub = [inf  inf  inf  inf];
[beta,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@Func,beta0,lb,ub,[],x,y,n,y0,w);          
ci = nlparci(beta,residual,jacobian);

% 參數辨識結果
fprintf('\nEstimated Parameters by Lsqnonlin():\n')
fprintf('\tk1 = %.4f ± %.4f\n',beta(1),ci(1,2)-beta(1))
fprintf('\tk2 = %.2f ± %.2f\n',beta(2),ci(2,2)-beta(2))
fprintf('\tk3 = %.4f ± %.4f\n',beta(3),ci(3,2)-beta(3))
fprintf('\tk4 = %.2f ± %.2f\n',beta(4),ci(4,2)-beta(4))
fprintf('  The sum of the squares is: %.1e\n\n',sum(residual.^2))


% ------------------------------------------------------------------
function f = Func(beta,x,y,n,y0,w)                  % 定義目標函數
tspan = [0  max(max(x))];
[tt yy] = ode45(@ModelEqs,tspan,y0,[],beta);   
for col = 1:length(n)
    yc(1:n(col),col) = spline(tt,yy(:,col),x(1:n(col),col));
end
f1 = y(:,1)-yc(:,1);
f2 = y(:,2)-yc(:,2);
f = [sqrt(w(1))*f1; sqrt(w(2))*f2];

% ------------------------------------------------------------------
function dydt = ModelEqs(t,y,beta)                 % Model equations
dydt = [beta(1)*y(1)*(1-y(1)/beta(2));
   beta(3)*y(1)-beta(4)*y(2)];