zhangqq72270

木蟲 (正式寫手)

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

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

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

回復此樓

» 猜你喜歡

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

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

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

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

查看全部散金貼

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

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

nanaxiao86

木蟲 (著名寫手)

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

祝你的問題早日解決啦~

贊一下

回復此樓

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

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

zhangqq72270

木蟲 (正式寫手)

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

引用回帖:

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

謝謝

回復此樓

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

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

sunshine819

金蟲 (小有名氣)

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

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

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

贊一下

回復此樓

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

zhangqq72270

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

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

wswsjz

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

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

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

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

贊一下

回復此樓

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

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

zhangqq72270

木蟲 (正式寫手)

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

送紅花一朵

引用回帖:

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

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

贊一下

回復此樓

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

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

zhangqq72270

木蟲 (正式寫手)

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

引用回帖:

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

一定

回復此樓

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

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

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

[考研] 304求調(diào)劑 +6	司空. 2026-03-18	6/300	2026-03-18 23:03 by 星空星月
[考研] 一志愿中海洋材料工程專碩330分求調(diào)劑 +5	小材化本科 2026-03-18	5/250	2026-03-18 21:53 by 無懈可擊111
[考研] 085600材料與化工調(diào)劑 324分 +8	llllkkkhh 2026-03-18	8/400	2026-03-18 21:01 by Catalysis25
[考研] 328求調(diào)劑，英語六級551，有科研經(jīng)歷 +3	生物工程調(diào)劑 2026-03-16	10/500	2026-03-18 20:41 by Wangjingyue
[考研] 化學工程321分求調(diào)劑 +15	大米飯！ 2026-03-15	18/900	2026-03-18 14:52 by haxia
[考研] 311求調(diào)劑 +6	26研0 2026-03-15	6/300	2026-03-18 14:43 by haxia
[考研] 一志愿西南交大，求調(diào)劑 +4	材化逐夢人 2026-03-18	4/200	2026-03-18 14:22 by 007_lilei
[考研] 304求調(diào)劑 +12	小熊joy 2026-03-14	13/650	2026-03-18 12:34 by Linda Hu
[考研] 278求調(diào)劑 +5	煙火先于春 2026-03-17	5/250	2026-03-18 08:43 by 星空星月
[考研] 085601求調(diào)劑 +4	Du.11 2026-03-16	4/200	2026-03-17 17:08 by ruiyingmiao
[考研] 290求調(diào)劑 +3	p asserby. 2026-03-15	4/200	2026-03-17 16:35 by wangkm
[考博] 26申博 +4	八旬速覽 2026-03-16	4/200	2026-03-17 13:00 by 輕松不少隨
[考研] 302求調(diào)劑 +4	小賈同學123 2026-03-15	8/400	2026-03-17 10:33 by 小賈同學123
[考研] 304求調(diào)劑 +3	曼殊2266 2026-03-14	3/150	2026-03-16 16:39 by houyaoxu
[考研] 321求調(diào)劑 +5	大米飯！ 2026-03-15	5/250	2026-03-16 16:33 by houyaoxu
[考研] 285求調(diào)劑 +6	ytter 2026-03-12	6/300	2026-03-16 15:05 by njzyff
[考研] 材料080500調(diào)劑求收留 +3	一顆meteor 2026-03-13	3/150	2026-03-14 10:54 by peike
[考研] 304求調(diào)劑 +6	Mochaaaa 2026-03-12	7/350	2026-03-13 22:18 by 星空星月
[考研] 材料與化工求調(diào)劑一志愿 985 總分 295 +8	dream…… 2026-03-12	8/400	2026-03-13 22:17 by 星空星月
[考博] 福州大學楊黃浩課題組招收2026年專業(yè)學位博士研究生，2026.03.20截止 +3	Xiangyu_ou 2026-03-12	3/150	2026-03-13 09:36 by duanwu655