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ù)此樓

» 猜你喜歡

281求調(diào)劑（0805）已經(jīng)有16人回復(fù)
304求調(diào)劑已經(jīng)有6人回復(fù)
材料工程專碩調(diào)劑已經(jīng)有6人回復(fù)
一志愿天大材料與化工（085600）總分338 已經(jīng)有4人回復(fù)
085700資源與環(huán)境308求調(diào)劑已經(jīng)有3人回復(fù)
求材料調(diào)劑已經(jīng)有8人回復(fù)
294求調(diào)劑材料與化工專碩已經(jīng)有5人回復(fù)
一志愿華中科技大學(xué)，080502，354分求調(diào)劑已經(jīng)有4人回復(fù)
一志愿吉林大學(xué)材料學(xué)碩321求調(diào)劑已經(jīng)有6人回復(fù)
085410人工智能專碩317求調(diào)劑（0854都可以）已經(jīng)有3人回復(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的回帖

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的回帖

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的回帖

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的回帖

zhangqq72270

木蟲 (正式寫手)

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

送紅花一朵

引用回帖:

4樓: Originally posted by sunshine819 at 2013-08-16 11:38:32
強非線性的是不是需要用其他算法，比如GA這類全局搜索的，這種問題如果給的初值不在平衡點的吸引域內(nèi)是不是就沒法獲得很好的效果。

現(xiàn)在就是很難找到合適的初始值，不知道你有沒有什么好的辦法

贊一下

回復(fù)此樓

6樓2013-08-16 14:49:19

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

zhangqq72270

木蟲 (正式寫手)

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

引用回帖:

5樓: Originally posted by wswsjz at 2013-08-16 11:48:16
努力，再努力！再創(chuàng)佳績！

一定

回復(fù)此樓

7樓2013-08-16 14:49:43

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

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

[考研] 能源材料化學(xué)課題組招收碩士研究生8-10名 +4	脫穎而出 2026-03-16	10/500	2026-03-18 20:59 by 脫穎而出
[考研] 材料專碩英一數(shù)二306 +4	z1z2z3879 2026-03-18	4/200	2026-03-18 20:00 by 楤哥
[考研] 26調(diào)劑/材料/英一數(shù)二/總分289/已過A區(qū)線 +7	步川酷紫123 2026-03-13	7/350	2026-03-18 17:12 by 盡舜堯1
[考研] 298-一志愿中國農(nóng)業(yè)大學(xué)-求調(diào)劑 +7	手機用戶 2026-03-17	7/350	2026-03-18 14:34 by vgtyfty
[考研] 工科材料085601 279求調(diào)劑 +6	困于星晨 2026-03-17	6/300	2026-03-18 10:21 by kkcoco25
[考研] 265求調(diào)劑 +3	梁梁校校 2026-03-17	3/150	2026-03-18 09:12 by zhukairuo
[考研] 268求調(diào)劑 +6	簡單點0 2026-03-17	6/300	2026-03-18 09:04 by 無際的草原
[碩博家園] 湖北工業(yè)大學(xué) 生命科學(xué)與健康學(xué)院-課題組招收2026級食品/生物方向碩士 +3	1喜春8 2026-03-17	5/250	2026-03-17 17:18 by ber川cool子
[考研] 275求調(diào)劑 +4	太陽花天天開心 2026-03-16	4/200	2026-03-17 10:53 by 功夫瘋狂
[考研] 302求調(diào)劑 +4	小賈同學(xué)123 2026-03-15	8/400	2026-03-17 10:33 by 小賈同學(xué)123
[考研] [導(dǎo)師推薦]西南科技大學(xué)國防/材料導(dǎo)師推薦 +3	尖角小荷 2026-03-16	6/300	2026-03-16 23:21 by 尖角小荷
[考研] 304求調(diào)劑 +5	素年祭語 2026-03-15	5/250	2026-03-16 17:00 by 我的船我的海
[考研] 一志愿211 0703方向310分求調(diào)劑 +3	努力奮斗112 2026-03-15	3/150	2026-03-16 16:44 by houyaoxu
[考研] 0856求調(diào)劑 +3	劉夢微 2026-03-15	3/150	2026-03-16 10:00 by houyaoxu
[考研] 327求調(diào)劑 +6	拾光任染 2026-03-15	11/550	2026-03-15 22:47 by 拾光任染
[考研] 308 085701 四六級已過求調(diào)劑 +7	溫喬喬喬喬 2026-03-12	14/700	2026-03-14 10:49 by JourneyLucky
[考研] 281求調(diào)劑 +9	Koxui 2026-03-12	11/550	2026-03-13 20:50 by Koxui
[考研] 26調(diào)劑/材料科學(xué)與工程/總分295/求收留 +9	2026調(diào)劑俠 2026-03-12	9/450	2026-03-13 20:46 by 18595523086
[考研] 考研調(diào)劑 +4	芬達46 2026-03-12	4/200	2026-03-13 16:04 by ruiyingmiao
[考研] 290求調(diào)劑 +7	ADT 2026-03-12	7/350	2026-03-13 15:17 by JourneyLucky