楚天湘水

應(yīng)助: 1 (幼兒園)
金幣: 191.5
散金: 676
紅花: 2
帖子: 123
在線(xiàn): 131.2小時(shí)
蟲(chóng)號(hào): 1167266
注冊(cè): 2010-12-10
專(zhuān)業(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
在線(xiàn): 242.1小時(shí)
蟲(chóng)號(hào): 700353
注冊(cè): 2009-02-12
性別: GG
專(zhuān)業(yè): 化學(xué)反應(yīng)工程

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

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

哦，看明白了！呵呵。最小二乘法可以直接擬合，你運(yùn)用非線(xià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
在線(xiàn): 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
在線(xiàn): 131.2小時(shí)
蟲(chóng)號(hào): 1167266
注冊(cè): 2010-12-10
專(zhuān)業(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
在線(xiàn): 76.5小時(shí)
蟲(chóng)號(hào): 1730549
注冊(cè): 2012-04-01
性別: GG
專(zhuān)業(yè): 爆炸與沖擊動(dòng)力學(xué)

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

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

樓主的擬合代碼最好貼出來(lái)，不然怎么知道錯(cuò)誤；另外，其實(shí)可以輸出散點(diǎn)數(shù)據(jù)，然后用Origin來(lái)擬合曲線(xiàn)，感覺(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
在線(xiàn): 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
在線(xiàn): 981.6小時(shí)
蟲(chóng)號(hào): 286920
注冊(cè): 2006-10-15
性別: GG
專(zhuān)業(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
在線(xiàn): 131.2小時(shí)
蟲(chóng)號(hào): 1167266
注冊(cè): 2010-12-10
專(zhuān)業(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
在線(xiàn): 107.7小時(shí)
蟲(chóng)號(hào): 599631
注冊(cè): 2008-09-11
性別: MM
專(zhuān)業(yè): 催化化學(xué)

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

贊一下

回復(fù)此樓

如滿(mǎn)月，射天狼

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
在線(xiàn): 242.1小時(shí)
蟲(chóng)號(hào): 700353
注冊(cè): 2009-02-12
性別: GG
專(zhuān)業(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ā)表

[考研] 085600，321分求調(diào)劑 +3	大饞小子 2026-03-31	3/150	2026-03-31 14:43 by syjjj0321
[考研] 290分調(diào)劑求助 +9	吉祥止止陳 2026-03-25	9/450	2026-03-31 14:16 by 六翼TMR
[考研] 334求調(diào)劑 +7	Trying] 2026-03-31	7/350	2026-03-31 12:33 by 無(wú)際的草原
[考研] 一志愿大連理工大學(xué)，機(jī)械工程學(xué)碩，341 +3	西瓜田的守望者 2026-03-30	3/150	2026-03-31 11:08 by asdfzly
[考研] 320分，材料與化工專(zhuān)業(yè)，求調(diào)劑 +10	一定上岸aaa 2026-03-27	14/700	2026-03-31 10:47 by foria
[考研] 085601一志愿中山大學(xué)深圳材料工程330求調(diào)劑 +5	pipiver 2026-03-30	5/250	2026-03-31 07:37 by JourneyLucky
[考研] 材料求調(diào)劑一志愿哈工大總分298分，前三科223分 +10	dongfang59 2026-03-27	10/500	2026-03-30 23:42 by 果果媽咪
[考研] 0703化學(xué)321分求調(diào)劑 +10	三dd. 2026-03-30	11/550	2026-03-30 19:24 by markhwc
[考研] 289求調(diào)劑 +16	新時(shí)代材料 2026-03-27	16/800	2026-03-30 19:04 by Wang200018
[考研] 0703 化學(xué) 求調(diào)劑，一志愿山東大學(xué) 342 分 +7	Shern—- 2026-03-28	7/350	2026-03-30 16:31 by nothing投稿中
[考研] 0703化學(xué)/290求調(diào)劑/本科經(jīng)歷豐富/工科也可 +13	丹青奶蓋 2026-03-26	15/750	2026-03-30 12:35 by fangnagu
[考研] 求調(diào)劑，一志愿南京航空航天大學(xué)大學(xué) ，080500材料科學(xué)與工程學(xué)碩 +6	@taotao 2026-03-26	7/350	2026-03-30 10:43 by 我是小康
[考研] 349求調(diào)劑 +6	李木子啊哈哈 2026-03-25	6/300	2026-03-29 12:47 by 無(wú)際的草原
[考研] 數(shù)一英一271專(zhuān)碩（085401）求調(diào)劑，可跨 +7	前行必有光 2026-03-28	8/400	2026-03-28 23:22 by 小木蟲(chóng)tim
[考研] 081200-314 +3	LILIQQ 2026-03-27	4/200	2026-03-28 09:41 by 保護(hù)地球你我做?/a>
[考研] 070300化學(xué)求調(diào)劑 +4	起個(gè)名咋這么難 2026-03-27	4/200	2026-03-27 21:39 by 83503孫老師
[考研] 求調(diào)劑 +4	零八# 2026-03-27	4/200	2026-03-27 18:07 by yu221
[考研] 292求調(diào)劑 +4	求求了收下我吧?/a> 2026-03-26	4/200	2026-03-27 10:37 by zhshch
[考研] 環(huán)境專(zhuān)碩324分求調(diào)劑推薦 +5	軒小寧—— 2026-03-26	5/250	2026-03-26 12:05 by i_cooler
[考研] 一志愿天津大學(xué)339材料與化工求調(diào)劑 +3	江往賣(mài)魚(yú) 2026-03-26	3/150	2026-03-26 09:42 by 王小欠i