24小时热门版块排行榜    

查看: 705  |  回复: 2

llgcty

金虫 (正式写手)

[求助] 在线等,请版主看看如何运算不了。总是报错

小弟实在是一窍不通对matlab,不要怪小弟提的问题比较傻,一个是报错中的问题,这个[k,fval,flag] = fmincon(@ObjFunc7Fmincon一句的@ObjFunc7Fmincon是怎么设定的?我看有的人写@ObjFunc4Fmincon.
然后麻烦大家看看这个我自己模拟的代码,是按照版主 @dbb627的代码写的,到底是哪里不对麻烦改下谢谢了。
function parafit
%  
% r1 = k(1)*C(1);
% r2 = k(2)*C(1);
% r3 = k(3)*C(2);
% r4 = k(4)*C(3);
% r5 = k(5)*C(2);
% r6 = k(6)*C(3);
%
% dCAdt = - r1 - r2;
% dCBdt = r1 + r4 - r3 - r5;
% dCCdt = r2 + r3 - r4 – r6 ;
% dCDdt = r5 + r6;
clear all
clc
%        t/min   CA-1     CB-2        CC-3   CD-4   / mol/L
  Kinetics=[0       100    0         0          0      
          10        0.2441  0.1390   0.28.31   0.0171   
          20        0.432  0.1629   0.2370   0.0192   
          30        0.756   0.1992  0.2509    0.0349   
          40        0.522   0.2091   0.2737    0.0531   
          50        0.426    0.1967   0.3231   0.0634   
          60        0.397    0.1803   0.3536   0.0822   
          70        0.431    0.1732   0.3330   0.0945  
          80        0.531   0.1645    0.3310  0.0987];
k0 = [0.0000000005  0.0000000005  0.0000000005  0.00000000005  0.00000005  0.00000005];         % 参数初值
lb = [0  0  0  0  0  0];                   % 参数下限
ub = [1 1  1  1  1  1];    % 参数上限
x0 = [100  0  0  0];
yexp = Kinetics;                  % yexp: 实验数据[x1        x4        x5        x6]
warning off
% 使用函数fmincon()进行参数估计
[k,fval,flag] = fmincon(@ObjFunc7Fmincon,k0,[],[],[],[],lb,ub,[],[],x0,yexp);
fprintf('\n使用函数fmincon()估计得到的参数值为:\n')
fprintf('\tk1 = %.11f\n',k(1))
fprintf('\tk2 = %.11f\n',k(2))
fprintf('\tk3 = %.11f\n',k(3))
fprintf('\tk4 = %.11f\n',k(4))
fprintf('\tk5 = %.11f\n',k(5))
fprintf('\tk6 = %.11f\n',k(6))
fprintf('  The sum of the squares is: %.1e\n\n',fval)
k_fm= k;
warning off
% 使用函数lsqnonlin()进行参数估计
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc7LNL,k0,lb,ub,[],x0,yexp);      
ci = nlparci(k,residual,jacobian);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.11f\n',k(1))
fprintf('\tk2 = %.11f\n',k(2))
fprintf('\tk3 = %.11f\n',k(3))
fprintf('\tk4 = %.11f\n',k(4))
fprintf('\tk5 = %.11f\n',k(5))
fprintf('\tk6 = %.11f\n',k(6))
fprintf('  The sum of the squares is: %.1e\n\n',resnorm)
k_ls = k;
output
warning off
% 以函数fmincon()估计得到的结果为初值,使用函数lsqnonlin()进行参数估计
k0 = k_fm;
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc7LNL,k0,lb,ub,[],x0,yexp);      
ci = nlparci(k,residual,jacobian);
fprintf('\n\n以fmincon()的结果为初值,使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.11f\n',k(1))
fprintf('\tk2 = %.11f\n',k(2))
fprintf('\tk3 = %.11f\n',k(3))
fprintf('\tk4 = %.11f\n',k(4))
fprintf('\tk5 = %.11f\n',k(5))
fprintf('\tk6 = %.11f\n',k(6))
fprintf('  The sum of the squares is: %.1e\n\n',resnorm)
k_fmls = k;
output
tspan = [0 : 10 : 80];
[t x] = ode45(@KineticEqs,tspan,x0,[],k_fmls);
figure;
plot(t,x(:,1),t,yexp(:,2),'*');legend('ca-pr','ca-real')
figure;plot(t,x(:,2:5));
hold on
plot(t,yexp(:,3:6),'o');legend('cb-pr','cc-pr','cd-pr','ce-pr','cb-real','cc-real','cd-real','ce-real')
% ------------------------------------------------------------------
function f = ObjFunc7Fmincon(k,x0,yexp)
tspan = [0 : 10 : 80];
[t x] = ode45(@KineticEqs,tspan,x0,[],k);   
y(:,2) = x(:,1);
y(:,3:6) = x(:,2:5);
f =  sum((y(:,2)-yexp(:,2)).^2) + sum((y(:,3)-yexp(:,3)).^2)   ...
    + sum((y(:,4)-yexp(:,4)).^2) + sum((y(:,5)-yexp(:,5)).^2)   ...
    + sum((y(:,6)-yexp(:,6)).^2) ;

