24小时热门版块排行榜

>论坛更新日志 (3900)
>文献求助 (405)
>虫友互识 (340)
>导师招生 (316)
>博后之家 (160)
>硕博家园 (154)
>考博 (142)
>休闲灌水 (141)
>论文投稿 (114)
>基金申请 (90)
>招聘信息布告栏 (54)
>考研 (54)
>教师之家 (51)
>公派出国 (43)
>绿色求助(高悬赏) (37)
>找工作 (37)

返回列表

当前只显示满足指定条件的回帖，点击这里查看本话题的所有回帖

楚天湘水

金虫 (小有名气)

应助: 1 (幼儿园)
金币: 191.5
散金: 676
红花: 2
帖子: 123
在线: 131.2小时
虫号: 1167266
注册: 2010-12-10
专业: 化工系统工程

[求助] matlab 拟合反应动力学参数结果很差。大家帮忙看一下

已知动力学模型为：
dA1dt = r1 + r2;
dA2dt =  - r1-r2;
dA3dt = -r1；
dA4dt = r1 - r2;
dA5dt = r2;

r1 = k(1)*((a(3)*a(2) - (1/2.20)*a(4)*a(1)));
r2 = k(2)*(a(4)*a(2) -  (1/0.40)*a(5)*a(1)));

已知的参数：
时间：tspan = [0 60 120 180 240 300 360 420 600 1200 6900 7500],
活度系数：KineticsData1
% 动力学数据
%  a1 a2 a3 a4 a5
ExpData=...
[0.0000 0.8122 0.3894 0.0000 0.0000
   0.1827 0.6353 0.1599 0.2163 0.0444
   0.2547 0.5668 0.0979 0.2485 0.0805
   0.2854 0.5390 0.0788 0.2415 0.1061
   0.2978 0.5258 0.0738 0.2359 0.1192
   0.3019 0.5210 0.0710 0.2309 0.1288
   0.3092 0.5143 0.0691 0.2264 0.1349
   0.3178 0.5062 0.0680 0.2223 0.1349
   0.3178 0.5062 0.0680 0.2186 0.1433
   0.3211 0.5030 0.0674 0.2154 0.1478
   0.3209 0.5036 0.0673 0.2140 0.1482
   0.3205 0.5045 0.0674 0.2136 0.1476 ]

下面是写的matlab程序：
function Kinetic2
% 动力学ODE方程模型的参数估计
clear all
clc

k0 = [900000 50000];       % 参数初值
lb = [0  0];                % 参数下限
ub = [+inf  +inf ]; % 参数上限
a0 = [0.0000 0.8122 0.3894 0.0000 0.0000];
KineticsData1;
yexp = ExpData(:,1:5); % yexp: 实验数据[a1 a2 a3 a4 a5]

% 使用函数fmincon()进行参数估计
[k,fval,flag] = fmincon(@ObjFunc4Fmincon,k0,[],[],[],[],lb,ub,[],[],a0,yexp);
fprintf('\n使用函数fmincon()估计得到的参数值为:\n')
fprintf('\tk1 = %.4f\n',k(1))
fprintf('\tk2 = %.4f\n',k(2))
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,[],a0,yexp);
ci = nlparci(k,residual,jacobian);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.4f\n',k(1))
fprintf('\tk2 = %.4f\n',k(2))
fprintf('  The sum of the squares is: %.1e\n\n',resnorm)
exitflag

% 以函数fmincon()估计得到的结果为初值，使用函数lsqnonlin()进行参数估计
k0 = k_fmincon;
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[1e-15],a0,yexp);
ci = nlparci(k,residual,jacobian);
fprintf('\n\n以fmincon()的结果为初值，使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk1 = %.4f\n',k(1))
fprintf('\tk2 = %.4f\n',k(2))
fprintf('  The sum of the squares is: %.1e\n\n',resnorm)
exitflag

function f = ObjFunc4Fmincon(k,a0,yexp)
tspan = [0 60 120 180 240 300 360 420 600 1200 6900 7500];
[t a] = ode23s(@KineticEqs,tspan,a0,[],k);
y(:,1:5) = a(:,1:5);
f = sum(((y(:,1)-yexp(:,1))/0.01).^2) + sum(((y(:,2)-yexp(:,2))/0.01).^2) ...
+ sum(((y(:,3)-yexp(:,3))/0.005).^2) + sum(((y(:,4)-yexp(:,4))/0.01).^2) ...
+ sum(((y(:,5)-yexp(:,5))/0.01).^2);

