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

zhangqq72270

木蟲 (正式寫手)

應(yīng)助: 2 (幼兒園)
金幣: 2102.1
帖子: 530
在線: 202.4小時(shí)
蟲號(hào): 1994269

[交流] 使用lsqnonlin函數(shù)優(yōu)化動(dòng)力學(xué)參數(shù)，總是得不到合理的結(jié)果

大家好:
      我最近在使用matlab回歸動(dòng)力學(xué)參數(shù)，首先是使用4階龍哥庫塔算法進(jìn)行積分，之后使用lsqnonlin進(jìn)行非線性最小二乘擬合，運(yùn)行已經(jīng)寫好的m文件，總是得不到合理的結(jié)果，不知道是不是自己程序編寫有問題，請(qǐng)大神指導(dǎo)，謝謝。
下面粘貼上自己的m文件：
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function zqq100
% 動(dòng)力學(xué)ODE方程模型的參數(shù)估計(jì)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
clear all
clc
kz0 = [1.73  -3  4.3  3.4  29.4  26.46  22.8  29.4];       % 參數(shù)初值
lb  = [];                                                 % 參數(shù)下限
ub  = [];                                                 % 參數(shù)上限
x0  = [72516.987800 7922.540990 2213.519110 2440.339360];
KineticsData;
yexp = ExpData(:,2:5);                   % yexp: 實(shí)驗(yàn)數(shù)據(jù)[x1 x2 x3 x4]

%再參數(shù)化
a = 0; b = 0.022146; d=0.000022146;
tt=a:d:b;
z = func(tt);
Tintegral = quadl(@func,a,b);             %  實(shí)際上應(yīng)該是(1000/T)integral
Taverage  = Tintegral/(b-a) ;             %  實(shí)際上應(yīng)該是(1000/T)average

% 優(yōu)化參數(shù)設(shè)置
options=optimset ('Display','iter','TolX',1e-25,'TolFun',1e-25,'LargeScale','off','MaxIter',2,'MaxFunEvals',inf);

% 使用函數(shù)lsqnonlin()進(jìn)行參數(shù)估計(jì)
[kz,resnorm,residual,exitflag,output,lambda,jacobian] = ...
lsqnonlin(@ObjFunc4LNL,kz0,lb,ub,options,x0,yexp,Taverage);
ci = nlparci(kz,residual,jacobian);

%把回歸參數(shù)轉(zhuǎn)化為自己需要的參數(shù)
k01 = exp(kz(1)+kz(5)*Taverage);
ke1 = exp(ci(1,2)-kz(1)+(ci(5,2)-kz(5))*Taverage);
Ea1 = kz(5)*1000*8.314;
Ee1 = (ci(5,2)-kz(5))*1000*8.314;

k02 = exp(kz(2)+kz(6)*Taverage);
ke2 = exp(ci(2,2)-kz(2)+(ci(6,2)-kz(6))*Taverage);
Ea2 = kz(6)*1000*8.314;
Ee2 = (ci(6,2)-kz(6))*1000*8.314;

k03 = exp(kz(3)+kz(7)*Taverage);
ke3 = exp(ci(3,2)-kz(3)+(ci(7,2)-kz(7))*Taverage);
Ea3 = kz(7)*1000*8.314;
Ee3 = (ci(7,2)-kz(7))*1000*8.314;

k04 = exp(kz(4)+kz(8)*Taverage);
ke4 = exp(ci(4,2)-kz(4)+(ci(8,2)-kz(8))*Taverage);
Ea4 = kz(8)*1000*8.314;
Ee4 = (ci(8,2)-kz(8))*1000*8.314;
% ------------------------------------------------------------------
% 輸出結(jié)果
fprintf('\n\n使用函數(shù)lsqnonlin()估計(jì)得到的參數(shù)值為:\n')
fprintf('\tk01 = %.4f ± %.4f\n',k01,ke1)
fprintf('\tk02 = %.4f ± %.4f\n',Ea1,Ee1)
fprintf('\tk03 = %.4f ± %.4f\n',k02,ke2)
fprintf('\tk04 = %.4f ± %.4f\n',Ea2,Ee2)
fprintf('\tEa1 = %.4f ± %.4f\n',k03,ke3)
fprintf('\tEa2 = %.4f ± %.4f\n',Ea3,Ee3)
fprintf('\tEa3 = %.4f ± %.4f\n',k04,ke4)
fprintf('\tEa4 = %.4f ± %.4f\n',Ea4,Ee4)
fprintf('  The sum of the squares is: %.1e\n\n',resnorm)

