|
EDA365欢迎您登录!
您需要 登录 才可以下载或查看,没有帐号?注册
x
代码如下;
4 e z5 I$ F! Z- d Y$ l$ ~1 Ffunction Kinetics43 J3 T: i a# k. P, E$ f( `1 z
clear all
8 B* v# C9 ]2 U' P% a y& p( P( K. O- Dclc
+ C4 v6 e( B0 s- w6 W1 f+ Yk0 = [2.159*10^8 4.734*10^6 174686 149892 0.99 0.7 0.5 0.78];%参数初值,其实是不知道的
1 U8 W" }7 _) h3 E, jlb = [-inf -inf -inf -inf -inf -inf -inf -inf]; % 参数下限
" Z" @( r6 m: U0 {8 Nub = [inf inf inf inf inf inf inf inf]; % 参数上限6 M3 Y2 _# u) {( ^! O' p0 h/ f
u0 = [0,0]; %y初值
; R* B+ A9 \5 y3 B+ ~8 Xa=[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];: V% Z6 x: V+ L5 O" U
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];
! m( H+ Q- W5 Q I- }+ X, Byexp=[a,b];
8 G1 G% ?& q+ r+ `# N1 t% 使用函数lsqnonlin()进行参数估计% E/ U; Q' A4 i/ k4 E
. e' y. w4 v: T: [
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...9 s+ n% T F" }1 V' ?5 ?
lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[],u0,yexp);
0 l% w* {, E. Ici = nlparci(k,residual,jacobian);
* A4 ]; ^% |' V- B3 Jfprintf('\tk1 = %.4f ± %.4f\n',k(1),ci(1,2)-k(1))! t3 D* P* e& n2 u$ e
fprintf('\tk2 = %.4f ± %.4f\n',k(2),ci(2,2)-k(2))
" k# h% f, |+ g+ Xfprintf('\tk3 = %.4f ± %.4f\n',k(3),ci(3,2)-k(3)) 4 I8 w! `3 L4 M
fprintf('\tk4 = %.4f ± %.4f\n',k(4),ci(4,2)-k(4))
. @! K L% x" G2 e" O+ }5 nfprintf('\tk5 = %.4f ± %.4f\n',k(5),ci(5,2)-k(5))
" K9 O, L% v1 e+ T. [. Y1 k0 Ufprintf('\tk6 = %.4f ± %.4f\n',k(6),ci(6,2)-k(6))$ {) _5 I, N* L& x8 m+ |& I" i3 M& J
fprintf('\tk7 = %.4f ± %.4f\n',k(7),ci(7,2)-k(7))
* @% B, ~* ~$ E& D. Gfprintf('\tk8 = %.4f ± %.4f\n',k(8),ci(8,2)-k(8))
" U% c, P1 A+ |9 G3 K* cfprintf(' The sum of the squares is: %.1e\n\n',resnorm)
( p& z+ G- n) L+ G% G" zfprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')$ y: r5 a/ X5 S2 L% w
%目标函数
( w/ g* u: f$ v5 H% D0 ~- _. ]function f =ObjFunc4LNL(k,u0,yexp)0 I1 z8 E1 { q3 s( z( O/ T8 {2 Y
global theta Pt Ac R T Fa0
6 O D, O1 U' Ltheta=zeros(5,1);
( A, n& K4 x6 p4 M0 r" cr0= 1.24*10^-3; % Outer diameter of the membrane, [m]) y: R" j% a. c9 \$ ~$ |9 z1 h
ri= 9.4*10^-4;% Inner diameter of the membrane, [m]" s7 k4 D4 I# b v
r=(r0-ri)/log(r0/ri);% equivalent diameter ,[m]
' L% {) M9 p) v2 C+ `Pi=3.14159;
9 D9 p1 d7 V4 e( F/ |3 x: Y( bL=0.05;% membrane length,[m]
! U7 H2 x5 P* J& ^' Q, ?1 dAc=Pi*r*L;% membrane area [m^2]
9 H7 o' a" d2 d. ]+ F3 eR=8.314; % [J/K.mol]+ p8 d. F, j' q
Pt=101.325;
2 ^1 }6 `) v+ x+ q8 X+ |Ftot=[2.678*10^-2;2.678*10^-2*1.5;2.678*10^-2*2];
# `2 j' z$ ?( q6 x) ~4 e7 T2 R x* jKmax=length(Ftot);+ Z( w2 O2 h b7 K& e) e4 C
# y5 {8 g9 m" }3 ~8 ]# h, H
T0=600+273;% Inlet temperature [K]
% v' [; K9 L* k! [nmax=9;% Z+ D% O. N/ Y3 e5 @
X1=zeros(nmax,Kmax);X2=zeros(nmax,Kmax);
2 B9 V1 G8 G7 Qfor K=1:Kmax
/ p5 b$ q" k5 B- bFa0=Ftot(K); % If Pt is the parameter/Matrix! F3 j* l' O3 `# H1 I$ u
%Pt=Pt0+(k-1)*1; % If Pt is an axis _; B! ^0 p* b8 }) F6 d2 U. x8 S1 O
%L=Lsize(k); % If L is the parameter/Matrix
! h/ m4 B" L) m4 q9 o/ H( T%theta(2)=H2Oratio+(k-1)*1; % If m is the parameter/Matrix0 g5 d. M& g4 [7 t0 F
theta(2)=3; % F0(H2O)/F0(CH4)
( F J! s- v" N( X6 Jtheta(3)=0; % F0(CO)/F0(CH4)9 K8 o, K, [& V
theta(4)=0; % F0(CO2)/F0(CH4)
+ Y( d. R( z) ?. T6 btheta(5)=0; % F0(H2)/F0(CH4)
+ ]; D1 H t$ c3 M' nfor n=1:nmax! c3 s$ ?& {1 ^" w
%delta=delta0+(n-1)*1e-6; % If delta is the parameter/Matrix) y- ^' @. U8 _8 z7 ^6 D' A
%rit=r0+delta; % inner tube radius [m]3 r2 m& ^+ ~1 q* I5 D- l% B
%sweep=s0+(n-1); % If sweep is the parameter/Matrix
& }8 O) a- \* L9 j9 _5 ?2 H' R%T=T0; % If T is a constant/ u9 I' F/ d/ \ Y' Y( e8 G
T=T0+(n-1)*50;0 I: [9 d. f* z% o. |
ksispan=(0:0.1:1); % precision G: S% ^3 y* F2 S0 J" X+ S1 m
[ksi,S1]=ode23s(@KineticEqs,ksispan,u0,[],k); % PBR reactor
+ w9 r- l$ l" T0 I& k) T0 IX1(n,Kmax)=real(S1(end,1)); % methane finale conversion6 q% w! w- f3 o' W% t/ G% L: Q
X2(n,Kmax)=real(S1(end,2)); % carbon dioxide finale conversion
. @0 ~: F5 z) u) e- dend$ v" I2 I* A5 ^" h; X
end
9 |3 b$ ]& X; R" `# U6 C) r* Uf1=X1(n,Kmax)-yexp(:,1);
# q: o; J. l+ J, _f2=X2(n,Kmax)-yexp(:,2);
T6 G$ z( Z8 Ff=[f1;f2];- R l! L7 X# }/ q
%T0=T0-273.15;T=T-273.15; % Temperature expressed in Celsius 4 m% j/ v2 [" ^9 G. d2 k" w
%Tspan=(T0:1:T);9 m. k4 r& B+ I y( |
%Pspan=(Pt0:1 t);
( |) b3 @6 U& K6 L%DSPan=(delta0*1e6:1:delta*1e6);
/ J% b" t( U' G%Sspan=(s0:1:sweep);
1 N: q$ g# J: z6 K: I6 Z3 |' V9 S%D(:, =X3(:, -X1(:, ; % Delta = XCH4(MR) - XCH4(PBR)
& T( E* H9 _6 k8 K
2 e0 O* C9 q r4 Y& }$ s( Lfunction diff=KineticEqs(ksi,u,k)) ^$ b( Q! S( ~% h7 V2 e3 n& x# M
global theta Pt Ac Fa0 T R. R* s# |" H0 j6 r$ f' c
S=Pt/(1+theta(2)+theta(3)+theta(4)+theta(5)+2*(u(1)+u(2)));$ F9 t: ^1 Q; h+ w7 M h# H8 T0 Q5 N
P(1)=(1-u(1)-u(2))*S;) V5 ^- Y8 k/ i. F
P(2)=(theta(2)-u(1)-2*u(2))*S;3 I$ n8 \$ I& B. l, u- m7 R# D7 B" F
P(3)=(theta(3)+u(1))*S;
7 s; m7 M7 E' lP(4)=(theta(4)+u(2))*S;
' M# g5 Y9 m' v: U* B( m: j3 p7 yP(5)=(theta(5)+3*(u(1)+u(2))+u(2))*S;
/ V. O! d/ x* R% Differential System# V6 U" n, |+ i+ A" b
diff1=Ac/Fa0*k(1)*exp(-k(3)/R/T)*P(1)^k(5)*P(2)^k(7);
e* t9 V: x: |0 }2 p7 l" b3 |( Ndiff2=Ac/Fa0*k(2)*exp(-k(4)/R/T)*P(1)^k(6)*P(2)^k(8);
$ r" Q. p2 W3 _. m0 ldiff=[diff1 diff2]';" F# h; J$ B; T6 ]; H6 n
7 u2 e& e- L# D& [
运行结果:0 ~: B; V$ R* }
In nlparci (line 104)
7 H% v- q7 C% X3 {8 n8 p, g In Kinetics4 (line 15)8 _$ R4 F2 A3 k
警告: 矩阵为奇异工作精度。* q2 y5 t' A0 X- P* t
k1 = 215900000.0000 ± NaN5 W) ^6 Y6 v4 H, o' Y8 C0 u
k2 = 4734000.0000 ± NaN$ m* e' b, Q4 d9 S7 e# Y
k3 = 174686.0000 ± NaN! g- ^% F4 I0 x4 r( W/ B: U
k4 = 149892.0000 ± NaN. G3 R& H4 a7 A* ^* |4 B
k5 = 1.2381 ± NaN
9 T& \; N# r1 X# w6 z k6 = 0.5426 ± NaN% B1 N1 A: u+ K/ k, r8 }; u
k7 = -0.3426 ± NaN
2 h6 N. g+ Q* M: V' Y. @; d, n k8 = 0.4455 ± NaN* |) f+ e" H# K2 T% _ s# t
The sum of the squares is: 2.0e+00
3 f9 h/ O! `: K. [* z/ X n: x+ s. S |
|