function f = ObjFunc4LNL(k,a0,yexp)
tspan = [0 60 120 180 240 300 360 420 600 1200 6900 7500];
[t a] = ode23s(@KineticEqs,tspan,a0,[],k);
y(:,1:5) = a(:,1:5);
f1 = (y(:,1) - yexp(:,1))/0.01;
f2 = (y(:,2) - yexp(:,2))/0.01;
f3 = (y(:,3) - yexp(:,3))/0.005;
f4 = (y(:,4) - yexp(:,4))/0.01;
f5 = (y(:,5) - yexp(:,5))/0.01;
f = [f1; f2; f3; f4; f5];
plot(tspan,yexp(:,1), 'b^',tspan,yexp(:,2), 'ro',tspan,yexp(:,3),'y*',...
tspan,yexp(:,4), 'go',tspan,yexp(:,5), 'cs')
hold on
plot(tspan,y(:,1), 'b-',tspan,yexp(:,2), 'r-',tspan,yexp(:,3),'y-', ...
tspan,yexp(:,4), 'g-',tspan,yexp(:,5), 'c-')
axis([0 800 0 1])

%动力学模型-----------------------------------------------------------------------
function dadt = KineticEqs(t,a ,k)
Keq=[2.20 0.40];
dadt=...
[((k(1)*(a(3)*a(2)-(1/Keq(1))*a(4)*a(1))) + (k(2)*(a(4)*a(2)-(1/Keq(2))*a(5)*a(1))));
-((k(1)*(a(3)*a(2)-(1/Keq(1))*a(4)*a(1))) + (k(2)*(a(4)*a(2)-(1/Keq(2))*a(5)*a(1))));
-(k(1)*(a(3)*a(2)-(1/Keq(1))*a(4)*a(1)));
  ((k(1)*(a(3)*a(2)-(1/Keq(1))*a(4)*a(1))) - (k(2)*(a(4)*a(2)-(1/Keq(2))*a(5)*a(1))));
(k(2)*(a(4)*a(2)-(1/Keq(2))*a(5)*a(1)))];

大家帮忙看一看，哪里有问题，从做出的图可以看出，处理的结果很差啊。
无论K初值给出怎么样，都是拟合出来的a1和实验值a1相差很远，而其他值都拟合很好。

[ Last edited by 楚天湘水 on 2012-5-1 at 00:44 ]

回复此楼

1楼2012-05-01 00:38:48

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

dingd

铁杆木虫 (职业作家)

计算强帖: 4
应助: 1641 (讲师)
金币: 15037.3
散金: 101
红花: 234
帖子: 3410
在线: 1223.5小时
虫号: 291104
注册: 2006-10-28

【答案】应助回帖

★ ★
dbb627: 金币+2, 谢谢应助！ 2012-05-02 14:06:30

见下：

CODE:

ConstStr r1 = k1*(a3*a2-(1/2.20)*a4*a1),

         r2 = k2*(a4*a2-(1/0.40)*a5*a1);

Variable t,a1,a2,a3,a4,a5;

ODEFunction A1' = r1 + r2;

            A2' =  -r1-r2;

            A3' = -r1;

            A4' = r1 - r2;

            A5' = r2;

Data;

0        0.0000        0.8122        0.3894        0.0000        0.0000

60        0.1827        0.6353        0.1599        0.2163        0.0444

120        0.2547        0.5668        0.0979        0.2485        0.0805

180        0.2854        0.5390        0.0788        0.2415        0.1061

240        0.2978        0.5258        0.0738        0.2359        0.1192

300        0.3019        0.5210        0.0710        0.2309        0.1288

360        0.3092        0.5143        0.0691        0.2264        0.1349

420        0.3178        0.5062        0.0680        0.2223        0.1349

600        0.3178        0.5062        0.0680        0.2186        0.1433

1200        0.3211        0.5030        0.0674        0.2154        0.1478

6900        0.3209        0.5036        0.0673        0.2140        0.1482

7500        0.3205        0.5045        0.0674        0.2136        0.1476