% ------------------------------------------------------------------
function f = ObjFunc4LNL(kz,x0,yexp,Taverage)
tspan = [0.000000 0.006021 0.009538 0.012340 0.014647 0.016594 0.018268 0.019730 0.020395 0.021021 0.021613 0.022146];
[t x] = ode45(@KineticEqs,tspan,x0,[],kz,Taverage);
y(:,1:4)=x(:,1:4);
f1 = y(:,1) - yexp(:,1);
f2 = y(:,2) - yexp(:,2);
f3 = y(:,3) - yexp(:,3);
f4 = y(:,4) - yexp(:,4);
f = [f1; f2; f3; f4];

% ------------------------------------------------------------------
function dxdt = KineticEqs(t,x,kz,Taverage)

T = 2E+09*t.^4 - 8E+07*t.^3 + 1E+06*t.^2 + 1593.6*t + 823.94;
D = 1000./T-Taverage;

k(1)=exp(kz(1)-kz(5)*D);
k(2)=exp(kz(2)-kz(6)*D);
k(3)=exp(kz(3)-kz(7)*D);
k(4)=exp(kz(4)-kz(8)*D);

dxdt =  ...
[ ( - k(1)*x(1) )
( - k(1)*x(1) - k(3)*x(3) + k(2)*x(2) )
( - k(1)*x(1) + k(3)*x(3) + k(4)*x(3) )
( - k(2)*x(2) - k(3)*x(3) - k(4)*x(3) )
];
% ------------------------------------------------------------------
function z = func(tt)
z =1000./( 2E+09*tt.^4 - 8E+07*tt.^3 + 1E+06*tt.^2 + 1593.6*tt + 823.94);

數(shù)據(jù)文件
% 動(dòng)力學(xué)數(shù)據(jù):
%    t       y1       y2       y3       y4
ExpData = ...
[ 0.000000 72516.987800 7922.540990       2213.519110    2440.339360
0.006021 42499.607900 22206.384100 3306.398400    3606.838600
0.009538 29000.355000 29930.619400 3600.077060    3836.127990
0.012340 20185.280200 35570.546400 3528.468400    3684.127360
0.014647 14159.635500 39715.404600 3262.616960    3307.931340
0.016594 9923.223120       42747.821100 2919.806950    2824.661990
0.018268 6889.657400       44873.256100 2561.423650    2313.346410
0.019730 4712.646380       46269.306400 2219.932610    1823.752130
0.020395 3871.540620       46727.452000 2061.574630    1592.788930
0.021021 3166.514550       47033.087200 1914.970250    1369.417740
0.021613 2580.686780       47184.935100 1784.008050    1149.425230
0.022146 2121.809490       47180.789300 1679.141140    937.562404
]

[ Last edited by zhangqq72270 on 2013-8-14 at 10:15 ]

回復(fù)此樓

» 猜你喜歡

309求調(diào)劑已經(jīng)有8人回復(fù)
一志愿南航 335分 | 0856 | GPA 4.07 | 有科研經(jīng)歷已經(jīng)有8人回復(fù)
求調(diào)劑已經(jīng)有5人回復(fù)
求化學(xué)調(diào)劑已經(jīng)有5人回復(fù)
085600，材料與化工321分求調(diào)劑已經(jīng)有10人回復(fù)
0703 化學(xué) 求調(diào)劑，一志愿山東大學(xué) 342 分已經(jīng)有5人回復(fù)
一志愿南開大學(xué)0710生物學(xué)359求調(diào)劑已經(jīng)有3人回復(fù)
085600，專業(yè)課化工原理，320分求調(diào)劑已經(jīng)有4人回復(fù)
材料與化工272求調(diào)劑已經(jīng)有16人回復(fù)
290求調(diào)劑已經(jīng)有3人回復(fù)

