请会matlab的大佬帮我看一下动力学模型的问题

啦啦啦1525 · 发表于 1970-1-1 08:00:00

我是看别人文献发现他用了lhhw的模型，因为我和他的体系相似所以借用了，想编写一下算法

function nihe_0316
% 微分方程的参数估计
%
clear all;
clc;
format long;
global t0;
% k10,E,Ka 对应k(1)、k(2)、k(3)
k0 =[5.33E9;93.06;0.0437]; % 参数初值，重要！，可将结果反复代入迭代计算！！
x0 = 0.63713;  %初始状态，没有就取第一行数据
lb = [0 0 0]; % 参数下限,根据拟合结果，适当调整
ub = [+inf +inf +inf]; % 参数上限
% 实验数据:
expdata= ...
[0       0.63713
60       0.79642
120       0.87377
240       0.90782
360       0.91735];
t0=expdata(:,1);
yexp =expdata(:,2);       % yexp: 对应实验数据[Xb]

% ------------------------------------------------------------------
% 使用函数lsqnonlin()进行参数估计
k=lsqnonlin(@ObjFunc4LNL,k0,lb,ub,[],x0,yexp);

fprintf('\n\n使用函数lsqnonlin()估计得到的参数值为:\n')
fprintf('\tk10 = %.8f\n',k(1))
fprintf('\tE = %.8f\n',k(2))
fprintf('\tKa = %.8f\n',k(3))

% ----------------------------------------------------
tspan = [t0(1) t0(end)];
[t,xx] = ode45(@KineticEqs,tspan,x0,[],k);
x=spline(t,xx,t0); %插值，使理论值和实验值个数一致

% 计算相关系数
A=corrcoef(yexp,x(:,1));
fprintf('\nXb的相关系数R为:\n')
fprintf('\tR = %.8f\n',A(1,2))
%
% 计算决定系数
r(1)=1-sum((yexp-x(:,1)).^2)./sum((yexp-mean(yexp)).^2);
% ********************************决定系数
fprintf('\nXb的决定系数R^2为:\n')
fprintf('\tR^2 = %.8f\n',r(1))

% 计算均方根误差
rmse(1)=sqrt(sum((yexp-x(:,1)).^2)/length(t0));
fprintf('\nXb的均方根误差RMSE为:\n')
fprintf('\tRMSE = %.8f\n',rmse(1))
%
% 计算残差平方和 SSE（Sum of Squares for Error）
ssr(1)=sum((yexp(:,1)-x(:,1)).^2);
fprintf('\nXb的残差平方和为:\n')
fprintf('\tSSE = %.8f\n',ssr(1))
% ----------------------------------------------------
figure(1)
plot(t,xx(:,1),'b-',t0,yexp,'ro'),legend('理论值','实验值','Location','best')
xlabel('时间');ylabel('实验值N');

% ------------------------------------------------------------------
function f = ObjFunc4LNL(k,x0,yexp)
tspan = [t0(1) t0(end)];
[t,xx] = ode45(@KineticEqs,tspan,x0,[],k);
x=spline(t,xx,t0);
f = x - yexp; % 理论值与实验值的差值，残差
end

% ------------------------------------------------------------------
function dXbdt = KineticEqs(t,Xb,k) % 微分方程, k10,E,Ka 对应k(1)、k(2)、k(3)
T=358.15;
m=3;
N0=0.1467;
Cb0=0.5564;
Ca=Cb0*(3-Xb);
Cb=Cb0*(1-Xb);
Cc=Cb0*Xb;
Keq=exp(4873/T-13.813);
Kb=5.2112*k(3);
Kc=2.2695*k(3);
dXbdt = (m*k(1)*exp(-k(2)/8.314*T)*(Ca*Cb-Cc/Keq))/(N0*(1+k(3)*Ca+Kb*Cb+Kc*Cc)^2);
end
% ------------------------------------------------------------------
end

结果算出来的X的计算值不随时间变化，理论上是一个dx/dt=f（x）的方程，我看书上dc/dt=f（c）的方程用RK法可以计算加拟合，不知道为什么我这个就不行，请大家帮我看一下

请会matlab的大佬帮我看一下动力学模型的问题

相关帖子