[已解决] 采用lsqnonlin进行数据拟合,有一组数据拟合的是折线效果很差,其余几组数据都拟合的很好

[复制链接]
秤砣97 发表于 2021-9-18 17:45:15
本帖最后由 TouAkira 于 2021-9-19 20:54 编辑

大佬们好,我又来了。。。最近在用matlab计算反应动力学参数,采用fmincon函数和lsqnonlin函数结合,也在本论坛大佬的帮助下顺利解出来了结果。但在随后的拟合过程中发现有一两组数据拟合不出来,出来的是折线,而我其余的所有数据拟合的都是平滑的曲线,且拟合的很好(如上传的两张图片),不知道是什么原因。代码如下(其中数据cq0和cexp以及最后的微分方程是拟合不出来的数据组):
  1. function KineticsEst5
  2. clear all; clc
  3. k0=[0.002 0.002];
  4. lb=[0 0]; ub=[1 1];
  5. cq0=[6.3662 0];
  6. tspan=[0;2;5;7;10;20;30;45;60];
  7. cexp=[6.3662
  8. 5.919802056
  9. 5.532036814
  10. 5.501160744
  11. 5.513511172
  12. 5.493012008
  13. 5.53286442
  14. 5.517012582
  15. 5.486454822]; %实验数据

  16. %使用函数fmincon()进行参数估计
  17. [k,fval,flag]=fmincon(@ObjFun4Fmincon,k0,[],[],[],[],lb,ub,[],[],cq0,cexp);
  18. fprintf('\tk=%.4f\n',k)
  19. fprintf('The sum of the squares is:%.1e\n\n',fval)
  20. k_fmincon=k;

  21. %以函数fmincon()估计得到的结果为初值,使用函数lsqnonlin()进行参数估计
  22. k0=k_fmincon;
  23. [k,resnorm,residual,exitflag,output,lambda,jacobian]=lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[],cq0,cexp);
  24. ci=nlparci(k,residual,jacobian);
  25. fprintf('\n\n以fmincon()的结果为初值,使用函数lsqnonlin()估计得到的参数值为:\n')
  26. ci
  27. k
  28. [t c]=ode45(@KineticEqs,tspan(1):1:tspan(end),cq0,[],k);
  29. plot(tspan,cexp,'ko',t,c(:,1),'k-',t,c(:,2),'k-');
  30. disp('t c q ='), disp([t c])
  31. end


  32. %---------------------------------------------------------------------
  33. function f=ObjFun4Fmincon(k,cq0,cexp)
  34. tspan=[0;2;5;7;10;20;30;45;60];
  35. [t c]=ode45(@KineticEqs,tspan,cq0,[],k);
  36. y(:,1)=c(:,1);
  37. f=sum(y(:,1)-cexp.^2);
  38. end

  39. %---------------------------------------------------------------------
  40. function f=ObjFunc4LNL(k,cq0,cexp)
  41. tspan=[0 2 5 7 10 20 30 45 60];
  42. [t c]=ode45(@KineticEqs,tspan,cq0,[],k);
  43. y(:,1)=c(:,1);
  44. f1=y(:,1)-cexp;
  45. f=[f1];
  46. end

  47. %---------------------------------------------------------------------
  48. function dcdt=KineticEqs(t,c,k)
  49. dcdt=[(-k(1)*c(1)+k(1)/0.0426*c(2))
  50.      (2*k(1)*c(1)-(k(2)+2*k(1)/0.0426)*c(2))];

  51. end
复制代码
我尝试了很多次,改初值以及上下限,实在没办法只好又来求助各位大佬,帮我看这是什么原因引起的,万分感谢!

出问题的拟合

出问题的拟合

正常的拟合

正常的拟合

最佳答案


TouAkira 发表于 2021-9-19 08:42:19
应该就是单纯的微分方程组公式不适用。

你目前这个
dcdt=[(-k(1)*c(1)+k(1)/0.0426*c(2))
     (2*k(1)*c(1)-(k(2)+2*k(1)/0.0426)*c(2))];
实际上换一下元,可以写成
dcdt=[ A * c(1)+ B * c(2)
     C * c(1)+ D * c(2) ];
的形式,而且c的初值都给定了,这就是个简单的线性微分方程组,是能求出解析解的
syms y1( t ) y2( t ) A B C D Y_10
Eqns = [ diff( y1, t ) == A * y1 + B * y2;
    diff( y2, t ) == C * y1 + D * y2 ];
Conds = [ y1( 0 ) == Y_10; y2( 0 ) == 0 ];
[ y1Sol( t ), y2Sol( t ) ] = dsolve( Eqns, Conds )
这样求出来的解析解(老长一大串)之后,再把
A = ( -k( 1 ) )
B = ( k( 1 ) / 0.0426 )
C = ( 2 * k( 1 ) )
D = ( - k( 2 ) - 2 * k( 1 ) / 0.0426 )

都代入y1Sol中去,等于你的第一个反应物的浓度随时间变化的函数表达式就求出来了,等于是对y1Sol(它是待拟合参数k(1) k(2) 时间t 以及两个初始浓度的函数,但实际上由于其他都是已知的,仅含有待拟合参数k(1) k(2),这个非线性拟合的计算量,比原先每次计算误差都要求解一遍微分方程组,要小得多得多)

