24小时热门版块排行榜    

CyRhmU.jpeg
查看: 1151  |  回复: 0

安之两隅

金虫 (小有名气)

[求助] Matlab 非线性拟合问题

Matlab小白,参照黄华江matlab在化学工程中的应用的教材来改写模型,但是运行的时候出现了,请大家帮忙看一下,麻烦了!

function kinetics3
clear all
clc
%实验数据
B = [1.319 1.3161 1.3112 1.2892 1.2695 1.2601 1.2485 1.2389 1.2245 1.2052 1.181 1.1456 1.1402 1.1134 1.0827 1.036 0.9634];
C= [0.0813 0.0957 0.1188 0.1532 0.1877 0.2089 0.2333 0.2621 0.2891 0.3323 0.3738 0.4275 0.4464 0.4916 0.5393 0.6141 0.7114];
r = [0.02286 0.02244 0.02202 0.0216 0.02118 0.02076 0.02034 0.01992 0.0195 0.01908 0.01866 0.01796 0.01726 0.01656 0.01586 0.01376 0.01166];

% 用多变量线性回归方法估计动力学参数
R =B./r;
y = R;
X = [ones(size(y)) B C];
b = X\y  % 或[b,bint] = regress(y,X,0.05);
KZ = 1/b(1);
KB = b(2)*KZ;
KC = b(3)*KZ;
% 用lsqnonlin()--求解非线性最小二乘法(非线性数据拟合)问题
beta0 = [KZ KB KC];
lb = [1 1 1];
ub = [+inf +inf +inf];
[beta,resnorm,resid,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFun,beta0,lb,ub,[],B,C,r);
ci = nlparci(beta,resid,jacobian);

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

% 残差关于拟合值的残差图
rc = RateEqs(beta,B,C);
plot(rc,resid,'*');
xlabel('反应速率拟合值');
ylabel('实验速率');
refline(0,0);

% 参数辨识结果
fprintf('\n\nEstimated Parameters:\n')
fprintf('\tKZ= %.4f ± %.4f\n',beta(1),ci(1,2)-beta(1))
fprintf('\tKB = %.2f ± %.2f\n',beta(2),ci(2,2)-beta(2))
fprintf('\tKC = %.2f ± %.2f\n',beta(3),ci(3,2)-beta(3))



% ------------------------------------------------------------------
function f = ObjFun(beta,B,C,r)
rc = RateEqs(beta,B,C);
f = r - rc;

% ------------------------------------------------------------------
function rc = RateEqs(beta,B,C)    % Rate equation
rc = beta(1)*B./(1+beta(2)*B+beta(3)*C);



错误使用 snls (line 47)
Objective function is returning undefined values at initial point. lsqnonlin cannot continue.

出错 lsqncommon (line 167)
            snls(funfcn,xC,lb,ub,flags.verbosity,options,defaultopt,initVals.F,initVals.J,caller, ...

出错 lsqnonlin (line 253)
   lsqncommon(funfcn,xCurrent,lb,ub,options,defaultopt,allDefaultOpts,caller,...

出错 Kinetics3 (line 23)
    lsqnonlin(@ObjFun,beta0,lb,ub,[],B,C,r);
已阅   回复此楼   关注TA 给TA发消息 送TA红花 TA的回帖
相关版块跳转 我要订阅楼主 安之两隅 的主题更新
信息提示
请填处理意见