| 查看: 1968 | 回复: 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 ] |
回复此楼» 猜你喜欢
售SCI文章，我:8O.5.5.1O.54,科目全,可十急
已经有4人回复
-
售SCI一区T0P文章，我:8O.55.1.O.54,科目全,可伽急
已经有4人回复
-
售SCI一区T0P文章，我:8.O.55.1.O54,科目全,可伽急
已经有4人回复
-
售SCI一区文章，我:8O5.5.1.O5.4,科目全,可伽急
已经有6人回复
-
售SCI一区文章，我:8O5.5.1.O5.4,科目全,可伽急
已经有6人回复
-
售一区SCI文章T0P，我:8O.551.O54,科目全,可十急
已经有4人回复
-
某top大学教授说“能够在市场中兑现的能力才是真能力”无比同意！
已经有13人回复
-
同一个分子，一条来自文献，一条来自AI——不告诉你答案，你会选哪条？
已经有7人回复
-
售SCI文章，我:8O.5.5.1O.54,科目全,可十急
已经有4人回复
-
售SCI一区T0P文章，我:8.O.55.1.O54,科目全,可伽急
已经有4人回复
» 本主题相关价值贴推荐，对您同样有帮助:
- 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人回复
» 抢金币啦！回帖就可以得到:
-
持续创业者，A8离异男征婚，事业遇到瓶颈，寻找生活伴侣
+1/479
-
专利代理人
+1/96
-
复旦大学生物与信息方向博士研究生等招聘
+1/96
-
急招哈尔滨工业大学（深圳）机械学科2027级"申请-考核"博士研究生1名
+1/88
-
南科大薛亚辉课题组诚聘离子输运、二维材料、原子力显微镜等方向“快响行动”博士生
+1/66
-
南方科技大学力航系急招收纳米材料、力学、机械等方向推免硕士研究生两名
+1/62
-
华南师范大学胡勇军教授课题组招收2027年推荐面试攻读研究生（推免生）
+1/31
-
澳大利亚昆士兰大学博士生招聘/气体发酵、酵母代谢工程与合成生物学方向
+2/28
-
哈工大深圳（机器人方向）-姜欣课题组招收2027年博士研究生（10月7日截止）
+1/16
-
中山大学(深圳校区)-新病毒及宿主交互课题组招收博士生
+1/13
-
【西北工业大学-985西安】材料学院-余泓课题组-诚招2027年博士
+1/12
-
中科院东莞材料所｜新能源材料研究部招聘研究员 / 副研究员（可中央事业编制）
+1/2
-
2027 International Conference on Healthcare, Medicine and Medical Science
+1/1
-
ICTCPS（纺织、着色与高分子科学）Textiles, Coloration and Polymer Science
+1/1
-
中科院东莞材料所二维材料与智能器件方向招聘
+1/1
-
从NSA网络攻击事件看漏洞智能分析——AI大模型本地部署与网络安全算力工作站配置方案
+1/1
-
武汉大学电信学院-梁国镇课题组-诚招硕士博士 - 深紫外光源方向 - 2027年入学
+1/1
-
2027 CoASS（天文学与空间科学）
+1/1
-
2027 IFPAM（基础物理与先进材料）
+1/1
-
浙江大学Kemal Celebi教授课题组招收2027年秋季博士生-材料/化学/化工
+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
