找回密码
 注册
关于网站域名变更的通知
查看: 653|回复: 3
打印 上一主题 下一主题

常微分方程组参数拟合问题

[复制链接]

该用户从未签到

跳转到指定楼层
1#
发表于 2021-3-2 10:28 | 只看该作者 回帖奖励 |倒序浏览 |阅读模式

EDA365欢迎您登录!

您需要 登录 才可以下载或查看,没有帐号?注册

x
代码如下;' S: I7 V6 c; g4 l1 P6 ~8 g! e
function Kinetics4
% x8 d8 a8 `+ `' D7 ?# wclear all
- l8 U( E% m# t3 Pclc9 _$ z) Q& B$ m. n( U. X1 @  o
k0 = [2.159*10^8 4.734*10^6 174686 149892 0.99 0.7 0.5 0.78];%参数初值,其实是不知道的  " d# J$ j8 |+ H1 W0 U6 u
lb = [-inf -inf -inf -inf -inf -inf -inf -inf];                   % 参数下限2 T1 t% h( a2 s& W
ub = [inf inf inf inf inf inf inf inf];    % 参数上限6 j7 i( s% x7 ~* y; j. f  z
u0 = [0,0]; %y初值* T" Q8 u$ ^1 U8 p& y
a=[0;0.06680583;0.192617534;0.301693824;0.419783943;0.527041818;0.598345407;0.65055108;0.6876854;0;0.0295189820000000;0.0760025420000000;0.143652564000000;0.303435528000000;0.419715167000000;0.525496977000000;0.565069320000000;0.743536298000000;0;0.0411725930000000;0.0608472370000000;0.0986899250000000;0.241132264000000;0.324760943000000;0.466871464000000;0.531688237000000;0.683614442000000];. I- w* o* X" o1 @1 I
b=[0;0.2378;0.3288;0.3556;0.39;0.3877;0.3766;0.3580;0.3384;0;0.1373;0.2123;0.2755;0.3561;0.3549;0.3475;0.3486;0.2518;0;0.115879961000000;0.178182567000000;0.219982937000000;0.315058656000000;0.339244190000000;0.348664305000000;0.351963292000000;0.313743618000000];8 W2 s0 _' y: n+ _' \
yexp=[a,b];& }! j/ A! L7 a  s
% 使用函数lsqnonlin()进行参数估计  h  B) n" J) X4 B$ |" \* R. x

5 O" z* a* P$ n& _, K& I# |4 j[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
1 ]! i6 [0 w* x8 E0 u    lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[],u0,yexp);      
% R/ [, j$ k2 p) n+ q, nci = nlparci(k,residual,jacobian);& g% b0 G6 K, i5 H& }
fprintf('\tk1 = %.4f ± %.4f\n',k(1),ci(1,2)-k(1))
* J' `8 j5 O# c* y4 _fprintf('\tk2 = %.4f ± %.4f\n',k(2),ci(2,2)-k(2))
: c6 ^" S/ n9 w. [" X4 x$ [fprintf('\tk3 = %.4f ± %.4f\n',k(3),ci(3,2)-k(3))  / T2 [$ l! D) F4 @* w2 v( y$ Q9 F+ y
fprintf('\tk4 = %.4f ± %.4f\n',k(4),ci(4,2)-k(4))1 T  q1 f) w9 r: V* a: M( b
fprintf('\tk5 = %.4f ± %.4f\n',k(5),ci(5,2)-k(5))
! g4 [( c  M3 f, {; Ffprintf('\tk6 = %.4f ± %.4f\n',k(6),ci(6,2)-k(6))$ W" _1 S) f2 J1 R# P. b9 y/ @
fprintf('\tk7 = %.4f ± %.4f\n',k(7),ci(7,2)-k(7)); @2 O) \# K0 L3 ^+ g
fprintf('\tk8 = %.4f ± %.4f\n',k(8),ci(8,2)-k(8))
+ ~/ p$ e( A! }. x; Gfprintf('  The sum of the squares is: %.1e\n\n',resnorm)
  i1 g( Y1 F, T" pfprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')6 j: B% Q2 S' d6 l
%目标函数
5 R; c4 w* j: l( h* t7 w. _# yfunction f =ObjFunc4LNL(k,u0,yexp)4 i7 `0 f9 W. W9 w& L8 ^
global theta Pt Ac R T Fa0- f/ N& ]0 N( L/ U
theta=zeros(5,1);
0 ]# H& a% u5 M( c9 P+ |+ b4 \/ Zr0= 1.24*10^-3; % Outer diameter  of the membrane, [m]
% j# w3 }+ J+ P$ ~/ c' ~" i3 Gri= 9.4*10^-4;% Inner diameter  of the membrane, [m]
; A! Y% ]1 V7 q' b7 D7 z5 Or=(r0-ri)/log(r0/ri);% equivalent diameter ,[m]
7 `5 A6 n1 U- X7 p; xPi=3.14159;9 v9 Q$ H9 s6 N/ ?' n/ |& G. U
L=0.05;% membrane length,[m]
+ o9 x% ?0 g7 T$ A: oAc=Pi*r*L;% membrane area [m^2]- o4 E0 [) }4 C  G4 G. b) c3 E+ c
R=8.314; % [J/K.mol]% g5 F& J' s5 J$ U) S, L+ t
Pt=101.325;
. F' T9 U- N1 j2 g7 t+ AFtot=[2.678*10^-2;2.678*10^-2*1.5;2.678*10^-2*2];, Y# B: \9 G- a% |$ i  V
Kmax=length(Ftot);
" ]- r7 v2 H9 M- c* D2 D9 V/ g; n
4 s' c4 Q- t4 S5 n+ o, \- vT0=600+273;% Inlet temperature [K]  x9 U. h. x* X$ [: r& c7 h1 O
nmax=9;
# q9 [3 Q" A  s: ~; D& c6 Z& uX1=zeros(nmax,Kmax);X2=zeros(nmax,Kmax);
  P3 x( ^" ^% K" F0 Jfor K=1:Kmax9 t/ _% \% c. n' k- s
Fa0=Ftot(K); % If Pt is the parameter/Matrix- l# X3 k) ], d3 l$ m1 L6 {
%Pt=Pt0+(k-1)*1; % If Pt is an axis
( z$ [: S, }) t0 y. k%L=Lsize(k); % If L is the parameter/Matrix1 T; K/ ^9 F3 J" A$ a
%theta(2)=H2Oratio+(k-1)*1; % If m is the parameter/Matrix
! r  r6 _8 m9 f$ i2 g4 stheta(2)=3; % F0(H2O)/F0(CH4)
! d7 A9 W9 ?$ p+ H2 Utheta(3)=0; % F0(CO)/F0(CH4)
) `# x( w: f  j; c. ntheta(4)=0; % F0(CO2)/F0(CH4)9 N2 l$ e1 {5 O$ b# }2 _
theta(5)=0; % F0(H2)/F0(CH4)
: J8 W' k: k9 W, r0 i- i7 U5 ^. I3 W/ ffor n=1:nmax
5 c" D$ ~- \% o  r, g, ?0 ~%delta=delta0+(n-1)*1e-6; % If delta is the parameter/Matrix
2 \) O( l, k3 |%rit=r0+delta; % inner tube radius [m]( ?; K2 D" v  M
%sweep=s0+(n-1); % If sweep is the parameter/Matrix. O+ E: I8 I" B. Z' I, F7 b" M
%T=T0; % If T is a constant
4 n- w2 |# S% t* @$ CT=T0+(n-1)*50;; O  U8 T3 n( ]9 B4 p$ W
ksispan=(0:0.1:1); % precision
/ h: f6 Y4 r" f( l4 `2 ^% N$ N[ksi,S1]=ode23s(@KineticEqs,ksispan,u0,[],k); % PBR reactor
2 p% |+ P4 m; M$ C0 D- Q1 X  CX1(n,Kmax)=real(S1(end,1)); % methane finale conversion
  `8 Y. h6 q$ J+ G/ S% I& mX2(n,Kmax)=real(S1(end,2)); % carbon dioxide finale conversion  ]6 s! w: ?* }' m9 z
end$ f* H; X# ]" [" y& r6 R
end7 O, y' U5 E; T6 q* h0 s
f1=X1(n,Kmax)-yexp(:,1);
. v) A; J: C9 k6 Uf2=X2(n,Kmax)-yexp(:,2);: u+ {- a% z3 K( D3 S( W! A7 R2 M4 r
f=[f1;f2];$ `  |( b. p9 j! D2 C% `3 ]
%T0=T0-273.15;T=T-273.15; % Temperature expressed in Celsius  
3 Y* {+ P/ B( G# w%Tspan=(T0:1:T);5 i, s4 V6 a3 R7 l' v2 J6 a9 a, s3 g/ ?+ ~
%Pspan=(Pt0:1t);
2 a7 Q- T/ ?2 D%DSPan=(delta0*1e6:1:delta*1e6);  K0 t8 M4 c( n1 b" [. S8 D7 y
%Sspan=(s0:1:sweep);2 ~* m3 G0 p  \' o
%D(:,=X3(:,-X1(:,; % Delta = XCH4(MR) - XCH4(PBR)
. \4 z8 p% P, g  N
1 m6 x  X7 u2 X. T# }  ufunction diff=KineticEqs(ksi,u,k); [' A0 n5 v1 a
global theta  Pt  Ac Fa0 T R
7 |5 @2 ]$ m* a2 }) A9 IS=Pt/(1+theta(2)+theta(3)+theta(4)+theta(5)+2*(u(1)+u(2)));& {% _) O. U5 T/ {- ~3 A- g2 {& g
P(1)=(1-u(1)-u(2))*S;7 o8 }9 i/ M$ N
P(2)=(theta(2)-u(1)-2*u(2))*S;) U! z, ?# Q! ?: I9 j3 c
P(3)=(theta(3)+u(1))*S;
* e  u- @7 L( J+ @" t) aP(4)=(theta(4)+u(2))*S;
" f- |) P. b9 t" QP(5)=(theta(5)+3*(u(1)+u(2))+u(2))*S;2 B* v2 L5 B* e9 b" s$ v
% Differential System* D* [9 U9 ^6 x4 }1 D  W( ?! b
diff1=Ac/Fa0*k(1)*exp(-k(3)/R/T)*P(1)^k(5)*P(2)^k(7);: v  n( F6 N7 |! `6 l( H& M
diff2=Ac/Fa0*k(2)*exp(-k(4)/R/T)*P(1)^k(6)*P(2)^k(8);
' y( q$ {# T8 rdiff=[diff1 diff2]';8 y/ X( J+ J1 C
4 n# v% w  r0 }0 q" d6 _
运行结果:
) s  v  J9 o' TIn nlparci (line 104)
& e& s' G* ]- W( n: ?( x. D  In Kinetics4 (line 15)
- h- W& {! o; r3 M/ X3 o% z3 T: ?警告: 矩阵为奇异工作精度。
2 k- q' a0 U5 q# \$ ?' c' e& R        k1 = 215900000.0000 ± NaN8 \9 P9 R6 i; _6 e! V) Q
        k2 = 4734000.0000 ± NaN
' ?6 b2 \& ?/ k/ W  M: s        k3 = 174686.0000 ± NaN
& H+ g: Y7 Y4 `% J1 u7 i  i# {        k4 = 149892.0000 ± NaN  Y: D8 P3 f" h& n% g# O: e& h
        k5 = 1.2381 ± NaN- U0 m/ Q, v+ v' J6 U; d' V
        k6 = 0.5426 ± NaN; }8 d0 v' ~# L. O; W
        k7 = -0.3426 ± NaN
, ?7 I5 k- N% z- u# k& X        k8 = 0.4455 ± NaN4 X1 v" O# d; [4 g3 _
  The sum of the squares is: 2.0e+00
5 U* ]/ `" R6 a% w( F+ v, ~

该用户从未签到

2#
发表于 2021-3-2 10:48 | 只看该作者
S=Pt/(1+theta(2)+theta(3)+theta(4)+theta(5)+2*(u(1)+u(2)));% x+ w; i& s5 j, d
你确定没写错??你自己文档里面写的分母可是( 1 + m + 2 * x_1 + x_2 )( G/ H) N" _' |/ `$ n

0 ?/ K) r6 E- ], Q) d: {然后,这个拟合太耗时,你给的参数上下界太大了,导致随机生成的数据,很多组代进微分方程组去运算起来非常非常慢,至少先去看文献也好查经验公式也好,把数量级、正负号之类的给限制一下,要不然很难做优化。: f  c& u1 k, ]& O/ e
& @: e) d( N3 Z- }, z
你给的原始参数的估计值的效果
- |1 k2 D# X# S' v3 D, }, }
1 c. i5 R6 i: g; D/ O: k
; ^8 m# d' ?6 f, X1 W, s
9 K' k8 u. ]" ?' G拟合得到的两组参数的效果18 Y8 j5 y& `2 g, k5 B. I( v: e3 ^3 h

4 {4 u- T1 n0 V1 ?
* M  Z1 E4 ^" [5 a! V拟合得到的两组参数的效果2+ J! K2 @4 M( U" m# \

该用户从未签到

3#
发表于 2021-3-2 11:26 | 只看该作者
S=Pt/(1+theta(2)+theta(3)+theta(4)+theta(5)+2*(u(1)+u(2)));6 S7 G( g: F, o$ l. V
你确定没写错??你自己文档里面写的分母可是( 1 + m + 2 * x_1 + x_2 )
6 c. `0 ~0 a: B5 o- r. ~& N: {
) S6 S8 C; ?+ k4 H% d然后,这个拟合太耗时,你给的参数上下界太大了,导致随机生成的数据,很多组代进微分方程组去运算起来非常非常慢,至少先去看文献也好查经验公式也好,把数量级、正负号之类的给限制一下,要不然很难做优化。
  h/ d2 F1 ?2 M' f. N: l* E8 j

该用户从未签到

4#
发表于 2021-3-2 15:26 | 只看该作者
来学习一下
您需要登录后才可以回帖 登录 | 注册

本版积分规则

关闭

推荐内容上一条 /1 下一条

EDA365公众号

关于我们|手机版|EDA365电子论坛网 ( 粤ICP备18020198号-1 )

GMT+8, 2025-10-30 23:33 , Processed in 0.140625 second(s), 23 queries , Gzip On.

深圳市墨知创新科技有限公司

地址:深圳市南山区科技生态园2栋A座805 电话:19926409050

快速回复 返回顶部 返回列表