| 查看: 1680 | 回复: 6 | |||
[交流]
使用lsqnonlin函数优化动力学参数，总是得不到合理的结果
|
|
大家好: 我最近在使用matlab回归动力学参数，首先是使用4阶龙哥库塔算法进行积分，之后使用lsqnonlin进行非线性最小二乘拟合，运行已经写好的m文件，总是得不到合理的结果，不知道是不是自己程序编写有问题，请大神指导，谢谢。 下面粘贴上自己的m文件： %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% function zqq100 % 动力学ODE方程模型的参数估计 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% clear all clc kz0 = [1.73 -3 4.3 3.4 29.4 26.46 22.8 29.4]; % 参数初值 lb = []; % 参数下限 ub = []; % 参数上限 x0 = [72516.987800 7922.540990 2213.519110 2440.339360]; KineticsData; yexp = ExpData(:,2:5); % yexp: 实验数据[x1 x2 x3 x4] %再参数化 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 % 优化参数设置 options=optimset ('Display','iter','TolX',1e-25,'TolFun',1e-25,'LargeScale','off','MaxIter',2,'MaxFunEvals',inf); % 使用函数lsqnonlin()进行参数估计 [kz,resnorm,residual,exitflag,output,lambda,jacobian] = ... lsqnonlin(@ObjFunc4LNL,kz0,lb,ub,options,x0,yexp,Taverage); ci = nlparci(kz,residual,jacobian); %把回归参数转化为自己需要的参数 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; % ------------------------------------------------------------------ % 输出结果 fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\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); 数据文件 % 动力学数据: % 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 ] |
回复此楼» 猜你喜欢
化学调剂0703
已经有8人回复
-
梁成伟老师课题组欢迎你的加入
已经有8人回复
-
环境工程调剂
已经有6人回复
-
326求调剂
已经有7人回复
-
一志愿985，本科211，0817化学工程与技术319求调剂
已经有5人回复
-
今年的国基金是打分制吗？
已经有3人回复
-
304求调剂
已经有5人回复
-
一志愿苏州大学材料工程（085601）专硕有科研经历三项国奖两个实用型专利一项省级立项
已经有5人回复
-
304求调剂
已经有4人回复
-
本人考085602 化学工程 专硕
已经有14人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- matlab的lsqnonlin函数怎么用
已经有7人回复
- 紧急求助，利用Matlab对实验数据进行拟合求解参数。
已经有27人回复
- 关于matlab的参数估计
已经有15人回复
- 怎样用一组参数同时拟合两个曲线--matlab
已经有5人回复
- Matlab最小二乘参数优化
已经有3人回复
- 帮忙给看看这个matlab优化函数 问题
已经有8人回复
- 动力学参数拟合
已经有26人回复
- MATLAB非线性优化拟合怎么改才正确
已经有3人回复
- matlab fmincon优化函数 怎样知道每次得到的迭代点
已经有3人回复
- 求助 用Oringin拟合转状态方程
已经有6人回复
- 请教matlab反应动力学参数估计遇到的问题，谢谢
已经有15人回复
- 拟合非线性方程系数的精度问题
已经有5人回复
- matlab拟合方程参数时初值的选择
已经有15人回复
- matlab拟合模型参数
已经有5人回复
- lsqnonlin函数拟合微分方程组参数拟合问题
已经有10人回复
- matlab非线性参数拟合问题
已经有7人回复
- 【求助】催化反应动力学matlab计算各基元反应的速率常数时，该如何避免较小量被忽略？
已经有3人回复
- 【求助】催化反应动力学
已经有5人回复
- 【求助】matlab 二次规划的优化的问题
已经有4人回复
- 【求助】使用Matlab拟合反应动力学方程问题
已经有7人回复
- 【求助】使用Matlab预估动力学方程问题
已经有13人回复
» 抢金币啦！回帖就可以得到:
-
山东征女友，坐标济南
+1/169
-
博士后出站经验请教
+1/78
-
哈工大本部-电化学新能源方向博士捡漏名额
+1/77
-
加拿大阿尔伯塔大学招收电磁、无线通信、机器学习方向全奖硕士/博士/博士后/访问学者
+1/56
-
何时使用 CODA™ 科里奥利质量流量仪表- 艾里卡特（Alicat）
+2/52
-
华南师范大学（211）博士招生- 电子、自动化、机械、生物学、物理相关专业
+2/48
-
澳大利亚西澳大学（UWA）张金强课题组招收全奖博士生（化工新能源方向）
+1/32
-
2026年工科硕士调剂-上海大学全国重点实验室团队-材料数据挖掘方向-研究生3-5人
+1/31
-
诚聘助理研究员 - 物理/材料/机械相关专业博士
+2/20
-
医学327分专硕调剂
+1/16
-
北京某研究院生物化学与分子生物学相关专业申请考核制博士招生（春季）
+1/15
-
香港大学David Srolovitz教授课题组招聘「香江学者」博士后
+1/13
-
加拿大阿尔伯塔大学招收电磁、无线通信、机器学习方向全奖硕士/博士/博士后/访问学者
+2/10
-
上海科技大学招收有机化学方向招收2026级工程博士
+1/8
-
上海师范大学化学与材料科学学院任新意副研究员招收调剂学生3-4名（有机化学专业）
+1/5
-
复旦大学集成电路学院程增光课题组急聘科研助理
+1/3
-
北京高校副校长团队招收机械类，环境类学硕和专硕
+1/3
-
085404 270求调剂，b区求导师收留接受跨专业
+1/2
-
温州大学招收2026年入学博士研究生（化学、材料、环境）
+1/1
-
339求调剂
+1/1
2楼2013-08-14 10:22:39
3楼2013-08-14 10:25:21
★
小木虫: 金币+0.5, 给个红包，谢谢回帖
小木虫: 金币+0.5, 给个红包，谢谢回帖
| 强非线性的是不是需要用其他算法，比如GA这类全局搜索的，这种问题如果给的初值不在平衡点的吸引域内是不是就没法获得很好的效果。 |
» 本帖已获得的红花（最新10朵）
zhangqq722704楼2013-08-16 11:38:32
5楼2013-08-16 11:48:16
6楼2013-08-16 14:49:19
7楼2013-08-16 14:49:43
