function BatchDistill
% 恒定滯液量的雙組分間歇精餾的回流比控制
%
% 假設：
% 1. 沸騰上升的蒸汽流量（V）恒定；
% 2. 下降的液體摩爾流量（L）恒定； 
% 3. 塔板上的滯液量（M）恒定；
% 4. 再沸器和塔頂冷凝器中的滯液量恒定。
%
%   Author: HUANG Huajiang
%   Copyright 2003 UNILAB Research Center, 
%   East China University of Science and Technology, Shanghai, PRC
%   $Revision: 1.0 $  $Date: 2003/01/24 $

clear all
global alpha M0 M V Rinit x0set KC Flag tc

alpha = 3.5;      % 相對揮發度 
M0 = 100;         % 塔頂冷凝器中的摩爾滯液量, kmol 
M =  5;           % 塔板上的摩爾滯液量, kmol
V = 10;           % 沸騰上升的蒸汽流量, kmol/h
Rinit = 1E10;     % 全回流條件（全回流-剛開始時）
charge = 2500;    % 蒸餾塔負荷（處理量）, kmol
xB = 0.8;         % 蒸餾塔釜進料組成, 摩爾分率

% 全回流時, Flag = 0
Flag = 0;

% 控制參數
x0set = 0.9;              % 餾出液組成的設定值
KC = 100;                 % 控制增益

dt = 0.01;
tf = 300;

% 初始值
MB = charge;
dist = 0;                       % 總的餾出液量
x0 = [zeros(1,8) xB MB dist];   % x0 x1 x2 x3 x4 x5 x6 x7 xB MB dist
tspan = [0 tf];
[t,x] = ode45(@DistilEqs,tspan,x0)
plot(t,x(:,1),'r-',t,x(:,2),'k-.',t,x(:,3),'b--',t,x(:,4),'k-',t,x(:,5),'b-.',t,x(:,6),'r--',t,x(:,7),'b-',t,x(:,8),'k--')
xlabel('時間 (h)')
ylabel('液相摩爾分率')
legend('x_0','x_1','x_2','x_3','x_4','x_5','x_6','x_7')
x0 = x(:,1);
Flag = ones(size(t,1),1);
Flag(find(t<tc)) = 0;    % t>=tc時，開始進入控制狀態（flag = 1）
Rc = KC*(x0set-x0).*Flag;
figure
plot(t,Rc)
xlabel('time (h)')
ylabel('Rc')

function dxdt = DistilEqs(t,x)
global alpha M0 M V Rinit x0set KC Flag tc

% x(1...11)表示[x0 x1 x2 x3 x4 x5 x6 x7 xB MB dist]
xB = x(9);
MB = x(10);
dist = x(11);

% 當塔頂餾分的組成超過設定值x0set時，進入控制狀態，記下此時的時刻tc
if x(1)>=x0set
    Flag = 1;
    tc = t;
end

%  控制方程
Rc = KC*(x0set-x(1))*Flag;
if Rc<0    % 限制Rc> = 0.0
    Rc = 0;
end

R = Rc+Rinit*(1-Flag);
D = V/(R+1);
L = V-D;

% 汽液平衡關系
i=2:9
y(i) = alpha*x(i)./(1+(alpha-1)*x(i));

% 塔頂冷凝器
dxdt(1) = (V*y(2)-(L+D)*x(1))/M0;      % d/dt(x0)

% 第1~7塊塔板的易揮發組分質量平衡方程
for i=2:8
    dxdt(i) = (L*(x(i-1)-x(i))+V*(y(i+1)-y(i)))/M; 
end

% 塔釜
dxdt(9) = (L*x(8)-V*y(9))/MB;       % d/dt(xB)   
dxdt(10) = L - V;                   % d/dt(MB)   

% 餾出液總量
dxdt(11) = D;                       % d/dt(dist)

dxdt = dxdt';