| 查看: 1885 | 回复: 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 ] |
回复此楼» 猜你喜欢
科研孤儿太难了
已经有9人回复
-
时间戳又变了
已经有15人回复
-
时间戳变了，能看出什么问题？
已经有23人回复
-
只有每年这种时候来逛逛小木虫
已经有22人回复
-
今天系统多次维护，明天很可能放榜！
已经有10人回复
-
放榜前的不淡定
已经有12人回复
-
言之凿凿，明天官网会通知放榜？
已经有16人回复
-
重要消息，中午系统在维护
已经有12人回复
-
哪位老哥知道今年的国自然具体哪一天放榜？
已经有20人回复
-
今天基金会出结果吗？20260819
已经有15人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- 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人回复
» 抢金币啦！回帖就可以得到:
-
散金祈福
+5/4480
-
告别流水线爽文！108万字完结仙侠《女娲石之灵石奇缘（仙界篇）》值得细品
+1/4065
-
每年这时候我都再仔细读一遍《范进中举》
+1/484
-
华西生物制药研究院招聘科研助理
+2/418
-
五大联赛看好哪支球队
+5/250
-
散金求好运
+1/203
-
上海交通大学张航课题组招聘博士后（电化学储能）
+1/93
-
山东征女友，坐标济南
+1/91
-
Polymers 期刊特刊 组稿
+1/90
-
哈尔滨工业大学（校本部）电子与信息工程学院2027级博士招生
+1/68
-
面上基金祈福
+1/52
-
祈祷老公中国基
+1/21
-
pcold 原核表达载体
+1/13
-
2027环境工程博士申请
+1/9
-
新审查指南施行半年，AI方向的专利申请是不是更看"技术味"了？
+1/5
-
南京医科大学-冯璋课题组-招收2027级生化与分子生物学专业博士和硕士研究生
+1/3
-
北理工集成电路杰青团队 | 诚招科助理
+1/2
-
清华大学环境课题组招聘环境工作专业客座生研究生2-3名（应用导向）
+1/2
-
河海大学海洋环境污染与生态毒理学研究组招收2027级硕士生 及博士研究生
+1/2
-
清华大学环境课题组招聘环境工作专业客座生研究生2-3名（应用导向）
+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
