zhangqq72270

實習(xí)版主

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

» 猜你喜歡

292分，材料與化工，申請調(diào)劑已經(jīng)有13人回復(fù)
316求調(diào)劑已經(jīng)有7人回復(fù)
材料專碩322分已經(jīng)有4人回復(fù)
286分調(diào)劑已經(jīng)有16人回復(fù)
337求調(diào)劑已經(jīng)有8人回復(fù)
一志愿北交大材料工程總分358 已經(jīng)有7人回復(fù)
362求調(diào)劑已經(jīng)有13人回復(fù)
332求調(diào)劑已經(jīng)有7人回復(fù)
342求調(diào)劑已經(jīng)有11人回復(fù)
319求調(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

引用回帖:

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

查看全部 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的回帖

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

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

[考研] 320分，材料與化工專業(yè)，求調(diào)劑 +11	一定上岸aaa 2026-03-27	15/750	2026-04-01 16:57 by oooqiao
[考研] 一志愿211，335分，0856，求調(diào)劑院校和導(dǎo)師 +11	傾____蕭 2026-03-27	12/600	2026-04-01 15:27 by oooqiao
[考研] 求生物學(xué)調(diào)劑 +7	15172915737 2026-04-01	7/350	2026-04-01 15:27 by 閑人終南山
[考研] 材料專碩調(diào)劑 +14	椰椰。 2026-03-29	14/700	2026-04-01 15:26 by oooqiao
[考研] 330分求調(diào)劑 +11	qzenlc 2026-03-29	11/550	2026-04-01 14:32 by chenqifeng666
[考研] 材料專業(yè)求調(diào)劑 +7	月月鳥木 2026-04-01	7/350	2026-04-01 11:36 by wangjy2002
[考研] 314求調(diào)劑 +6	1xiaojun23 2026-03-31	6/300	2026-04-01 11:28 by 今天星期八LQ
[考研] 化學(xué)0703 調(diào)劑 306分一志愿211 +12	26要上岸 2026-03-28	12/600	2026-04-01 11:10 by chemdavid
[碩博家園] 博一被送出聯(lián)培感覺不適應(yīng)怎么辦 +3	全村的狗 2026-03-31	3/150	2026-04-01 10:44 by 328838485
[考研] 358求調(diào)劑 +3	王向陽花 2026-03-31	3/150	2026-04-01 09:56 by zzchen2000
[考研] 求調(diào)劑，一志愿南京航空航天大學(xué) ，080500材料科學(xué)與工程學(xué)碩，總分289分 +10	@taotao 2026-03-29	10/500	2026-04-01 09:30 by oooqiao
[考研] 085602 307分求調(diào)劑 +10	不知道叫什么！ 2026-03-26	10/500	2026-03-31 19:53 by Dyhoer
[考研] 346求調(diào)劑一志愿070303有機化學(xué) +11	蘿卜燉青菜 2026-03-28	12/600	2026-03-31 17:54 by xhai2011
[考研] 調(diào)劑求院校招收 +7	鶴鯨鴿 2026-03-28	7/350	2026-03-31 11:21 by oooqiao
[考研] 279求調(diào)劑 +4	蝶舞輕繞 2026-03-29	4/200	2026-03-29 09:45 by laoshidan
[考研] 339求調(diào)劑，想調(diào)回江蘇 +6	烤麥芽 2026-03-27	8/400	2026-03-28 10:40 by 烤麥芽
[考研] 081200-314 +3	LILIQQ 2026-03-27	4/200	2026-03-28 09:41 by 保護地球你我做?/a>
[考研] 308求調(diào)劑 +7	墨墨漠 2026-03-27	7/350	2026-03-28 07:43 by 熱情沙漠
[考研] 321求調(diào)劑 +6	wasdssaa 2026-03-26	6/300	2026-03-26 20:57 by sanrepian
[考研] 中國科學(xué)院深圳先進技術(shù)研究院-光纖傳感課題組招生-中國科學(xué)院大學(xué)、深圳理工大學(xué)聯(lián)培 +5	YangTyu1 2026-03-26	5/250	2026-03-26 18:27 by 貓咪貓咪呀