当前位置: 首页 > 计算模拟 >Matlab 用多变量线性回归方法估计动力学参数

Matlab 用多变量线性回归方法估计动力学参数

作者 l1003785517
来源: 小木虫 200 4 举报帖子
+关注

各位虫友好!
       我在用多变量线性回归方法估计动力学参数的时候写了如下代码(根据书上例题进行修改的),但是程序运行不了,说是目标函数返回了未定义参数。
有哪位帮忙看一下是哪里写错了吗?我自己一直没检查出来哪里错了,谢谢大家!

function KineticsEst
% 动力学参数估计

clear all
clc
C_DC=[6.31317 6.29412 6.27507 6.26237 6.24713 6.23316 6.221095 6.210935];
C_TT=[0.029302 0.04459 0.059878 0.069433 0.080899 0.091091 0.099372 0.107016];
C_NP=[0.007644 0.011466 0.015288 0.018473 0.022295 0.026117 0.029939 0.032487];
r1=[0.0574 0.0407 0.0241 0.0278 0.0315 0.0268 0.022 0.0172];
% 用多变量线性回归方法估计动力学参数
R = C_DC/r1;
y = R;
X = [ones(size(y)) C_DC C_TT C_NP];
b=X\y;    % 或b = X\y
k1= 1/b(1);
k2 = b(2)*k1;
k3 = b(3)*k1;
k4 = b(4)*k1;

% 用lsqnonlin()--求解非线性最小二乘法(非线性数据拟合)问题
beta0 = [k1 k2 k3 k4];
lb = [0  0  0  0];
ub = [+inf  +inf  +inf  +inf];
[beta,resnorm,resid,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFun,beta0,lb,ub,[], C_DC,C_TT,C_NP,r1);
ci = nlparci(beta,resid,jacobian);

% 模型适定性判别
Ne = length(r1);
Np = length(beta);
[rho2, F] = rho2_F(k1, r1, resnorm, Ne, Np);

% 残差关于拟合值的残差图
rc = RateEqs(beta,C_DC,C_TT,C_NP);
plot(rc,resid,'*')
xlabel('反应速率拟合值, mol/(min g催化剂)')
ylabel('残差R, mol/(min g催化剂)')
refline(0,0)

% 参数辨识结果
fprintf('\n\nEstimated Parameters:\n')
fprintf('\tk1 = %.4f ± %.4f\n',beta(1),ci(1,2)-beta(1))
fprintf('\tk3 = %.4f ± %.4f\n',beta(2),ci(2,2)-beta(2))
fprintf('\tk4 = %.4f ± %.4f\n',beta(3),ci(3,2)-beta(3))
fprintf('\tk5 = %.4f ± %.4f\n',beta(4),ci(4,2)-beta(4))
fprintf('  决定性指标ρ^2: %.3f\n',rho2)
fprintf('  F比: %.3f\n\n',F)


% ------------------------------------------------------------------
function f = ObjFun(beta,C_DC,C_TT,C_NP,r1)
rc =RateEqs(beta,C_DC,C_TT,C_NP);
f = r1 - rc;

% ------------------------------------------------------------------
function rc =RateEqs(beta,C_DC,C_TT,C_NP)    % Rate equation
rc = beta(1)*C_DC./(1+beta(2)*C_DC+beta(3)*C_TT+beta(4)*C_NP); 返回小木虫查看更多

今日热帖
  • 精华评论
  • hzlhm

    你这里还却  rho2_F( )自定义 函数,

猜你喜欢
下载小木虫APP
与700万科研达人随时交流
  • 二维码
  • IOS
  • 安卓