» 本主題相關(guān)價(jià)值貼推薦，對(duì)您同樣有幫助:

matlab的lsqnonlin函數(shù)怎么用已經(jīng)有7人回復(fù)
緊急求助，利用Matlab對(duì)實(shí)驗(yàn)數(shù)據(jù)進(jìn)行擬合求解參數(shù)。已經(jīng)有27人回復(fù)
關(guān)于matlab的參數(shù)估計(jì) 已經(jīng)有15人回復(fù)
怎樣用一組參數(shù)同時(shí)擬合兩個(gè)曲線--matlab 已經(jīng)有5人回復(fù)
Matlab最小二乘參數(shù)優(yōu)化已經(jīng)有3人回復(fù)
幫忙給看看這個(gè)matlab優(yōu)化函數(shù) 問題已經(jīng)有8人回復(fù)
動(dòng)力學(xué)參數(shù)擬合已經(jīng)有26人回復(fù)
MATLAB非線性優(yōu)化擬合怎么改才正確已經(jīng)有3人回復(fù)
matlab fmincon優(yōu)化函數(shù) 怎樣知道每次得到的迭代點(diǎn) 已經(jīng)有3人回復(fù)
求助用Oringin擬合轉(zhuǎn)狀態(tài)方程已經(jīng)有6人回復(fù)
請(qǐng)教matlab反應(yīng)動(dòng)力學(xué)參數(shù)估計(jì)遇到的問題，謝謝已經(jīng)有15人回復(fù)
擬合非線性方程系數(shù)的精度問題已經(jīng)有5人回復(fù)
matlab擬合方程參數(shù)時(shí)初值的選擇已經(jīng)有15人回復(fù)
matlab擬合模型參數(shù) 已經(jīng)有5人回復(fù)
lsqnonlin函數(shù)擬合微分方程組參數(shù)擬合問題已經(jīng)有10人回復(fù)
matlab非線性參數(shù)擬合問題已經(jīng)有7人回復(fù)
【求助】催化反應(yīng)動(dòng)力學(xué)matlab計(jì)算各基元反應(yīng)的速率常數(shù)時(shí)，該如何避免較小量被忽略？已經(jīng)有3人回復(fù)
【求助】催化反應(yīng)動(dòng)力學(xué) 已經(jīng)有5人回復(fù)
【求助】matlab 二次規(guī)劃的優(yōu)化的問題已經(jīng)有4人回復(fù)
【求助】使用Matlab擬合反應(yīng)動(dòng)力學(xué)方程問題已經(jīng)有7人回復(fù)
【求助】使用Matlab預(yù)估動(dòng)力學(xué)方程問題已經(jīng)有13人回復(fù)

» 搶金幣啦！回帖就可以得到:

查看全部散金貼

1樓 2013-08-14 09:53:53

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

sunshine819

金蟲 (小有名氣)

應(yīng)助: 32 (小學(xué)生)
金幣: 1170.9
帖子: 137
在線: 42.6小時(shí)
蟲號(hào): 2254567

★
小木蟲: 金幣+0.5, 給個(gè)紅包，謝謝回帖

強(qiáng)非線性的是不是需要用其他算法，比如GA這類全局搜索的，這種問題如果給的初值不在平衡點(diǎn)的吸引域內(nèi)是不是就沒法獲得很好的效果。

贊一下

回復(fù)此樓

» 本帖已獲得的紅花（最新10朵）

zhangqq72270

4樓2013-08-16 11:38:32

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

查看全部 7 個(gè)回答

nanaxiao86

木蟲 (著名寫手)

應(yīng)助: 31 (小學(xué)生)
金幣: 1561.7
帖子: 1555
在線: 134.8小時(shí)
蟲號(hào): 1448832

祝你的問題早日解決啦~

贊一下

回復(fù)此樓

2樓2013-08-14 10:22:39

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

zhangqq72270

木蟲 (正式寫手)

應(yīng)助: 2 (幼兒園)
金幣: 2102.1
帖子: 530
在線: 202.4小時(shí)
蟲號(hào): 1994269

