楚天湘水

應(yīng)助: 1 (幼兒園)
金幣: 191.5
散金: 676
紅花: 2
帖子: 123
在線: 131.2小時(shí)
蟲(chóng)號(hào): 1167266
注冊(cè): 2010-12-10
專業(yè): 化工系統(tǒng)工程

[求助] matlab 擬合反應(yīng)動(dòng)力學(xué)參數(shù)結(jié)果很差。大家?guī)兔匆幌?

已知?jiǎng)恿W(xué)模型為：
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)));

已知的參數(shù)：
時(shí)間：tspan = [0 60 120 180 240 300 360 420 600 1200 6900 7500],
活度系數(shù)：KineticsData1
% 動(dòng)力學(xué)數(shù)據(jù)
%  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 ]

下面是寫(xiě)的matlab程序：
function Kinetic2
% 動(dòng)力學(xué)ODE方程模型的參數(shù)估計(jì)
clear all
clc

k0 = [900000 50000];       % 參數(shù)初值
lb = [0  0];                % 參數(shù)下限
ub = [+inf  +inf ]; % 參數(shù)上限
a0 = [0.0000 0.8122 0.3894 0.0000 0.0000];
KineticsData1;
yexp = ExpData(:,1:5); % yexp: 實(shí)驗(yàn)數(shù)據(jù)[a1 a2 a3 a4 a5]

% 使用函數(shù)fmincon()進(jìn)行參數(shù)估計(jì)
[k,fval,flag] = fmincon(@ObjFunc4Fmincon,k0,[],[],[],[],lb,ub,[],[],a0,yexp);
fprintf('\n使用函數(shù)fmincon()估計(jì)得到的參數(shù)值為:\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;

% 使用函數(shù)lsqnonlin()進(jìn)行參數(shù)估計(jì)
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[],a0,yexp);
ci = nlparci(k,residual,jacobian);
fprintf('\n\n使用函數(shù)lsqnonlin()估計(jì)得到的參數(shù)值為:\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

% 以函數(shù)fmincon()估計(jì)得到的結(jié)果為初值，使用函數(shù)lsqnonlin()進(jìn)行參數(shù)估計(jì)
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()的結(jié)果為初值，使用函數(shù)lsqnonlin()估計(jì)得到的參數(shù)值為:\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])

%動(dòng)力學(xué)模型-----------------------------------------------------------------------
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)))];

大家?guī)兔匆豢�，哪里有�?wèn)題，從做出的圖可以看出，處理的結(jié)果很差啊。
無(wú)論K初值給出怎么樣，都是擬合出來(lái)的a1和實(shí)驗(yàn)值a1相差很遠(yuǎn)，而其他值都擬合很好。

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

回復(fù)此樓

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

已閱回復(fù)此樓關(guān)注TA 給TA發(fā)消息送TA紅花 TA的回帖

回帖支持 ( 顯示支持度最高的前 50 名 )

feng198800

金蟲(chóng) (正式寫(xiě)手)

應(yīng)助: 15 (小學(xué)生)
金幣: 860.2
散金: 234
紅花: 1
帖子: 415
在線: 242.1小時(shí)
蟲(chóng)號(hào): 700353
注冊(cè): 2009-02-12
性別: GG
專業(yè): 化學(xué)反應(yīng)工程

【答案】應(yīng)助回帖

★ ★
zhangguangping: 金幣+2, 謝謝應(yīng)助！ 2012-05-29 18:54:39

哦，看明白了！呵呵。最小二乘法可以直接擬合，你運(yùn)用非線性約束函數(shù)fmincon尋找初值也只是該區(qū)域的最小值點(diǎn)，和自己估計(jì)初值直接用最小二乘法擬合擬合的結(jié)果是一樣的（見(jiàn)樓主的擬合結(jié)果都是The sum of the squares is: 3.3e+003），有可能都是局部最優(yōu)的，可以不用fmincon，結(jié)合全局最優(yōu)化方法尋找初值，例如遺傳算法或者PSO算法，個(gè)人觀點(diǎn)，僅供參考，歡迎大家指點(diǎn)！

贊一下(1人)

回復(fù)此樓

10樓2012-05-29 11:50:35

已閱回復(fù)此樓關(guān)注TA 給TA發(fā)消息送TA紅花 TA的回帖

普通回帖

dingd

鐵桿木蟲(chóng) (職業(yè)作家)

計(jì)算強(qiáng)帖: 4
應(yīng)助: 1641 (講師)
金幣: 15037.3
散金: 101
紅花: 234
帖子: 3410
在線: 1223.7小時(shí)
蟲(chóng)號(hào): 291104
注冊(cè): 2006-10-28

【答案】應(yīng)助回帖

