EDA365电子论坛网
标题:
常微分方程组参数拟合问题
[打印本页]
作者:
pTDbn25
时间:
2021-3-2 10:28
标题:
常微分方程组参数拟合问题
代码如下;
5 \5 c* P: u8 R% ~
function Kinetics4
8 {9 }& t0 w( k0 D! Y8 u1 S4 X
clear all
" s1 v" J/ z: O+ \
clc
; k4 `8 V$ z6 Y- f7 ?1 G
k0 = [2.159*10^8 4.734*10^6 174686 149892 0.99 0.7 0.5 0.78];%参数初值，其实是不知道的
' `+ N1 h7 U( {, P+ {8 A
lb = [-inf -inf -inf -inf -inf -inf -inf -inf]; % 参数下限
- T8 b6 r* \( K
ub = [inf inf inf inf inf inf inf inf]; % 参数上限
: s& N9 z; Q7 B
u0 = [0,0]; %y初值
" [8 }0 m1 S& O
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];
$ E3 Y+ I2 X; h( 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];
/ {; |3 J- y9 h, h* ^, n8 i" C
yexp=[a,b];
: Y) _: L4 K* M! E3 ~! O% \ g
% 使用函数lsqnonlin()进行参数估计
9 A9 a2 ^7 j b8 Q/ g1 ~( Y' g
; f# J. b6 h, g+ G- n4 {5 R5 L2 a
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
/ z- Q. U2 D8 F1 w
lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[],u0,yexp);
8 i* I$ m( M6 T
ci = nlparci(k,residual,jacobian);
" r1 I% r; D8 M; ^8 u I
fprintf('\tk1 = %.4f ± %.4f\n',k(1),ci(1,2)-k(1))
6 w% W0 s& P; T3 t) C& L
fprintf('\tk2 = %.4f ± %.4f\n',k(2),ci(2,2)-k(2))
9 a$ ~+ E7 |4 g- v7 K1 x
fprintf('\tk3 = %.4f ± %.4f\n',k(3),ci(3,2)-k(3))
8 s4 P5 x, [5 ]( ^5 P0 ]! A5 Q
fprintf('\tk4 = %.4f ± %.4f\n',k(4),ci(4,2)-k(4))
. y6 S# t8 {& t/ m5 m
fprintf('\tk5 = %.4f ± %.4f\n',k(5),ci(5,2)-k(5))
2 M& d+ m1 S; X* t
fprintf('\tk6 = %.4f ± %.4f\n',k(6),ci(6,2)-k(6))
' l4 F% O6 }; J
fprintf('\tk7 = %.4f ± %.4f\n',k(7),ci(7,2)-k(7))
0 V; u- Y, M8 D1 l2 w0 E
fprintf('\tk8 = %.4f ± %.4f\n',k(8),ci(8,2)-k(8))
, D; y4 m" M* A9 f. y* V
fprintf(' The sum of the squares is: %.1e\n\n',resnorm)
$ z% I/ B6 |- f5 e: J" R5 W
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
7 x: ` d5 R& u* T7 a4 V$ n
%目标函数
" J0 w! I+ n, G9 r4 e4 z2 E
function f =ObjFunc4LNL(k,u0,yexp)
& f9 e: c/ r" K1 e. O
global theta Pt Ac R T Fa0
2 z" H) `/ |8 l" x4 S* b: W
theta=zeros(5,1);
, C) ?/ G/ J! h
r0= 1.24*10^-3; % Outer diameter of the membrane, [m]
; `& R3 v4 a) C% G3 J3 I
ri= 9.4*10^-4;% Inner diameter of the membrane, [m]
2 C: y0 B2 f5 d( T- i
r=(r0-ri)/log(r0/ri);% equivalent diameter ,[m]
2 U: }" {( ~8 ^8 h, W2 T6 `
Pi=3.14159;
$ ^! L2 C. o& D# v+ c1 h# v% ?
L=0.05;% membrane length,[m]
4 H2 r- Q" w( e; r; a* m' c- Y" Y
Ac=Pi*r*L;% membrane area [m^2]
# \% X, l9 E: X2 u- j7 I$ G
R=8.314; % [J/K.mol]
. h2 y$ w6 ^) T9 z9 R
Pt=101.325;
! J- A9 C- S3 `! Q5 m4 u& d6 X" J
Ftot=[2.678*10^-2;2.678*10^-2*1.5;2.678*10^-2*2];
* y" V8 E- Z: E
Kmax=length(Ftot);
2 Y- _' y) I& c" y
7 B3 w( {. j- S
T0=600+273;% Inlet temperature [K]
% U, [; n0 J$ T0 T5 k; _* L
nmax=9;
% K8 u- Z+ a7 ~! ]$ h: h f
X1=zeros(nmax,Kmax);X2=zeros(nmax,Kmax);
% |& Z E( O# e* u1 E9 Z
for K=1:Kmax
+ A4 s! \# @# @5 h; N% D
Fa0=Ftot(K); % If Pt is the parameter/Matrix
6 {2 _4 B9 z) p% t
%Pt=Pt0+(k-1)*1; % If Pt is an axis
" M% W: q! W" p: ^
%L=Lsize(k); % If L is the parameter/Matrix
( t4 x9 M* H' ?
%theta(2)=H2Oratio+(k-1)*1; % If m is the parameter/Matrix
. ?1 f" m3 y# C& e
theta(2)=3; % F0(H2O)/F0(CH4)
8 n- f/ d; R3 G4 G$ B
theta(3)=0; % F0(CO)/F0(CH4)
+ f$ a! I7 x- V/ O) O0 l
theta(4)=0; % F0(CO2)/F0(CH4)
5 H4 n* F- O& T" J; u# w! G
theta(5)=0; % F0(H2)/F0(CH4)
: O( M$ K4 d5 L! C3 E7 h0 v, G
for n=1:nmax
* j& D% `4 d* A+ z
%delta=delta0+(n-1)*1e-6; % If delta is the parameter/Matrix
3 U, |, s3 F$ v4 B5 i
%rit=r0+delta; % inner tube radius [m]
7 G& a. n" e, H. a4 o
%sweep=s0+(n-1); % If sweep is the parameter/Matrix
6 P% N2 `: ~ P% D
%T=T0; % If T is a constant
$ x4 I$ J2 q ^/ l& W
T=T0+(n-1)*50;
/ ^+ Y9 F; i" f0 X
ksispan=(0:0.1:1); % precision
7 \" O; s4 x0 }" O' _
[ksi,S1]=ode23s(@KineticEqs,ksispan,u0,[],k); % PBR reactor
5 e1 G2 Z" w8 V( E1 f* S9 V4 r: U
X1(n,Kmax)=real(S1(end,1)); % methane finale conversion
$ N) r! T! t% I2 }+ P6 w1 `9 u
X2(n,Kmax)=real(S1(end,2)); % carbon dioxide finale conversion
2 i; I6 o! c, Z' D2 X
end
$ k( Y4 a+ E8 o. T! l5 {: {
end
5 u" O. J. ]" `8 y0 }0 f' h/ Y
f1=X1(n,Kmax)-yexp(:,1);
! |+ k6 @6 x+ h" U1 X: F7 ~
f2=X2(n,Kmax)-yexp(:,2);
2 | n; U/ D, ?
f=[f1;f2];
; V" c" J- G% g
%T0=T0-273.15;T=T-273.15; % Temperature expressed in Celsius
, A; H6 r! C% J
%Tspan=(T0:1:T);
, C4 f5 e8 |, d$ y
%Pspan=(Pt0:1
t);
/ u, X" A9 s# m* L! Q: N# z) [* F
%Dspan=(delta0*1e6:1:delta*1e6);
$ a# F& k, y5 d9 T: O
%Sspan=(s0:1:sweep);
: ^4 p. }% A+ E
%D(:,
=X3(:,
-X1(:,
; % Delta = XCH4(MR) - XCH4(PBR)
! s* [4 o/ y, j9 K6 d+ H) s R; T
# e- y# A8 F- F* b2 s0 r
function diff=KineticEqs(ksi,u,k)
' A% M; V4 \1 ~7 n$ i) v/ d
global theta Pt Ac Fa0 T R
! a. l1 e8 ~; j/ ?
S=Pt/(1+theta(2)+theta(3)+theta(4)+theta(5)+2*(u(1)+u(2)));
$ q- Q# R( Z6 s7 F4 }4 ^. o
P(1)=(1-u(1)-u(2))*S;
: |9 O* C1 c- G; u# T% E
P(2)=(theta(2)-u(1)-2*u(2))*S;
2 M3 c6 @0 A t6 h% Z- [
P(3)=(theta(3)+u(1))*S;
5 u4 J, u* l5 L6 c
P(4)=(theta(4)+u(2))*S;
; s3 k3 i. N" u5 E8 a2 t# c a i
P(5)=(theta(5)+3*(u(1)+u(2))+u(2))*S;
/ f3 N* U- M W+ q9 |
% Differential System
' [( Z$ M' D2 C5 y9 o
diff1=Ac/Fa0*k(1)*exp(-k(3)/R/T)*P(1)^k(5)*P(2)^k(7);
' l0 ~2 `" i# U
diff2=Ac/Fa0*k(2)*exp(-k(4)/R/T)*P(1)^k(6)*P(2)^k(8);
/ m: t3 U7 s( P1 F3 i# O
diff=[diff1 diff2]';
' w. |. o1 O; l% D5 s5 \
4 j3 Y! o. j5 \) M9 Y/ j$ P" B
运行结果：
7 J# p7 Q8 u0 V( @% g4 u
In nlparci (line 104)
( A, Y H `* j& p. U% R5 H. k
In Kinetics4 (line 15)
4 w$ ^$ ^* _( z
警告: 矩阵为奇异工作精度。
8 }1 t, }$ M0 L) M s
k1 = 215900000.0000 ± NaN
6 H* i8 C- y" ^ y8 R, c5 A4 {
k2 = 4734000.0000 ± NaN
" n3 _$ \! `4 |& j' y
k3 = 174686.0000 ± NaN
# ^1 `; X4 Z7 p3 w
k4 = 149892.0000 ± NaN
6 c- U) ~, `; s* x# E/ m# H, a6 j
k5 = 1.2381 ± NaN
$ @; O( a( p* T# J* c: r4 G6 t# L
k6 = 0.5426 ± NaN
$ q& q0 c+ E; r6 g- u
k7 = -0.3426 ± NaN
6 T- T; \3 T- ~2 Q
k8 = 0.4455 ± NaN
& c8 i3 }# u+ g8 c5 d) E
The sum of the squares is: 2.0e+00
3 `$ K7 @* h5 S3 y5 ]: o
作者:
cichishia
时间:
2021-3-2 10:48
S=Pt/(1+theta(2)+theta(3)+theta(4)+theta(5)+2*(u(1)+u(2)));
( z) r; G; S' G; V* w( s. x3 O& R
你确定没写错？？你自己文档里面写的分母可是( 1 + m + 2 * x_1 + x_2 )
2 X' U# c$ K2 |4 |& J0 f" j( d
( q$ X9 e" c" B& A" K1 k; D1 d( O
然后，这个拟合太耗时，你给的参数上下界太大了，导致随机生成的数据，很多组代进微分方程组去运算起来非常非常慢，至少先去看文献也好查经验公式也好，把数量级、正负号之类的给限制一下，要不然很难做优化。
8 j$ p8 R% h9 e2 b
4 \7 ~8 t0 l0 H7 ~
你给的原始参数的估计值的效果
; Q+ P& U( h. ^ g
/ ^6 k, k: j- `1 U
[attach]317568[/attach]
: l! |- t, u% U* d2 a. C
+ S0 `+ k* `( u/ X* {2 }
拟合得到的两组参数的效果1
, D9 V6 V$ [$ X/ s9 [8 `
[attach]317566[/attach]
% t% @+ |! R0 z) ^
2 U7 v- L2 D, }- X- \# j
拟合得到的两组参数的效果2
% S8 G5 }" ]9 }6 c7 R* A
[attach]317567[/attach]
作者:
shuddkk
时间:
2021-3-2 11:26
S=Pt/(1+theta(2)+theta(3)+theta(4)+theta(5)+2*(u(1)+u(2)));
! ~5 L& C i- B8 @5 Y
你确定没写错？？你自己文档里面写的分母可是( 1 + m + 2 * x_1 + x_2 )
' L( f4 v* I9 Z9 M" x" p, x# x
( `1 ^/ Z, X. J4 D4 i& j# a
然后，这个拟合太耗时，你给的参数上下界太大了，导致随机生成的数据，很多组代进微分方程组去运算起来非常非常慢，至少先去看文献也好查经验公式也好，把数量级、正负号之类的给限制一下，要不然很难做优化。
! z$ l% P: v) R2 h
作者:
zaiyiaaaa
时间:
2021-3-2 15:26
来学习一下
欢迎光临 EDA365电子论坛网 (https://bbs.eda365.com/)
Powered by Discuz! X3.2
t);/ u, X" A9 s# m* L! Q: N# z) [* F
=X3(:,-X1(:,; % Delta = XCH4(MR) - XCH4(PBR)