% ------------------------------------------------------------------
function f = ObjFunc7LNL(k,x0,yexp)
tspan = [0.0 : 10 : 80];
[t x] = ode45(@KineticEqs,tspan,x0,[],k);   
y(:,2) = x(:,1);
y(:,3:6) = x(:,2:5);
f1 = y(:,2) - yexp(:,2);
f2 = y(:,3) - yexp(:,3);
f3 = y(:,4) - yexp(:,4);
f4 = y(:,5) - yexp(:,5);
f5 = y(:,6) - yexp(:,6);
f = [f1; f2; f3; f4; f5];

% ------------------------------------------------------------------
function dxdt = KineticEqs(t,x,k)
r1 = k(1)*x(1);
r2 = k(2)*x(1);
r3 = k(3)*x(2);
r4 = k(4)*x(3);
r5 = k(5)*x(2);
r6 = k(6)*x(3);

dCAdt = - r1 - r2;
dCBdt = r2 + r4 - r3 - r5;
dCCdt = r2 + r3 - r4 - r6;
dCDdt = r5 + r6;

dxdt = [dCAdt; dCBdt; dCCdt; dCDdt];
回复此楼

» 本帖@通知

» 猜你喜欢

已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

llgcty

金虫 (正式写手)

快点来吧,帮忙一下!!!
2楼2015-11-29 21:52:01
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖

llgcty

金虫 (正式写手)

各位大神帮帮忙
3楼2015-11-30 11:21:13
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 llgcty 的主题更新
最具人气热帖推荐 [查看全部] 作者 回/看 最后发表
[考博] 售SCI一区文章,我:8.O.55.1.O.54,科目齐全,可伽急 +3 jCd0dEvKHShX 2026-08-29 3/150 2026-08-29 23:14 by ZPa0EcMwuECS
[基金申请] 面上意见出来了 +8 黄鸟于飞Chao 2026-08-29 16/800 2026-08-29 20:45 by 黄鸟于飞Chao
[论文投稿] 售SCI一区文章,我:8.O.551.O.5.4,科目全,可伽急 +3 gy1nBQXYQJqL 2026-08-29 3/150 2026-08-29 18:26 by 4FFAWE8HcgUD
[基金申请] 有没有仍没收到信息的 +5 德尚中行 2026-08-27 6/300 2026-08-29 18:10 by manplx
[考博] 售SCI一区T0P文章,我:8.O.55.1.O54,科目全,可伽急 +5 ASdOkHsho7FD 2026-08-28 8/400 2026-08-29 17:22 by 4FFAWE8HcgUD
[硕博家园] 售SCI一区文章,我:8.O.55.1.O.54,科目齐全,可伽急 +5 ASdOkHsho7FD 2026-08-28 9/450 2026-08-29 14:03 by jCd0dEvKHShX
[教师之家] 售SCI一区T0P文章,我:8.O.55.1.O54,科目全,可伽急 +4 ASdOkHsho7FD 2026-08-28 5/250 2026-08-29 11:49 by jCd0dEvKHShX
[考研] 售SCI一区文章,我:8.O.55.1.O.54,科目齐全,可伽急 +6 ASdOkHsho7FD 2026-08-28 9/450 2026-08-29 11:18 by jCd0dEvKHShX
[基金申请] 我就是申请一个面上项目而已,这评审意见是按照杰青的条件评的吧? +6 gouxfjh 2026-08-28 9/450 2026-08-29 09:59 by jklily
[基金申请] 基金系统什么内容也没有 30+4 winsaint 2026-08-27 9/450 2026-08-28 11:06 by maolC
[基金申请] 怎么查啊 +6 huang1991js 2026-08-26 6/300 2026-08-28 08:42 by winsaint
[文学芳草园] 梦想 +3 myrtle 2026-08-26 3/150 2026-08-27 10:01 by angelyueyi
[基金申请] 我不理解! +15 Edward_pc 2026-08-26 23/1150 2026-08-26 20:34 by zzuzxg
[基金申请] 为什么国自然不能直接公布 +4 bjdxyxy 2026-08-26 4/200 2026-08-26 13:12 by qingmu1201
[基金申请] 国际合作可查了,中了面上 (EPI+1)(金币+50) +18 Ldrop2023 2026-08-26 18/900 2026-08-26 11:15 by cmrandy
[基金申请] 系统进不去 +4 yanglien 2026-08-26 5/250 2026-08-26 11:10 by wenfengw83
[基金申请] 今天务委会开完了,明天出结果吗 +19 angus9576 2026-08-25 23/1150 2026-08-26 10:03 by zp519
[基金申请] 国合现在查不到了吗? +10 chengyan1220 2026-08-24 20/1000 2026-08-26 08:57 by peasantsprig
[基金申请] 在坚冰还盖着北海的时候,我看到了怒放的梅花。 (金币+10) +6 ziyangfang 2026-08-25 9/450 2026-08-25 20:26 by huagongfeihu
[基金申请] 没有任何消息-是不是就凉了 +9 图啦图啦 2026-08-24 10/500 2026-08-25 11:59 by 南海小哥
信息提示
请填处理意见