% ********  EXAMPLE 7.2  THREE-TANK MIXING PROCESS  ****
% ********             from  MARLIN               ******
% COPYRIGHT McMASTER UNIVERSITY, 1994
% WRITTEN BY : T. MARLIN
% FILENAME: EXAM7_2.M
% LAST UPDATED: November 1994
% VERSION: 2.0
% DESCRIPTION OF FILE:
%  this m-file calculates the transient response 
%  for closed-loop three-tank mixing system 
%  to a step change input
%  An exercise is to add the feedback control to this file
% *******************************************************

clear ; clf ; axis('auto') ; clg

% *******************************************************
%  PROCESS DATA
% *******************************************************

xain=100. ;         % 100% A in stream A
xabin=1   ;         %   1% A in stream B
fb=6.9  ;           % feed flow
v=35.0  ;           % volume of each tank
vp = 50;            % initialize valve at 50 % open
fa=.0028*vp ;       % initial flow of A

xa0=(fa*xain+fb*xabin)/(fa+fb) ;
xa1=xa0 ; xa2=xa0 ; xa3=xa0 ;     % initialize tanks

% disturbance data
dist1 = 70        ;  % disturbed value of valve position
tdist1 = 20.   ; 
%
% ***************************************************
% CONTROLLER  --  CONTROLLER IN THIS EXERCISE
% ***************************************************
%kc =  ;
%ti =  ;
%td =  ;

%xa3m1 = xa3 ; xa3m2 = xa3 ;  % past values for d-mode
% set point data

%  spinit= xa3  ;
%  spnew= 2.96      ;
%  tsp= 10      ;
%  sp = spinit;

% **************************************************
% SIMULATION PARAMETERS & INITIALIZATION
% **************************************************

tstart=0.    ;
tend=200.   ;
delt=0.1       ;  % delta t
%                          initialization
%sp = xa0  ;
%xa3m1=xa3;  xa3m2=xa3 ;
pierror=0.; pierrm1=0;

iae = 0.  ;
ise = 0.  ;
itae = 0. ;

cmax=round(tend/delt)  ;
kount = 100.; % initialize to store first sample
nkount = 0. ;
% *************************************************
% SIMULATION
% *************************************************

for cnt=1:1:cmax
    
if cnt*delt > tdist1         %  introduce disturbances
            xabin = 1.8 ;    %  

    end

%  execute the digital PID control algorithm
%   pierror = sp - xa3;
%   vp =  ;
%   
%   
%   

%  simulate the plant

    fa=.0028*vp  ;
    xa0 = (fa*xain+fb*xabin)/(fa+fb) ;

    xa1 = xa1 + delt* ((fa+fb)/v)*(xa0-xa1) ;
    xa2 = xa2 + delt* ((fa+fb)/v)*(xa1-xa2) ;
    xa3 = xa3 + delt* ((fa+fb)/v)*(xa2-xa3) ;

% *************************************************
%  SAVE AND PLOT RESULTS
% *************************************************

    if kount >=2   % skip store
          kount = 0;
          nkount = nkount+1 ;
          Y(1,nkount) = xa3   ;
          MV(1,nkount)=  vp    ;
          T(1,nkount) =  cnt*delt   ;
          DIST(1,nkount)  = xabin   ;
    end
    
    kount=kount+1 ;

%  **********************
%  CALCULATE PERFORMANCE MEASURES
% ***********************

   iae = iae + delt * abs(pierror);
   ise = ise + delt * (pierror*pierror);

end
subplot (2,2,1), plot (T(1,:),Y(1,:))
xlabel ('time')
ylabel ('tank 3 concentration of A')

%
subplot (2,2,3), plot (T(1,:),MV(1,:))
xlabel('time')
ylabel('valve % open')

subplot (2,2,4), plot (T(1,:),DIST(1,:))
xlabel ('time')
ylabel ('Stream B concentration of A')

figure (1)
pause

close all