zhangqq72270

木蟲 (正式寫手)

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

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

大家好:
      我最近在使用matlab回歸動力學(xué)參數(shù)，首先是使用4階龍哥庫塔算法進行積分，之后使用lsqnonlin進行非線性最小二乘擬合，運行已經(jīng)寫好的m文件，總是得不到合理的結(jié)果，不知道是不是自己程序編寫有問題，請大神指導(dǎo)，謝謝。
下面粘貼上自己的m文件：
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function zqq100
% 動力學(xué)ODE方程模型的參數(shù)估計
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
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ù)據(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);             %  實際上應(yīng)該是(1000/T)integral
Taverage  = Tintegral/(b-a) ;             %  實際上應(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()進行參數(shù)估計
[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()估計得到的參數(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ù)文件
% 動力學(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ù)此樓

» 猜你喜歡

考研調(diào)劑已經(jīng)有3人回復(fù)
317分一志愿南理工材料工程本科湖工大求調(diào)劑已經(jīng)有9人回復(fù)
297求調(diào)劑已經(jīng)有10人回復(fù)
一志愿武漢理工，總分321，英一數(shù)二，求老師收留。已經(jīng)有10人回復(fù)
343求調(diào)劑已經(jīng)有5人回復(fù)
085601材料工程找調(diào)劑已經(jīng)有7人回復(fù)
求調(diào)劑，一志愿南京航空航天大學(xué) ，080500材料科學(xué)與工程學(xué)碩，總分289分已經(jīng)有6人回復(fù)
一志愿南航 335分 | 0856 | GPA 4.07 | 有科研經(jīng)歷已經(jīng)有5人回復(fù)
340求調(diào)劑已經(jīng)有5人回復(fù)
一志愿北化求調(diào)劑已經(jīng)有4人回復(fù)

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

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

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

查看全部散金貼

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

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

zhangqq72270

木蟲 (正式寫手)

應(yīng)助: 2 (幼兒園)
金幣: 2102.1
帖子: 530
在線: 202.4小時
蟲號: 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的回帖

查看全部 7 個回答

nanaxiao86

木蟲 (著名寫手)

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

祝你的問題早日解決啦~

贊一下

回復(fù)此樓

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

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

sunshine819

金蟲 (小有名氣)

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

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

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

贊一下

回復(fù)此樓

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

zhangqq72270

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

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

wswsjz

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

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

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

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

贊一下

回復(fù)此樓

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

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

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

[考研] 考研調(diào)劑 +3	小蠟新筆 2026-03-29	3/150	2026-03-29 11:47 by JourneyLucky
[考研] 一志愿雙一流機械285分求調(diào)劑 +3	幸運的三木 2026-03-29	4/200	2026-03-29 10:46 by 想上岸了嗚嗚
[考研] 一志愿鄭州大學(xué)，080500學(xué)碩，總分317分求調(diào)劑 +6	舉個栗子oi 2026-03-24	7/350	2026-03-29 10:39 by 唐沐兒
[考研] 283求調(diào)劑 +3	A child 2026-03-28	3/150	2026-03-28 15:41 by ms629
[考研] 394求調(diào)劑 +3	好事多磨靜候佳?/a> 2026-03-26	5/250	2026-03-28 14:24 by 唐沐兒
[考研] 299求調(diào)劑 +7	嗯嗯嗯嗯2 2026-03-27	7/350	2026-03-28 13:09 by 唐沐兒
[考研] 085701環(huán)境工程，267求調(diào)劑 +16	minht 2026-03-26	16/800	2026-03-28 12:16 by zllcz
[考研] 311求調(diào)劑 +3	希望上岸阿小楊 2026-03-23	3/150	2026-03-28 07:57 by 熱情沙漠
[考研] 085404求調(diào)劑，總分309，本科經(jīng)歷較為豐富 +4	來財aa 2026-03-25	4/200	2026-03-28 07:41 by 棒棒球手
[考研] 張芳銘-中國農(nóng)業(yè)大學(xué)-環(huán)境工程專碩-298 +4	手機用戶 2026-03-26	4/200	2026-03-28 07:17 by mmm just
[考研] 07化學(xué)280分求調(diào)劑 +10	722865 2026-03-23	10/500	2026-03-27 15:51 by Plutoqq
[考研] 復(fù)試調(diào)劑，一志愿南農(nóng)083200食品科學(xué)與工程 +5	XQTJZ 2026-03-26	5/250	2026-03-27 14:49 by 狂炫麥當(dāng)當(dāng)
[考研] 321求調(diào)劑 +6	Ymlll 2026-03-24	6/300	2026-03-26 20:50 by 不吃魚的貓
[考研] 化學(xué)工程085602 305分求調(diào)劑 +17	RichLi_ 2026-03-25	17/850	2026-03-26 19:44 by plmuchong
[考研] 一志愿南京郵電大學(xué) 288分材料考研求調(diào)劑 +3	jl0720 2026-03-26	3/150	2026-03-26 13:39 by zzll406
[考研] 生物技術(shù)與工程 +3	1294608413 2026-03-25	4/200	2026-03-25 18:02 by 1294608413
[考研] 一志愿北化315 求調(diào)劑 +3	akrrain 2026-03-24	3/150	2026-03-24 19:35 by 了了了了。。
[考研] 化工專碩求調(diào)劑 +3	question挽風(fēng) 2026-03-24	3/150	2026-03-24 18:48 by jhhcooi
[考研] 292求調(diào)劑 +4	鵝鵝鵝額額額額?/a> 2026-03-24	4/200	2026-03-24 16:41 by peike
[考研] 335求調(diào)劑 +4	yuyu宇 2026-03-23	5/250	2026-03-23 23:49 by Txy@872106