當(dāng)前只顯示滿足指定條件的回帖，點(diǎn)擊這里查看本話題的所有回帖

楚天湘水

金蟲 (小有名氣)

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

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

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

下面是寫的matlab程序：
function Kinetic2
% 動力學(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])

%動力學(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ī)兔匆豢�，哪里有問題，從做出的圖可以看出，處理的結(jié)果很差啊。
無論K初值給出怎么樣，都是擬合出來的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的回帖

hitzhengwei

木蟲 (著名寫手)

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

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

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

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

贊一下

回復(fù)此樓

品味人生棋局，聆聽世界心靈

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

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

查看全部 15 個回答

dingd

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

計(jì)算強(qiáng)帖: 4
應(yīng)助: 1641 (講師)
金幣: 15037.3
散金: 101
紅花: 234
帖子: 3410
在線: 1223.7小時
蟲號: 291104
注冊: 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的回帖

楚天湘水

金蟲 (小有名氣)

應(yīng)助: 1 (幼兒園)
金幣: 191.5
散金: 676
紅花: 2
帖子: 123
在線: 131.2小時
蟲號: 1167266
注冊: 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

能不能把代碼貼出來啊，

贊一下

回復(fù)此樓

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

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

dingd

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

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

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

★ ★
dbb627: 金幣+2, 謝謝應(yīng)助！ 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

贊一下

回復(fù)此樓

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

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

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

[考研] 080200專業(yè)277分，求帶走！ +3	瓶子PZ 2026-03-31	3/150	2026-03-31 18:19 by 夜空空哦
[考研] 一志愿華東師范大學(xué)有機(jī)化學(xué)專業(yè)，初試351分，復(fù)試被刷求調(diào)劑! +9	真名有冰 2026-03-29	10/500	2026-03-31 18:01 by xhai2011
[考研] 環(huán)境工程297分求調(diào)劑一志愿杭高院 +6	GENJIOW 2026-03-31	6/300	2026-03-31 18:01 by 544594351
[考研] 290分調(diào)劑求助 +10	吉祥止止陳 2026-03-25	10/500	2026-03-31 17:54 by 544594351
[考研] 生醫(yī)工0831調(diào)劑求推薦 +4	小熊睿睿_s 2026-03-27	6/300	2026-03-31 17:12 by 記事本2026
[考研] 086000生物與醫(yī)藥初試274求調(diào)劑 +4	小叮當(dāng)來了 2026-03-30	4/200	2026-03-31 16:48 by shengliu165
[考研] 362求調(diào)劑 +9	西南交材料專碩3 2026-03-31	9/450	2026-03-31 16:25 by w蟲蟲123
[考研] 0703化學(xué) +20	妮妮ninicgb 2026-03-27	20/1000	2026-03-31 13:33 by 無際的草原
[考研] 081200-11408-276學(xué)碩求調(diào)劑 +4	崔wj 2026-03-31	4/200	2026-03-31 11:56 by jp9609
[考研] 269求調(diào)劑 +4	我想讀研11 2026-03-31	4/200	2026-03-31 10:04 by cal0306
[考研] 一志愿食品科學(xué)與工程083200求調(diào)劑 +4	XQTJZ 2026-03-30	4/200	2026-03-31 04:10 by fmesaito
[考研] 福建理工大學(xué)材料學(xué)院先進(jìn)合金團(tuán)隊(duì)招收考研調(diào)劑學(xué)生 +3	大華金商都 2026-03-30	4/200	2026-03-31 01:04 by 方英俊602
[考研] 食品工程專碩一志愿中海洋309求調(diào)劑 +5	小張zxy張 2026-03-26	10/500	2026-03-31 00:29 by jp9609
[考研] 一志愿北京化工大學(xué)材料與化工（085600）296求調(diào)劑 +25	稻妻小編 2026-03-26	25/1250	2026-03-30 20:11 by 滴滴上岸呀
[考研] 材料與化工304求B區(qū)調(diào)劑 +4	邱gl 2026-03-26	7/350	2026-03-30 08:39 by 探123
[考研] 275求調(diào)劑 +10	jjjjjjjjjjl 2026-03-27	10/500	2026-03-27 23:47 by barnett0632
[考研] 298調(diào)劑 +3	jiyingjie123 2026-03-27	3/150	2026-03-27 11:57 by wxiongid
[考研] 調(diào)劑推薦 +5	清酒714 2026-03-26	6/300	2026-03-27 11:12 by 不吃魚的貓
[論文投稿] Journal of Mechanical Science and Technology +3	Russ_ss 2026-03-25	5/250	2026-03-27 10:49 by 陸小果畫大餅
[考研] 調(diào)劑 +4	柚柚yoyo 2026-03-26	4/200	2026-03-26 20:43 by fmesaito