function KineticsEst

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,output] = 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;

% ------------------------------------------------------------------

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

];