CODE:

function KineticsEst5 % 動力學ODE方程模型的參數估計

%

%

%

% The variables y here are y(1)=xB, y(2)=xoNB, y(3)=xmNB,y(4)=xpNB,y(5)=xDNB .

clear all

clc

k0 = [5 5 5 5 5]; % 參數初值

lb = [0 0 0 0 0]; % 參數下限

ub = [+inf +inf +inf +inf +inf]; % 參數上限

x0 = [0 0 0 0 0 0];

KineticsData;

yexp = ExpData(:,2:6); % 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 : 5.03 : 15.09];

[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 : 5.03 : 15.09];

[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)

M=1.003;

dxdt = ...

[ ( (k(1)+k(2)+k(3))*(1-x(1))*M*(1-x(2)-x(3)-x(4)-2*x(5)) )

( M*(1-x(2)-x(3)-x(4)-2*x(5))*(k(1)*(1-x(1))-k(4)*x(2)) )

( k(2)*(1-x(1))*M*(1-x(2)-x(3)-x(4)-2*x(5)) )

( M*(1-x(2)-x(3)-x(4)-2*x(5))*(k(3)*(1-x(1))-k(5)*x(4)) )

( M*(1-x(2)-x(3)-x(4)-2*x(5))*(k(4)*x(2)+k(5)*x(4)) )

];