| 查看: 1967 | 回复: 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 ] |
回复此楼» 猜你喜欢
某top大学教授说“能够在市场中兑现的能力才是真能力”无比同意！
已经有13人回复
-
同一个分子，一条来自文献，一条来自AI——不告诉你答案，你会选哪条？
已经有7人回复
-
售SCI文章，我:8O.5.5.1O.54,科目全,可十急
已经有4人回复
-
售SCI一区T0P文章，我:8.O.55.1.O54,科目全,可伽急
已经有4人回复
-
售SCI一区文章，我:8.O.55.1.O.54,科目齐全,可伽急
已经有4人回复
-
售SCI一区T0P文章，我:8.O.55.1.O.54,科目齐全,可+急
已经有4人回复
-
售SCI-T0P文章，我:8O.5.5.1.O.54,科目齐全,可+急
已经有3人回复
-
售SCI文章，我:8O.5.5.1O.54,科目全,可十急
已经有3人回复
-
售一区SCI文章T0P，我:8O.551.O54,科目全,可十急
已经有3人回复
-
售SCI文章，我:8O5.5.1.O.54,科目齐全,可+急
已经有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人回复
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
