function CSTRs_Closedloop
% 三個串聯(lián)的CSTR等溫反應器（閉環(huán)系統(tǒng)--比例反饋控制器）
%
%   Author: HUANG Huajiang
%   Copyright 2002 UNILAB Research Center, 
%   East China University of Science and Technology, Shanghai, PRC
%   $Revision: 1.0 $  $Date: 2002/07/20 $

clear all
clc

global k Kc CA0 CAD tau CA3set taui

% initial conditions
CA = 1.8;           % kmol/m3
CA1 = 0.4;          % kmol/m3
CA2 = 0.2;          % kmol/m3
CA3 = 0.1;          % kmol/m3
CA3set = 0.1;
EI = 0;

% disturbance
CAD = 0.2;

% parameters
k = 0.5;            % min-1
tau = 2;
Kc = 30;
taui = 5;
stoptime = 8.5;     % min

% Solve ODEs
[t,y] = ode45(@Equations,[0:0.5:stoptime],[CA1 CA2 CA3 EI]);

% Output results
CA1 = y(:,1);
CA2 = y(:,2);
CA3 = y(:,3);
EI = y(:,4);
E = CA3set - CA3;
CAM = 0.8+Kc*(E+EI/taui);
disp(' Results:')
disp('       t        CA1       CA2       CA3      CAM')
disp([t y(:,1:3) CAM])
plot(t,CA1,'b--',t,CA2,'m-.',t,CA3,'r-',t, CAM,'k:')
axis([0 t(end)+0.5 min(CA3)-0.01 max(CAM)])
xlabel('Time (min)')
ylabel('C_A_1, C_A_2, C_A_3, C_A_M (kmol/m^3)')
legend('C_A_1','C_A_2','C_A_3','C_A_M')

% ------------------------------------------------------------------
function dydt = Equations(t,y)
global k Kc CA0 CAD tau CA3set taui
CA1 = y(1);
CA2 = y(2);
CA3 = y(3);
EI = y(4);

% feedback controller
E = CA3set - CA3;
CAM = 0.8+Kc*(E+EI/taui);
CA0 = CAM+CAD;

% Model equations
dCA1dt = (CA0-CA1)/tau - k*CA1;
dCA2dt = (CA1-CA2)/tau - k*CA2;
dCA3dt = (CA2-CA3)/tau - k*CA3;
dEIdt = E;
dydt = [dCA1dt; dCA2dt; dCA3dt;dEIdt];