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

2楼: Originally posted by chendequan at 2017-02-19 10:02:37
>> KineticsEst5

Local minimum possible. Constraints satisfied.

fmincon stopped because the size of the current step is less than
the default value of the step size tolerance and constrai ...

2楼: Originally posted by chendequan at 2017-02-19 10:02:37
>> KineticsEst5

Local minimum possible. Constraints satisfied.

fmincon stopped because the size of the current step is less than
the default value of the step size tolerance and constrai ...