kineticsest_ident2.m

本代码为黄华江编著《实用化工计算机模拟—MATLAB在化学工程中的应用》的配套车程序

function KineticsEst_Ident2
% 动力学参数估计及模型辨识---判别四个模型中哪个模型最好 
%
%   Author: HUANG Huajiang
%   Copyright 2003 UNILAB Research Center, 
%   East China University of Science and Technology, Shanghai, PRC
%   $Revision: 1.0 $  $Date: 2003/03/23 $
%
%   [Ref] H. Scott Fogler, Elements of Chemical Reaction Engineering, Prentice Hall PTR, 
%   Upper Saddle River, New Jersey, 1999

clear all
clc

global rA Pe Pea Ph2

rA = [1.04 3.13 5.21 3.82 4.19 2.391 3.867 2.199 0.75];  
Pe = [1 1 1 3 5 0.5 0.5 0.5 0.5];
Pea = [1 1 1 1 1 1 0.5 3 5];
Ph2 = [1 3 5 3 3 3 5 3 1]; 

IDmodel = 'a';
x0 = [1 1];
[x,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc,x0,zeros(size(x0)),inf*ones(size(x0)),[],IDmodel);
plot(residual,'bo')
xlabel('实验序号')
refline(0,0)
ylabel('反应速率残差，mol/(kg cat s)')
title('模型(a)的残差图')
refline(0,0)
fprintf('\n Model (a):\n')
fprintf('  k = %.3f  KE = %.3f\n',x(1),x(2))
fprintf('  The sum of the squares is: %.3f\n\n',resnorm)

IDmodel = 'b';
x0 = [1 1 1];
[x,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc,x0,zeros(size(x0)),inf*ones(size(x0)),[],IDmodel);
figure
plot(residual,'bo')
xlabel('实验序号')
ylabel('反应速率残差，mol/(kg cat s)')
title('模型(b)的残差图')
refline(0,0)
fprintf('\n Model (b):\n')
fprintf('  k = %.3f  KE = %.3f  KA = %.3f\n',x(1),x(2),x(3))
fprintf('  The sum of the squares is: %.3f\n\n',resnorm)

IDmodel = 'c';
x0 = [1 1];  
[x,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc,x0,zeros(size(x0)),inf*ones(size(x0)),[],IDmodel);
figure
plot(residual,'bo')
xlabel('实验序号')
ylabel('反应速率残差，mol/(kg cat s)')
title('模型(c)的残差图')
refline(0,0)
fprintf('\n Model (c):\n')
fprintf('  k = %.3f  KE = %.3f\n',x(1),x(2))
fprintf('  The sum of the squares is: %.3f\n\n',resnorm)

IDmodel = 'd';
x0 = [1 1 1];  
[x,resnorm,residual,exitflag,output,lambda,jacobian] = ...
    lsqnonlin(@ObjFunc,x0,zeros(size(x0)),inf*ones(size(x0)),[],IDmodel);
figure
plot(residual,'bo')
xlabel('实验序号')
ylabel('反应速率残差，mol/(kg cat s)')
title('模型(d)的残差图')
refline(0,0)
fprintf('\n Model (d):\n')
fprintf('  k = %.3f  a = %.3f  b = %.3f\n',x(1),x(2),x(3))
fprintf('  The sum of the squares is: %.3f\n\n',resnorm)

% ------------------------------------------------------------------
function f = ObjFunc(x,IDmodel)
global rA Pe Pea Ph2
switch IDmodel
case 'a'        % Model (a), x(1)=k; x(2)=Ke
    f = rA - x(1)*Pe.*Ph2./(1+x(2)*Pe);             
case 'b'        % Model (b), x(1)=k; x(2)=Ke; x(3)=Ka
    f = rA - (x(1)*Pe.*Ph2)./(1+x(3)*Pea+x(2)*Pe); 
case 'c'        % Model (c), x(1)=k; x(2)=Ke
    f = rA - x(1)*Pe.*Ph2./(1+x(2)*Pe).^2;
case 'd'        % Model (d), x(1)=k; x(2)=a; x(3)=b;
    f = rA - x(1)*Pe.^x(2).*Ph2.^x(3);
otherwise
    disp('Unknown model.')
end