马上注册,结交更多好友,享用更多功能,让你轻松玩转社区。
您需要 登录 才可以下载或查看,没有账号?我要加入
x
最近本人应用Matlab计算反应动力学参数,遇到了点问题。特此向各位高人请教。
本人的程序如下:
clear all
clc
k0 = [0.5 0.5 0.5 0.5 0.5 0.5];
lb = [0 0 0 0 0 0];
ub = [+inf +inf +inf +inf +inf +inf];
x0 = [0.7584 0.0000 0.0000 0.0000 6.8256 0.0000];
KineticsData65;
yexp = ExpData(:,2:7);
[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('\tk6 = %.4f\n',k(6))
fprintf(' The sum of the squares is: %.1e\n\n',fval)
k_fmincon = k;
[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
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 : 0.01 : 0.20];
[t x] = ode45(@KineticEqs,tspan,x0,[],k);
y(:,1) = x(:,1);
y(:,2:5) = x(:,5:8);
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)+ sum((y(:,5)-yexp(:,5)).^2);
% ------------------------------------------------------------------
function f = ObjFunc4LNL(k,x0,yexp)
tspan = [0.00 : 0.01 : 0.20];
[t x] = ode45(@KineticEqs,tspan,x0,[],k);
y(:,1) = x(:,1);
y(:,2:5) = x(:,5:8);
f1 = y(:,1) - yexp(:,1);
f2 = y(:,2) - yexp(:,2);
f3 = y(:,3) - yexp(:,3);
f4 = y(:,4) - yexp(:,4);
f5 = y(:,5) - yexp(:,5);
f = [f1; f2; f3; f4; f5];
% ------------------------------------------------------------------
function dxdt = KineticEqs(t,x,k)
dxdt = ...
[ ( -k(1)*x(1)*x(5)+k(2)*x(2)*x(4) )
( k(1)*x(1)*x(5)-k(2)*x(2)*x(4)-K(3)*x(2)*x(5)+K(4)*x(3)*x(4) )
( K(3)*x(2)*x(5)-K(4)*x(3)*x(4)-k(5)*x(3)*x(5)+k(6)*x(4)*x(6) )
( K(1)*x(1)*x(5)-K(2)*x(2)*x(4)+K(3)*x(2)*x(5)-k(4)*x(3)*x(4)+k(5)*x(3)*x(5)-k(6)*x(4)*x(6) )
( -K(1)*x(1)*x(5)+K(2)*x(2)*x(4)-K(3)*x(2)*x(5)+k(4)*x(3)*x(4)-k(5)*x(3)*x(5)+k(6)*x(4)*x(6) )
( k(5)*x(3)*x(5)-k(6)*x(4)*x(6) )
];
出现的错误提示为:
Warning: Large-scale (trust region) method does not currently solve this type of problem,
switching to medium-scale (line search).
> In fmincon at 260
In KineticsEst65 at 25
??? Error using ==> fmincon
FMINCON cannot continue because user supplied objective function failed with the following error:
Undefined command/function 'K'.
Error in ==> KineticsEst65 at 25
[k,fval,flag] = fmincon(@ObjFunc4Fmincon,k0,[],[],[],[],lb,ub,[],[],x0,yexp);
请哪位高人帮我修改一下程序,非常感谢! |