引用回帖:

2樓: Originally posted by nanaxiao86 at 2013-08-14 10:22:39
祝你的問題早日解決啦~

謝謝

回復(fù)此樓

3樓2013-08-14 10:25:21

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

wswsjz

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

應(yīng)助: 66 (初中生)
金幣: 3265.5
帖子: 3032
在線: 506.1小時(shí)
蟲號(hào): 511383

★
小木蟲: 金幣+0.5, 給個(gè)紅包，謝謝回帖

努力，再努力！再創(chuàng)佳績！

贊一下

回復(fù)此樓

5樓2013-08-16 11:48:16

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

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

[考研] 一志愿河北工業(yè)大學(xué)0817化工278分求調(diào)劑 +10	jhybd 2026-03-23	15/750	2026-03-29 21:49 by nanaliuyun
[考研] 一志愿：西北大學(xué)，英一數(shù)一408-284分求調(diào)劑 +4	12.27 2026-03-27	4/200	2026-03-29 14:40 by zhshch
[考研] 一志愿鄭州大學(xué)，080500學(xué)碩，總分317分求調(diào)劑 +8	舉個(gè)栗子oi 2026-03-24	9/450	2026-03-29 13:08 by peike
[考研] 348求調(diào)劑 +5	小懶蟲不懶了 2026-03-28	5/250	2026-03-29 10:34 by 唐沐兒
[考研] 305求調(diào)劑 +8	RuiFairyrui 2026-03-28	8/400	2026-03-29 08:22 by fmesaito
[考研] 一志愿華理，數(shù)一英一285求A區(qū)調(diào)劑 +8	AZMK 2026-03-25	12/600	2026-03-28 18:15 by AZMK
[考研] 266分，求材料冶金能源化工等調(diào)劑 +7	哇呼哼呼哼 2026-03-27	9/450	2026-03-28 12:22 by zllcz
[考研] 求調(diào)劑推薦材料 304 +15	荷包蛋hyj 2026-03-26	15/750	2026-03-28 04:13 by fmesaito
[考研] 330一志愿中國海洋大學(xué) 化學(xué)工程 085602 有讀博意愿求調(diào)劑 +3	wywy.. 2026-03-27	4/200	2026-03-28 03:32 by fmesaito
[考研] 復(fù)試調(diào)劑，一志愿南農(nóng)083200食品科學(xué)與工程 +5	XQTJZ 2026-03-26	5/250	2026-03-27 14:49 by 狂炫麥當(dāng)當(dāng)
[考研] 085600，材料與化工321分，求調(diào)劑 +9	大饞小子 2026-03-27	9/450	2026-03-27 14:30 by mmm just
[考研] 一志愿鄭大085600，310分求調(diào)劑 +5	李瀟可 2026-03-26	5/250	2026-03-27 11:14 by 不吃魚的貓
[考研] 0703化學(xué)一志愿南京師范大學(xué)303求調(diào)劑 +3	zzffylgg 2026-03-24	3/150	2026-03-27 10:42 by shangxh
[考研] 打過很多競(jìng)賽，085406控制工程300分，求調(diào)劑 +3	askeladz 2026-03-26	3/150	2026-03-26 09:08 by 給你你注意休息
[考研] 332求調(diào)劑 +6	032500 2026-03-25	6/300	2026-03-25 22:45 by 418490947
[考研] 各位老師您好：本人初試372分 +5	jj涌77 2026-03-25	6/300	2026-03-25 14:15 by mapenggao
[考研] 285求調(diào)劑 +3	AZMK 2026-03-24	3/150	2026-03-25 12:23 by userper
[有機(jī)交流] 有機(jī)合成求助 20+3	FENGSHUJEI 2026-03-23	5/250	2026-03-24 19:31 by 88817753
[考研] 材料專碩找調(diào)劑 +5	哈哈哈吼吼吼哈 2026-03-23	5/250	2026-03-24 19:07 by 了了了了。。
[考研] 344求調(diào)劑 +3	desto 2026-03-24	3/150	2026-03-24 10:09 by 搏擊518