★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★ ★
感謝參與，應(yīng)助指數(shù) +1
楚天湘水: 金幣+20, ★有幫助 2012-06-02 23:38:27

1stOpt計(jì)算結(jié)果：

贊一下

回復(fù)此樓

2樓2012-05-02 10:22:42

已閱回復(fù)此樓關(guān)注TA 給TA發(fā)消息送TA紅花 TA的回帖

楚天湘水

金蟲(chóng) (小有名氣)

應(yīng)助: 1 (幼兒園)
金幣: 191.5
散金: 676
紅花: 2
帖子: 123
在線: 131.2小時(shí)
蟲(chóng)號(hào): 1167266
注冊(cè): 2010-12-10
專業(yè): 化工系統(tǒng)工程

引用回帖:

2樓: Originally posted by dingd at 2012-05-02 10:22:42:
1stOpt計(jì)算結(jié)果：
ed/e1/291104_1335925306_702.jpg

能不能把代碼貼出來(lái)啊，

贊一下

回復(fù)此樓

3樓2012-05-02 10:42:50

已閱回復(fù)此樓關(guān)注TA 給TA發(fā)消息送TA紅花 TA的回帖

hitzhengwei

木蟲(chóng) (著名寫(xiě)手)

應(yīng)助: 34 (小學(xué)生)
金幣: 1530.2
散金: 100
紅花: 1
帖子: 1921
在線: 76.5小時(shí)
蟲(chóng)號(hào): 1730549
注冊(cè): 2012-04-01
性別: GG
專業(yè): 爆炸與沖擊動(dòng)力學(xué)

【答案】應(yīng)助回帖

感謝參與，應(yīng)助指數(shù) +1

樓主的擬合代碼最好貼出來(lái)，不然怎么知道錯(cuò)誤；另外，其實(shí)可以輸出散點(diǎn)數(shù)據(jù)，然后用Origin來(lái)擬合曲線，感覺(jué)比matlab的好看很多。

贊一下

回復(fù)此樓

品味人生棋局，聆聽(tīng)世界心靈

4樓2012-05-02 10:57:15

已閱回復(fù)此樓關(guān)注TA 給TA發(fā)消息送TA紅花 TA的回帖

dingd

鐵桿木蟲(chóng) (職業(yè)作家)

計(jì)算強(qiáng)帖: 4
應(yīng)助: 1641 (講師)
金幣: 15037.3
散金: 101
紅花: 234
帖子: 3410
在線: 1223.7小時(shí)
蟲(chóng)號(hào): 291104
注冊(cè): 2006-10-28

【答案】應(yīng)助回帖

★ ★
dbb627: 金幣+2, 謝謝應(yīng)助！ 2012-05-02 14:06:30

見(jiàn)下：

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

贊一下

回復(fù)此樓

5樓2012-05-02 10:57:30

已閱回復(fù)此樓關(guān)注TA 給TA發(fā)消息送TA紅花 TA的回帖

zt_chem

木蟲(chóng) (正式寫(xiě)手)

應(yīng)助: 12 (小學(xué)生)
金幣: 2147.8
紅花: 6
帖子: 872
在線: 981.6小時(shí)
蟲(chóng)號(hào): 286920
注冊(cè): 2006-10-15
性別: GG
專業(yè): 催化化學(xué)

請(qǐng)教樓主：在用fmincon優(yōu)化時(shí)，k初始值確定為k0 = [900000 50000]有何根據(jù)？是根據(jù)文獻(xiàn)還是自己猜的值？

贊一下

回復(fù)此樓

6樓2012-05-02 21:49:24

已閱回復(fù)此樓關(guān)注TA 給TA發(fā)消息送TA紅花 TA的回帖

楚天湘水

金蟲(chóng) (小有名氣)

應(yīng)助: 1 (幼兒園)
金幣: 191.5
散金: 676
紅花: 2
帖子: 123
在線: 131.2小時(shí)
蟲(chóng)號(hào): 1167266
注冊(cè): 2010-12-10
專業(yè): 化工系統(tǒng)工程

引用回帖:

6樓: Originally posted by zt_chem at 2012-05-02 21:49:24:
請(qǐng)教樓主：在用fmincon優(yōu)化時(shí)，k初始值確定為k0 = 有何根據(jù)？是根據(jù)文獻(xiàn)還是自己猜的值？

文獻(xiàn)值，應(yīng)該就在附近，但我換初值也不行啊，都無(wú)語(yǔ)了，幫忙看一看吧

贊一下

回復(fù)此樓

7樓2012-05-03 00:10:49

已閱回復(fù)此樓關(guān)注TA 給TA發(fā)消息送TA紅花 TA的回帖

arielroy

銅蟲(chóng) (正式寫(xiě)手)

