| 查看: 1709 | 回复: 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 ] |
回复此楼» 猜你喜欢
085600材料与化工301分求调剂院校
已经有5人回复
-
081700，311，求调剂
已经有15人回复
-
一志愿北京化工085600 310分求调剂
已经有18人回复
-
材料与化工371求调剂
已经有14人回复
-
336材料与化工085600求调剂
已经有7人回复
-
材料334求调剂
已经有18人回复
-
331求调剂
已经有8人回复
-
332求调剂
已经有17人回复
-
一志愿南京航空航天大学 材料与化工329分求调剂
已经有4人回复
-
材料专硕322
已经有7人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- 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/760
-
武汉纺织大学全国重点实验室招收材料、化学、纺织相关研究生
+1/64
-
一志愿前4位代码0703均可报名-河北大学分析化学专业招收调剂考生
+1/47
-
安徽工业大学工程研究院接受调剂
+1/39
-
欢迎材料、化工、林业工程等相关方向学生报考/调剂 福建农林大学
+1/37
-
五邑大学马爱军课题组招收药学、生物与医药调剂生
+1/37
-
抗病毒性传染病创新药物全国重点实验室—接收化学、生物、药学、生物与医药调剂考生
+1/22
-
三峡大学国家青年人才课题组招2026年化学/材料/化工类硕士调剂生
+1/20
-
青科大接收材料类、化工化学类调剂研究生（须考数学）
+1/18
-
齐鲁工业大学-材料方向招收两名专硕调剂
+1/16
-
广东工业大学生态环境与海洋学院张轩教授课题组招收分离膜领域2026年入学博士生1人
+1/9
-
南通大学-交通与土木工程学院 招调剂学生
+1/8
-
江苏大学化学化工学院2026年硕士研究生招生调剂公告
+1/5
-
招收2026级研究生
+1/5
-
哈工程烟台院材化专业2026调剂信息
+1/4
-
武汉纺织大学全国重点实验室吕永钢教授课题组招收材料、化学、生物相关专业硕士研究生
+1/4
-
南昌大学工程建设学院2026年招聘2名科研助理
+1/4
-
调剂招生化学化工（0817、0856）
+1/3
-
哈尔滨医科大学李老师课题组招收生物信息学方向博士、硕士研究生
+1/3
-
济南大学2026年硕士研究生招生预调剂公告
+1/2
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
