版块导航: 正在加载中...

登录注册

应《网络安全法》要求，自2017年10月1日起，未进行实名认证将不得使用互联网跟帖服务。为保障您的帐号能够正常使用，请尽快对帐号进行手机号验证，感谢您的理解与支持！

24小时热门版块排行榜

返回列表

当前主题已经存档。

liuyong2007

金虫 (小有名气)

应助: 0 (幼儿园)
金币: 840
红花: 1
帖子: 169
在线: 82.4小时
虫号: 400205
注册: 2007-06-12
性别: GG
专业: 化学反应工程

[交流] 【求助】用Matlab求解动力学方程参数的估计问题,程序运行有错误【已完成】

哪位帮帮我,谢谢了.我用龙阁库塔积分和非线性最小二乘估计动力学参数,但是有错误.请高手指点一下,谢谢.程序如下:

function KineticsEst1_int
% 动力学方程为rA=dCA/dt=k1*CA*CB-k2*CC
clear all
clc
global CAm CBm CCm
t = [0 20  40  60  120  180  300];
CAm = [10  8  6  5  3  2  1]';
CBm = [8  7  6  5  3  2  1]';
CCm = [0 5  7  8 10  12  14]';
% 非线性拟合
beta0 = [0.0053 1.39];
tspan = [0 20  40  60  120  180  300];
CA0 = 10;
[beta,resnorm,resid,exitflag,output,lambda,jacobian] = ...
lsqnonlin(@OptObjFunc,beta0,[],[],[],tspan,CA0,CAm,CBm,CCm)
ci = nlparci(beta,resid,jacobian)
% 拟合效果图(实验与拟合的比较)
[t4plot CA4plot] = ode45(@KineticsEqs,[tspan(1)  tspan(end)],CA0,[],beta);
plot(t,CAm,'bo',t4plot,CA4plot,'k-')
legend('Exp','Model')
xlabel('时间t, s')
ylabel('浓度C_A, mol/L')
% 参数辨识结果
fprintf('Estimated Parameters:\n')
fprintf('\tk1 = %.4f ± %.4f\n',beta(1),ci(1,2)-beta(1))
fprintf('\tk2 = %.4f ± %.4f\n',beta(2),ci(2,2)-beta(2))
% ------------------------------------------------------------------
function f = OptObjFunc(beta,tspan,CA0,CAm)
[t CAc] = ode45(@KineticsEqs,tspan,CA0,[],beta);
f = CAc - CAm;
% ------------------------------------------------------------------
function dCAdt = KineticsEqs(t,CA,CB,CC,beta)
dCAdt = beta(1)*CA*CB-beta(2)*CC;          % k1= beta(1), k2= beta(2)

[ Last edited by nono2009 on 2009-9-23 at 21:15 ]

回复此楼

» 收录本帖的淘帖专辑推荐

matlab

» 猜你喜欢

1楼 2008-12-11 09:26:44

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

allenhero1228

金虫 (小有名气)

应助: 0 (幼儿园)
金币: 1617.5
帖子: 201
在线: 29.7小时
虫号: 420486
注册: 2007-07-10
性别: GG
专业: 化工/二甲醚催化剂研究

★ ★
woshilsh(金币+2,VIP+0):thanks for your help!

你的目标函数是KineticsEqs(t,CA,CB,CC,beta)，beta是你要拟合的参数，剩下的是形参，需要你输入。以你现在的顺序在计算过程中可能是这样，tspan赋给t，CA0赋给CA，CAM赋给CB，CBm赋给CC，CCm赋给beta，这样程序里面就没有你要拟合的参数，程序就没法进行计算。
好像是应该把beta放在最前面
以上是我的想法，不知道对不对，希望后面有牛人给你正确详细的回答
另外建议你把matlab中给你的错误提示也发上来

赞一下(10人)

回复此楼

2楼2008-12-11 12:33:19

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

liuyong2007

金虫 (小有名气)

应助: 0 (幼儿园)
金币: 840
红花: 1
帖子: 169
在线: 82.4小时
虫号: 400205
注册: 2007-06-12
性别: GG
专业: 化学反应工程

mat;ab提示错误信息,希望高手给我改改

??? Error using ==> optim\private\lsqncommon
User supplied function ==> OptObjFunc
failed with the following error:

