function TempOpt_BatchReactor
% 對(duì)間歇反應(yīng)器中的連串-平行復(fù)雜反應(yīng)系統(tǒng)，計(jì)算使產(chǎn)物B收率最大的最優(yōu)反應(yīng)溫度
%
%   Author: Huang Huajiang
%   Copyright 2003 UNILAB Research Center, 
%   East China University of Science and Technology, Shanghai, PRC
%   $Revision: 1.0 $  $Date: 2003/08/11 $

clear all
clc

R = 8.31434;	                        % Gas constant, kJ/kmol K
k0 = [5.78052E+10  3.92317E+12  1.64254E+4  6.264E+8];  % Arrhenius constant, 1/s
Ea = [124670  150386  77954  111528];   % Activation energy, kJ/kmol

C0 = [1  0  0  0  0];   % 初始濃度C0(i), kmol/m^3
tspan = [0 1e4];
T0 = 200+273;           % Reactor temperature, Kelvin
T = fminsearch(@ObjFunc,T0,[],tspan,C0,k0,Ea,R);
fprintf('\t使產(chǎn)物B收率最大的最優(yōu)反應(yīng)溫度為%3.1f℃',T - 273)


% ------------------------------------------------------------------
function f = ObjFunc(T,tspan,C0,k0,Ea,R)
[t,C] = ode45(@MassEquations, tspan, C0,[],k0,Ea,R,T);
CBmax = max(C(:,2));
f = - CBmax;

% ------------------------------------------------------------------
function dCdt = MassEquations(t,C,k0,Ea,R,T)
% Reaction rate constants, 1/s
k = k0.*exp(-Ea/(R*T));
k(5) = 2.16667E-04;  

% Reaction rates, kmoles/m3 s
rA = -(k(1)+k(2))*C(1);
rB = k(1)*C(1)-k(3)*C(2);
rC = k(2)*C(1)-k(4)*C(3);
rD = k(3)*C(2)-k(5)*C(4);
rE = k(4)*C(3)+k(5)*C(4);

% Mass balances
dCdt = [rA; rB; rC; rD; rE];