24小时热门版块排行榜

返回列表

tailei

银虫 (初入文坛)

应助: 0 (幼儿园)
金币: 561
帖子: 9
在线: 11小时
虫号: 1673274
注册: 2012-03-07
专业: 催化化学

[求助] 用matlab求动力学参数有错误，求指导

程序是：

tspan = [78.666667 59.000000 47.200000 39.333333  33.714286 29.500000 ];
x0=[0; 0; 0 ];
k0 = [0 0 0 0 0 0 0 0 0];%速率常数，平衡常数初值
lb = [0 0 0 0 0 0 0 0 0];
ub = [+inf +inf +inf +inf +inf +inf +inf +inf +inf];
data=...
[
78.666667 0.957100 0.797360 0.162324 70 2.5
59.000000 0.859400 0.792711 0.083877 70 2.5
47.200000 0.833500 0.764903 0.068597 70 2.5
39.333333 0.757000 0.709688 0.047313 50 1.5
33.714286 0.739600 0.699588 0.040012 50 1.5
29.500000 0.721000 0.686608 0.034392 50 1.5
];
yexp=data(:,2:4);%取矩阵2-6列数据
% 使用函数lsqnonlin()进行参数估计
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
lsqnonlin(@ObjFunc,k0,lb,ub,[],tspan,x0,yexp);
ci = nlparci(k,residual,jacobian);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.4f ± %.4f\n',k(1),ci(1,2)-k(1))
fprintf('\tk2 = %.4f ± %.4f\n',k(2),ci(2,2)-k(2))
fprintf('\tk3 = %.4f ± %.4f\n',k(3),ci(3,2)-k(3))
fprintf('\tk4 = %.4f ± %.4f\n',k(4),ci(4,2)-k(4))
fprintf('\tk5 = %.4f ± %.4f\n',k(5),ci(5,2)-k(5))
fprintf('\tk6 = %.4f ± %.4f\n',k(6),ci(6,2)-k(6))
fprintf('\tk7 = %.4f ± %.4f\n',k(7),ci(7,2)-k(7))
fprintf('\tk8 = %.4f ± %.4f\n',k(8),ci(8,2)-k(8))
fprintf('\tk9 = %.4f ± %.4f\n',k(9),ci(9,2)-k(9))
fprintf('  The sum of the squares is: %.1e\n\n',resnorm)

function f = ObjFunc(k,tspan,x0,yexp)          % 目标函数
[t Xsim] = ode45(@KineticsEqs,tspan,x0,[],k);
ysim(:,1) = Xsim(2:end,1);
ysim(:,2) = Xsim(2:end,2);
ysim(:,3) = Xsim(2:end,3);
f =(sum((ysim(:,1)-yexp(:,1))^2)+sum((ysim(:,2)-yexp(:,2))^2)
+sum((ysim(:,3)-yexp(:,3))^2))/72;

function dCdt = KineticsEqs(t,C,k)             % ODE模型方程,t为空时,C转化率、收率
P1=((1/(C(5)+1))*(1-C(2))*C(6))/(1-((1/(C(5)+1))*C(2))-((1/(C(5)+1))*C(4)));
P2=((C(5)/(C(5)+1))-2*(1/(C(5)+1)*C(3)-4*(1/(C(5)+1)*C(4)/(1-((1/(C(5)+1))*C(2))-((1/(C(5)+1))*C(4)));
P3=((1/(C(5)+1)*C(3))/(1-((1/(C(5)+1))*C(2))-((1/(C(5)+1))*C(4)));
P4=((1/(C(5)+1)*C(4))/(1-((1/(C(5)+1))*C(2))-((1/(C(5)+1))*C(4)));
P5=(((1/(C(5)+1)*C(3))+2*((1/(C(5)+1)*C(4)))/(1-((1/(C(5)+1))*C(2))-((1/(C(5)+1))*C(4)));
% k(1) = Kdmo,k(1)=Kh;k(3)=Kmg;k(4)=Keg,k(5)=Kme; k(6)=k1，k(7)=k2, k(8)=Kp1，k(9)=Kp2,
%P(1)=Pdmo,P(2)=Ph,P3=Pmg,P4=Peg,P5=Pme;
theB = 1+k(4)*P4+k5*P5+((k(1)*P3*P5*P5)/(k(8)*(P2^3));
theC =k(7)*(P3-(P3*P5)/(k(9)*(P2^2));
r1 = theA/theB;
r2 = theC/theB;

dCAdt = - r1;
dCBdt = r1 - r2 ;

dCdt = [dCAdt; dCBdt];

??? Error using ==> feval
Undefined function or method 'ObjFunc' for input arguments of type 'double'.

Error in ==> lsqnonlin at 203
         initVals.F = feval(funfcn{3},xCurrent,varargin{:});

Caused by:
Failure in initial user-supplied objective function evaluation. LSQNONLIN cannot continue.

求大神修正代码！！不胜感激！

回复此楼