應(yīng)助: 3 (幼兒園)
金幣: 263
紅花: 8
帖子: 967
在線: 107.7小時(shí)
蟲(chóng)號(hào): 599631
注冊(cè): 2008-09-11
性別: MM
專業(yè): 催化化學(xué)

很感興趣，幫樓主頂下帖。

贊一下

回復(fù)此樓

如滿月，射天狼

8樓2012-05-11 10:40:39

已閱回復(fù)此樓關(guān)注TA 給TA發(fā)消息送TA紅花 TA的回帖

feng198800

金蟲(chóng) (正式寫(xiě)手)

應(yīng)助: 15 (小學(xué)生)
金幣: 860.2
散金: 234
紅花: 1
帖子: 415
在線: 242.1小時(shí)
蟲(chóng)號(hào): 700353
注冊(cè): 2009-02-12
性別: GG
專業(yè): 化學(xué)反應(yīng)工程

【答案】應(yīng)助回帖

請(qǐng)問(wèn)樓主活度系數(shù)：KineticsData1
% 動(dòng)力學(xué)數(shù)據(jù)
%  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 ]

是怎么回事？程序語(yǔ)句中沒(méi)有體現(xiàn)你輸入實(shí)驗(yàn)數(shù)據(jù)！

贊一下

回復(fù)此樓

9樓2012-05-29 10:46:57

已閱回復(fù)此樓關(guān)注TA 給TA發(fā)消息送TA紅花 TA的回帖

最具人氣熱帖推薦 [查看全部]	作者	回/看	最后發(fā)表

[考研] 一志愿華東理工大學(xué)，080500學(xué)碩，317分，求調(diào)劑 +9	s1145 2026-03-31	9/450	2026-03-31 18:47 by JourneyLucky
[考研] 285求調(diào)劑 +3	FZAC123 2026-03-30	3/150	2026-03-31 17:49 by 熱情沙漠
[考研] 085601英二數(shù)二求調(diào)劑總分325 +4	余航航 2026-03-31	4/200	2026-03-31 17:38 by 唐沐兒
[考研] 求調(diào)劑 +9	akdhjs 2026-03-31	10/500	2026-03-31 17:14 by zhyzzh
[考研] 315求調(diào)劑 +6	akie... 2026-03-28	7/350	2026-03-31 16:48 by asdfzly
[考研] 材料科學(xué)與工程調(diào)劑 +13	深V宿舍吧 2026-03-30	14/700	2026-03-31 16:28 by hypershenger
[考研] 本科211安全工程，初試290分，求調(diào)劑 +3	2719846834 2026-03-28	3/150	2026-03-31 13:52 by 熱情沙漠
[考研] 320分，材料與化工專業(yè)，求調(diào)劑 +10	一定上岸aaa 2026-03-27	14/700	2026-03-31 10:47 by foria
[考研] 南京大學(xué)化學(xué)調(diào)劑 +11	景隨風(fēng) 2026-03-29	16/800	2026-03-31 10:14 by herarysara
[考研] 310求調(diào)劑 +10	爭(zhēng)取九點(diǎn)睡 2026-03-30	10/500	2026-03-30 16:45 by ztnimte
[考研] 071010 323 分求調(diào)劑 +3	Baekzhy 2026-03-27	3/150	2026-03-30 14:24 by andresqi
[考研] 085701求調(diào)劑初試286分 +5	secret0328 2026-03-28	5/250	2026-03-30 12:54 by fangnagu
[考研] 086502化學(xué)工程342求調(diào)劑 +6	阿姨復(fù)古不過(guò) 2026-03-27	6/300	2026-03-28 07:06 by wangy0907
[考研] 265求調(diào)劑 +8	小木蟲(chóng)085600 2026-03-27	8/400	2026-03-27 22:16 by 無(wú)際的草原
[考研] 一志愿上海理工能源動(dòng)力（085800）310分求調(diào)劑 +3	zhangmingc 2026-03-27	4/200	2026-03-27 19:01 by 給你你注意休息
[考研] 求調(diào)劑 +4	零八# 2026-03-27	4/200	2026-03-27 18:07 by yu221
[考研] 348求調(diào)劑 +4	小懶蟲(chóng)不懶了 2026-03-27	5/250	2026-03-27 12:47 by 果果媽咪
[考研] 一志愿鄭大085600，310分求調(diào)劑 +5	李瀟可 2026-03-26	5/250	2026-03-27 11:14 by 不吃魚(yú)的貓
[考研] 【2026考研調(diào)劑】制藥工程 284分求相關(guān)專業(yè)調(diào)劑名額 +4	袁奐奐 2026-03-25	8/400	2026-03-25 14:32 by lbsjt
[考研] 各位老師您好：本人初試372分 +5	jj涌77 2026-03-25	6/300	2026-03-25 14:15 by mapenggao