Error using ==> KineticsEst1_int22>OptObjFunc
Too many input arguments.

Error in ==> lsqnonlin at 163
[x,Resnorm,FVAL,EXITFLAG,OUTPUT,LAMBDA,JACOB] = ...

Error in ==> KineticsEst1_int22 at 17
[beta,resnorm,resid,exitflag,output,lambda,jacobian] = ...

赞一下

回复此楼

3楼2008-12-11 12:41:34

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

allenhero1228

金虫 (小有名气)

应助: 0 (幼儿园)
金币: 1617.5
帖子: 201
在线: 29.7小时
虫号: 420486
注册: 2007-07-10
性别: GG
专业: 化工/二甲醚催化剂研究

★
coldwind042(金币+1,VIP+0):谢谢帮助！

你试试把beta放在最前面，然后在global中把CA CB CC加上，而且我看你想把CA0,CAm,CBm,CCm都赋给函数，但函数里面只有三个形参，是不是应该把CA0,CAm,CBm,CCm中的某个变量给去了

赞一下(10人)

回复此楼

4楼2008-12-11 13:06:55

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

snipher950

木虫 (正式写手)

应助: 2 (幼儿园)
金币: 3024.6
红花: 3
帖子: 797
在线: 22.5小时
虫号: 522854
注册: 2008-03-12
性别: MM
专业: 化工系统工程

你自定义函数之间的参数传递有问题。给你改改……

赞一下

回复此楼

5楼2008-12-11 14:20:27

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

snipher950

木虫 (正式写手)

应助: 2 (幼儿园)
金币: 3024.6
红花: 3
帖子: 797
在线: 22.5小时
虫号: 522854
注册: 2008-03-12
性别: MM
专业: 化工系统工程

★ ★ ★ ★ ★
woshilsh(金币+5,VIP+0):感谢回答！支持热心虫友，常来！呵呵！

可以运行了，改了很多地方，你自己琢磨一下。只不过结果不对，那是因为你的实验数据是虚假的，Matlab拟和不出来好的结果。

function KineticsEst1_int
% 动力学方程为rA=dCA/dt=k1*CA*CB-k2*CC
clear all
clc
global CAm CBm CCm
t = [0 20  40  60  120  180  300];
CAm = [10  8  6  5  3  2  1]';
CBm = [8  7  6  5  3  2  1]';
CCm = [0 5  7  8 10  12  14]';
% 非线性拟合
beta0 = [0.0000053 10];
tspan = [0 20  40  60  120  180  300];
CA0 = 10;
CB0 = 8;
CC0 = 0;
[beta,resnorm,resid,exitflag,output,lambda,jacobian] = ...
lsqnonlin(@OptObjFunc,beta0,[0 0],[],[],tspan,[CA0 CB0 CC0])
ci = nlparci(beta,resid,jacobian)
% 拟合效果图(实验与拟合的比较)
[t4plot CA4plot] = ode45(@KineticsEqs,[tspan(1)  tspan(end)],[CA0 CB0 CC0],[],beta);
plot(tspan,CAm,'bo',t4plot,CA4plot(:,1),'k-')
legend('Exp','Model')
xlabel('时间t, s')
ylabel('浓度C_A, mol/L')
% 参数辨识结果
fprintf('Estimated Parameters:\n')
fprintf('\tk1 = %.4f ± %.4f\n',beta(1),ci(1,2)-beta(1))
fprintf('\tk2 = %.4f ± %.4f\n',beta(2),ci(2,2)-beta(2))
% ------------------------------------------------------------------
function f = OptObjFunc(beta,tspan,CA0)
global CAm CBm CCm
[t Y] = ode45(@KineticsEqs,tspan,CA0,[],beta);
if length(t)==length(tspan)
f = Y - [CAm CBm CCm];
else
f=100000;
end
% ------------------------------------------------------------------
function dy = KineticsEqs(t,Y,beta)
CA=Y(1);
CB=Y(2);
CC=Y(3);
dCAdt = beta(1)*CA*CB-beta(2)*CC;          % k1= beta(1), k2= beta(2)
dCBdt = beta(1)*CA*CB-beta(2)*CC;
dCCdt = beta(2)*CC;
dy=[dCAdt;dCBdt;dCCdt];

赞一下(10人)

回复此楼

6楼2008-12-11 14:35:07

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