然而由于你这ABCD四个系数彼此“不独立”,或者说ABC三者直接就线性相关,这有可能导致y1Sol只能为某一系列特定曲线,不幸的是,这一系列特定曲线全都跟你的实验数据趋势不怎么匹配

这种也很常见,由于各种因素导致预测公式不适用,比如最简单的,反应级数真的是1级(这个直接影响微分方程组右侧c(1) c(2)的幂次)么,等等。
回复此楼

6 条回复


TouAkira 发表于 2021-9-19 08:42:19
应该就是单纯的微分方程组公式不适用。

你目前这个
dcdt=[(-k(1)*c(1)+k(1)/0.0426*c(2))
     (2*k(1)*c(1)-(k(2)+2*k(1)/0.0426)*c(2))];
实际上换一下元,可以写成
dcdt=[ A * c(1)+ B * c(2)
     C * c(1)+ D * c(2) ];
的形式,而且c的初值都给定了,这就是个简单的线性微分方程组,是能求出解析解的
syms y1( t ) y2( t ) A B C D Y_10
Eqns = [ diff( y1, t ) == A * y1 + B * y2;
    diff( y2, t ) == C * y1 + D * y2 ];
Conds = [ y1( 0 ) == Y_10; y2( 0 ) == 0 ];
[ y1Sol( t ), y2Sol( t ) ] = dsolve( Eqns, Conds )
这样求出来的解析解(老长一大串)之后,再把
A = ( -k( 1 ) )
B = ( k( 1 ) / 0.0426 )
C = ( 2 * k( 1 ) )
D = ( - k( 2 ) - 2 * k( 1 ) / 0.0426 )

都代入y1Sol中去,等于你的第一个反应物的浓度随时间变化的函数表达式就求出来了,等于是对y1Sol(它是待拟合参数k(1) k(2) 时间t 以及两个初始浓度的函数,但实际上由于其他都是已知的,仅含有待拟合参数k(1) k(2),这个非线性拟合的计算量,比原先每次计算误差都要求解一遍微分方程组,要小得多得多)

然而由于你这ABCD四个系数彼此“不独立”,或者说ABC三者直接就线性相关,这有可能导致y1Sol只能为某一系列特定曲线,不幸的是,这一系列特定曲线全都跟你的实验数据趋势不怎么匹配

这种也很常见,由于各种因素导致预测公式不适用,比如最简单的,反应级数真的是1级(这个直接影响微分方程组右侧c(1) c(2)的幂次)么,等等。
回复此楼

秤砣97 发表于 2021-9-19 14:13:08
TouAkira 发表于 2021-9-19 08:42
应该就是单纯的微分方程组公式不适用。

你目前这个

感谢大佬第二次帮助! 我按照您给的这个代码得到了解析解,然后您的提示把原程序做了修改,不知道我理解的对不对啊,不过确实也能解出来,跟之前解微分方程的方法得到的结果是一样的,但是我之前那个拟合不出来的数据用这个求解析解的方法还是拟合不出来平滑的曲线,所以确实可能是大佬说的这个模型微分方程组不适用于我这一组数据,不过这也没关系,这一组数据不适用也能够解释。再次感谢大佬的帮助指点!

TouAkira 发表于 2021-9-20 08:53:08
论坛的"基本属性"中含有四项,分别是"回帖仅作者可见"、"回帖倒序排列"、"接收回复通知"、"使用个人签名",请不要滥用"回帖仅作者可见"这个功能!

回答者可能自己辛苦编写了代码,不想透露智力成果的细节,这还能够理解;提问者设置"回帖仅作者可见",今后如果再有其他坛友遇到相似问题、点击进来想要查看如何解决却发现什么有效信息都看不到,那论坛还有什么意义

一方面想要"给代码治病",另一方面又不想让其他人知道病是怎么治疗的,太自私了!

如果有需要保密的数据,那就自己生成能复现问题的模拟数据再发帖!

如果有敏感不能公开的信息,那就私下花钱去咨询专家!

秤砣97 发表于 2021-11-3 13:39:08
TouAkira 发表于 2021-9-20 08:53
论坛的"基本属性"中含有四项,分别是"回帖仅作者可见"、"回帖倒序排列"、"接收回复通知"、"使用个人签名", ...

不好意思才看到消息,好像已经改回去了,实在抱歉,下次一定注意!

秤砣97 发表于 2021-11-3 13:48:20
TouAkira 发表于 2021-9-20 08:53
论坛的"基本属性"中含有四项,分别是"回帖仅作者可见"、"回帖倒序排列"、"接收回复通知"、"使用个人签名", ...

当时是考虑怕被导师责怪泄露代码就设置了这个,最近文章已经投出去了就上线来解除这个限制,发现已经被管理员解除了,还被警告了。。感谢管理员提醒,今后一定注意,也希望管理员理解。

秤砣97 发表于 2021-11-3 15:07:53
TouAkira 发表于 2021-9-20 08:53
论坛的"基本属性"中含有四项,分别是"回帖仅作者可见"、"回帖倒序排列"、"接收回复通知"、"使用个人签名", ...

之前是因为怕导师责怪泄露代码所以设置了回帖仅作者可见,最近文章准备可以投了就上线把这个限制给取消掉,发现管理员已经取消了,也被警告了。。感谢管理员提醒,下次一定注意,也希望管理员理解!
您需要登录后才可以回帖 登录 | 注册

本版积分规则

相关帖子
相关文章
热门教程
站长推荐
快速回复 返回顶部 返回列表