% Fil dist5.m 
% Continuous column with variable holdup
% 
% states are:
%   First NT states: molde fraction light component on stage i
%   state NT+1: holdup condenser (MD)
%   state NT+2: holdup reboiler  (MB)
%   state NT+3: integrated error for PI-controller 1
%   state NT+4: integrated error for PI-controller 2
%

% In this version of the program we assume 
%   Bottom: MB is constant and V is constant
%   Top: MB controlled with PI-controller using reflux
%        xD controlled with PI-controller using D
%
% -----------------------------------
% Example of Use of program:
%
% Initial states (xb=0, x=0.04 on the other stages, xd=0.9, MD=0.1, MB=1
% Xinit = [0 0.04*ones(1,34) 0.9 0.001 1 0 0];
% [T,X]=ode15s('dist5',[0 5],Xinit);  % run 5 hours
% len=size(X);[T,X]=ode15s('dist5',[0 2],X(len(1),:));% run another 2h
% 
% Empty condenser and run another 2 hours; reset integrators
% len=size(X); Xinit2= X(len(1),:); 
% NT=36; Xinit2(NT+1)=0.001; Xinit2(NT+3)=0; Xinit2(NT+4)=0;
% [T,X]=ode15s('dist5',[0 2],Xinit2)
%
% plot(T,X(:,37)); %plot condenser level
% plot(T,X(:,36)); %plot condenser composition, 
% plot(T,X(:,1));  %plot reboiler composition, 
% etc..
%
% Comment: The tunings below do result in poor operation!
%
% To plot flows (after integration is finished)
% NT=36;
% 
% % Plotte LT
% MD=X(:,NT+1);  MDs = 0.08;
% err = MD - MDs; interr = X(:,NT+3);
% VB=1; KC=-5;taui=0.05;L0I=0.8*VB;
% LT = L0I - KC*err - KC/taui*interr; % PI-regulator
% plot(T,LT)
% 
% % Plotte D
% xNT = X(:,NT); XDs=0.95;
% errxd =  xNT-XDs; interrxd = X(:,NT+4);
% Kcd=100; tauid=1; D0=0;
% D = D0 + Kcd*errxd + Kcd/tauid*interrxd; % PI-regulator
% plot(T,LT)
% 
% END: Example of Use of program:
% -----------------------------------

function xdot = dist3(t,x)

% print time
t

% Constant relative volatility assumed. Binary mixture
alfa = 3;

% Total no. of stages including total condenser
NT=36;

% Feed stage (numbered from bottom)
NF=31;

F=1;  % mol/h  (may set = 0 to get closed operation)
qF= 1.24; % fraction liquid
zF=0.04;  % feed composition
VB=1;
%L0=1; 
D=0;
%B=0;
MD=x(NT+1); % holdup in condenser = state NT+1 
MB=x(NT+2);
MK=0.0001;
%KC=0;
%XFs=0.5;

% PI controller for level using reflux
KC=-5;
taui=0.05;
MDs = 0.08;
err = MD - MDs; 
L0I=0.8*VB; 
interr = x(NT+3);
LT = L0I - KC*err - KC/taui*interr; % PI-regulator
if LT < 0
  LT=0.001;
  err=0; % turn off integration to avoid windup
end

% PI-controller for XD using D
XDs=0.95;
errxd =  x(NT)-XDs; 
interrxd = x(NT+4);
Kcd=100; tauid=1; D0=0;
D = D0 + Kcd*errxd + Kcd/tauid*interrxd; % PI-regulator
if D < 0
  D=0.0001;
  errxd=0; % turn off integration to avoid windup
end


% Damp-vaeske likevekt
i=1:NT-1;    
y(i)=alfa*x(i)./(1+(alfa-1)*x(i));

% Dampstroem 
i=1:NT-1;    
V(i)=VB*ones(1,NT-1);
i=NF:NT-1;
V(i)=V(i)+(1-qF)*F;

% Vaeskestroem
i=2:NT;
L(i)=LT*ones(1,NT-1);
i=2:NF;
L(i)=L(i)+qF*F;

% Total kondensator
y(NT)=x(NT);

% Holdup i kolonnen (konstant)
i=2:NT-1;
M(i)=MK*ones(1,NT-2);

% Holdup i kokeren
% Anita: Antar konstant holdup
B = L(2) - V(1);
MBdot=L(2)-V(1)-B;

% Holdup i kondensatoren
MDdot=V(NT-1)-L(NT)-D;

% Dynamiske ligninger (komponentbalanser)

% Kokeren
Mxdot(1) = L(2)*x(2)-V(1)*y(1)-B*x(1);
xdot(1)=(Mxdot(1)-x(1)*MBdot)/MB;

% Kolonnen
i=2:NF-1;
xdot(i)=(L(i+1).*x(i+1)' - L(i).*x(i)' + V(i-1).*y(i-1) - V(i).*y(i))./M(i);
i=NF+1:NT-1;
xdot(i)=(L(i+1).*x(i+1)' - L(i).*x(i)' + V(i-1).*y(i-1) - V(i).*y(i))./M(i);

% Foedetrinn
xdot(NF) =(L(NF+1)*x(NF+1)+V(NF-1)*y(NF-1)-L(NF)*x(NF)-V(NF)*y(NF)+F*zF)/MK;

% Kondensatoren
Mxdot(NT) = V(NT-1)*y(NT-1)-L(NT)*x(NT)-D*x(NT);
xdot(NT)=(Mxdot(NT)-x(NT)*MDdot)/MD;

% Funksjonen f(t,x) maa returnere en kolonnevektor
xdot(NT+1)=MDdot;
xdot(NT+2)=MBdot;

%Integrator in controller
xdot(NT+3)=err;
xdot(NT+4)=errxd;

xdot=xdot(:);






