常微分方程组参数求解，利用最小二乘法优化。自己的程序总是索引超过矩阵维度，求助路过的大神施以援手，予以指正，小弟刚入论坛，金币微薄，仅表心意，万分感激。
未知参数仅有三个，但已知参数较多，所以程序看起来有点复杂，程序如下：
function Kinetics1
clear all
clc

k0 = [1.336*10^15  1.955*10^7  3.226*10^4];         % 参数初值
lb = [0  0  0];                   % 参数下限
ub = [+inf  +inf  +inf];    % 参数上限
x0 = [0.97,0.17];
yexp=[0.97 0.17;0.91 0.3;0.87 0.32;0.8 0.35;0.7 0.4;0.6 0.45;0.5 0.48];

% 使用函数lsqnonlin()进行参数估计
[k,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc4,k0,lb,ub,[],x0,yexp);      
ci = nlparci(k,residual,jacobian);
fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
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('  The sum of the squares is: %.1e\n\n',resnorm)
%目标函数
function f = ObjFunc4(k,x0,yexp)
E=zeros(3,1);Kar=zeros(5,1);DH=zeros(5,1);Tref=ones(5,1);theta=zeros(5,1);Eta=zeros(3,1);
% CONSTANTES
% inchangeables divers
Pi=3.14159;
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]
delta=1.5*10^-4;% membrane thickness [m]
L=0.05;% membrane length,[m]
Ac=Pi*r*L;% membrane area [m^2]
Eta(1)=0.01; % Oveall effectiveness factor for r(1)
Eta(2)=0.3; % Oveall effectiveness factor for r(2)
Eta(3)=0.01; % Oveall effectiveness factor for r(3)
Pp=1; % Permeate pressure [bar]
R=8.314; % [J/K.mol]
Tref0=298.15; % Ref temperature for Delta r H calculations [K]
% inchangeables cinetiques
E(1)=240100; % [J/mol]
E(2)=67100; % [J/mol]
E(3)=243900; % [J/mol]
Kar(3)=40.9; % [1/bar]
Kar(5)=2.9e-2; % [1/bar]
Kar(1)=1.8e-1; % [1/bar], Tref=823K
Kar(2)=0.4; % Tref=823K
DH(3)=-70700; % [J/mol]
DH(5)=-82900; % [J/mol]
DH(1)=-38300; % [J/mol]
DH(2)=88700; % [J/mol]
% Gibbs Energy coefficients
Gc=[-75.3 7.6e-2 1.9e-5
-241.7 4.2e-2 7.4e-6
-109.9 -9.2e-2 1.5e-6
-393.4 3.8e-3 1.3e-6];
% Rk: G(H2)=0 because pure specie
% stochiometric coefficients
% 1 2 3
mu=[-1 0 -1 % CH4
-1 -1 -2 % H2O
1 -1 0 % CO
0 1 1 % CO2
3 1 4]; % H2
% Reference Temperature for kinetics parameters, different from Tref0=298K
Tref=648*Tref;
for i=1:2
Tref(i)=823; % [K]
end
Fa=[1.34*10^-5:1.34*10^-5:1.072*10^-4];
kmax=length(Fa);
for k=1:kmax   
Fa0=Fa(k); % CH4 initial feed, [kmol/h]
theta(2)=3; % F0(H2O)/F0(CH4)
theta(3)=0; % F0(CO)/F0(CH4)
theta(4)=0; % F0(CO2)/F0(CH4)
theta(5)=0.001; % F0(H2)/F0(CH4)
Pt=1;% Total pressure [bar]
C1=Ac/Fa0;
tspan = [0.18 : 0.01 : 0.86];
[z x] = ode45(@KineticEqs,tspan,x0,[],k);
x1(k) = real(x(end,1));
x2(k) = real(x(end,2));
y(:,1) = x1;
y(:,2) = x2;
f1 = y(:,1) - yexp(:,1);
f2 = y(:,2) - yexp(:,2);
f = [f1; f2];
end
%常微分方程组
function diff = KineticEqs(z,x,k)
global theta Kar  E DH Gc C1 Pt R mu  Tref  Eta
kr=zeros(3,1);Ka=zeros(5,1);Ke=zeros(3,1);G0=zeros(5,1);k=zeros(3,1);
T=-2.7*10^3*z^2+2.87*10^3*z+37;
%-reaction rate constants
for i=1:3
kr(i)=k(i)*exp(-E(i)/(R*T)); % same units than krr
end
% - Adsorption constants
for i=1:5
Ka(i)=Kar(i)*exp(DH(i)*(1/Tref(i)-1/T)/R); % same units than Kar
end
% - Gibbs energie constants
for i=1:4
G0(i)=Gc(i,1)+Gc(i,2)*T+Gc(i,3)*T^2;
end
DGrxn=zeros(3,1);
for i=1:3
for j=1:4
    DGrxn(i)=DGrxn(i)+mu(j,i)*G0(j);
end
% Equilibrium constants
Ke(i)=exp(-DGrxn(i)/(R*T*1e-3)); % [-] no unit
end
% Expression of the partial pressures, functions of methane and carbon
% dioxide conversions, theta and Pt.
S=Pt/(1+theta(2)+theta(3)+theta(4)+theta(5)+2*x(1));
P(1)=(1-x(1))*S; % [bar]
P(2)=(theta(2)-x(1)-x(2))*S; % [bar]
P(3)=(theta(3)+x(1)-x(2))*S; % [bar]
P(4)=(theta(4)+x(2))*S; % [bar]
P(5)=(theta(5)+3*x(1)+x(2))*S; % [bar]
% Expression of the reaction rates, functions of Kinetic coefficients,
% partial pressures and conversions.
DEN2=(1+Ka(1)*P(1)+Ka(3)*P(3)+Ka(5)*P(5)+Ka(2)*P(2)/P(5))^2;
r(1)=Eta(1)*(kr(1)/(P(5)^2.5))*(P(1)*P(2)-(P(5)^3)*P(3)/Ke(1))/DEN2; %[kmol/(kg.h)]
r(2)=Eta(2)*(kr(2)/P(5))*(P(3)*P(2)-P(5)*P(4)/Ke(2))/DEN2; %[kmol/(kg.h)]
r(3)=Eta(3)*(kr(3)/(P(5)^3.5))*(P(1)*P(2)^2-(P(5)^4)*P(4)/Ke(3))/DEN2; %[kmol/(kg.h)]

% Differential system
diff1=C1*(r(1)+r(3)); % methane conversion [-]
diff2=C1*(r(2)+r(3)); % Carbon dioxide conversion [-]
diff=[diff1 diff2 ]';
如有大神看我程序麻烦，也可以直接看我上传的文件，帮忙写一个程序，必有重谢！+QQ：1550866621

 返回小木虫查看更多