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

pTDbn25 · 发表于 2021-3-2 10:28

EDA365欢迎您登录！

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

x

代码如下;
function Kinetics4
clear all
clc
k0 = [2.159*10^8 4.734*10^6 174686 149892 0.99 0.7 0.5 0.78];%参数初值，其实是不知道的
lb = [-inf -inf -inf -inf -inf -inf -inf -inf];                % 参数下限
ub = [inf inf inf inf inf inf inf inf]; % 参数上限
u0 = [0,0]; %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];
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];
yexp=[a,b];
% 使用函数lsqnonlin()进行参数估计

[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[],u0,yexp);
ci = nlparci(k,residual,jacobian);
fprintf('\tk1 = %.4f ± %.4f\n',k(1),ci(1,2)-k(1))
fprintf('\tk2 = %.4f ± %.4f\n',k(2),ci(2,2)-k(2))
fprintf('\tk3 = %.4f ± %.4f\n',k(3),ci(3,2)-k(3))
fprintf('\tk4 = %.4f ± %.4f\n',k(4),ci(4,2)-k(4))
fprintf('\tk5 = %.4f ± %.4f\n',k(5),ci(5,2)-k(5))
fprintf('\tk6 = %.4f ± %.4f\n',k(6),ci(6,2)-k(6))
fprintf('\tk7 = %.4f ± %.4f\n',k(7),ci(7,2)-k(7))
fprintf('\tk8 = %.4f ± %.4f\n',k(8),ci(8,2)-k(8))
fprintf('  The sum of the squares is: %.1e\n\n',resnorm)
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
%目标函数
function f =ObjFunc4LNL(k,u0,yexp)
global theta Pt Ac R T Fa0
theta=zeros(5,1);
r0= 1.24*10^-3; % Outer diameter  of the membrane, [m]
ri= 9.4*10^-4;% Inner diameter  of the membrane, [m]
r=(r0-ri)/log(r0/ri);% equivalent diameter ,[m]
Pi=3.14159;
L=0.05;% membrane length,[m]
Ac=Pi*r*L;% membrane area [m^2]
R=8.314; % [J/K.mol]
Pt=101.325;
Ftot=[2.678*10^-2;2.678*10^-2*1.5;2.678*10^-2*2];
Kmax=length(Ftot);

T0=600+273;% Inlet temperature [K]
nmax=9;
X1=zeros(nmax,Kmax);X2=zeros(nmax,Kmax);
for K=1:Kmax
Fa0=Ftot(K); % If Pt is the parameter/Matrix
%Pt=Pt0+(k-1)*1; % If Pt is an axis
%L=Lsize(k); % If L is the parameter/Matrix
%theta(2)=H2Oratio+(k-1)*1; % If m is the parameter/Matrix
theta(2)=3; % F0(H2O)/F0(CH4)
theta(3)=0; % F0(CO)/F0(CH4)
theta(4)=0; % F0(CO2)/F0(CH4)
theta(5)=0; % F0(H2)/F0(CH4)
for n=1:nmax
%delta=delta0+(n-1)*1e-6; % If delta is the parameter/Matrix
%rit=r0+delta; % inner tube radius [m]
%sweep=s0+(n-1); % If sweep is the parameter/Matrix
%T=T0; % If T is a constant
T=T0+(n-1)*50;
ksispan=(0:0.1:1); % precision
[ksi,S1]=ode23s(@KineticEqs,ksispan,u0,[],k); % PBR reactor
X1(n,Kmax)=real(S1(end,1)); % methane finale conversion
X2(n,Kmax)=real(S1(end,2)); % carbon dioxide finale conversion
end
end
f1=X1(n,Kmax)-yexp(:,1);
f2=X2(n,Kmax)-yexp(:,2);
f=[f1;f2];
%T0=T0-273.15;T=T-273.15; % Temperature expressed in Celsius
%Tspan=(T0:1:T);
%Pspan=(Pt0:1

t);
%DSPan=(delta0*1e6:1:delta*1e6);
%Sspan=(s0:1:sweep);
%D(:,

=X3(:,

-X1(:,

; % Delta = XCH4(MR) - XCH4(PBR)

function diff=KineticEqs(ksi,u,k)
global theta  Pt  Ac Fa0 T R
S=Pt/(1+theta(2)+theta(3)+theta(4)+theta(5)+2*(u(1)+u(2)));
P(1)=(1-u(1)-u(2))*S;
P(2)=(theta(2)-u(1)-2*u(2))*S;
P(3)=(theta(3)+u(1))*S;
P(4)=(theta(4)+u(2))*S;
P(5)=(theta(5)+3*(u(1)+u(2))+u(2))*S;
% Differential System
diff1=Ac/Fa0*k(1)*exp(-k(3)/R/T)*P(1)^k(5)*P(2)^k(7);
diff2=Ac/Fa0*k(2)*exp(-k(4)/R/T)*P(1)^k(6)*P(2)^k(8);
diff=[diff1 diff2]';

运行结果：
In nlparci (line 104)
  In Kinetics4 (line 15)
警告: 矩阵为奇异工作精度。
      k1 = 215900000.0000 ± NaN
      k2 = 4734000.0000 ± NaN
      k3 = 174686.0000 ± NaN
      k4 = 149892.0000 ± NaN
      k5 = 1.2381 ± NaN
      k6 = 0.5426 ± NaN
      k7 = -0.3426 ± NaN
      k8 = 0.4455 ± NaN
  The sum of the squares is: 2.0e+00

zaiyiaaaa · 发表于 2021-3-2 15:26

来学习一下

shuddkk · 发表于 2021-3-2 11:26

S=Pt/(1+theta(2)+theta(3)+theta(4)+theta(5)+2*(u(1)+u(2)));
你确定没写错？？你自己文档里面写的分母可是( 1 + m + 2 * x_1 + x_2 )

然后，这个拟合太耗时，你给的参数上下界太大了，导致随机生成的数据，很多组代进微分方程组去运算起来非常非常慢，至少先去看文献也好查经验公式也好，把数量级、正负号之类的给限制一下，要不然很难做优化。

cichishia · 发表于 2021-3-2 10:48

S=Pt/(1+theta(2)+theta(3)+theta(4)+theta(5)+2*(u(1)+u(2)));
你确定没写错？？你自己文档里面写的分母可是( 1 + m + 2 * x_1 + x_2 )

然后，这个拟合太耗时，你给的参数上下界太大了，导致随机生成的数据，很多组代进微分方程组去运算起来非常非常慢，至少先去看文献也好查经验公式也好，把数量级、正负号之类的给限制一下，要不然很难做优化。

你给的原始参数的估计值的效果

拟合得到的两组参数的效果1

拟合得到的两组参数的效果2

帐号		自动登录	找回密码
密码			注册

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

EDA365欢迎您登录！

推荐内容 /1