% Copyright 1997 Thomas Mc Avoy
% Routine to tune the reactor cooling water PI controller for the Tennessee 
% Eastman process.  Since the process is open loop unstable, this loop is 
% tuned closed loop. This m-file calls the Tennessee Eastman
% simulation, TEFUNC.  TEFUNC is a mex file that allows Fortran to be called.
% returns the measurement vector xmeas(41), the state vector y(50), 
% and the derivative of y, yp(50).  It is called as follows
% [xmeas,yp,y]=tefunc(time,IDV,y,iflag).  IDV is the disturbance
% vector. iflag=1 initializes all of the constants in the common blocks
% as well as the state vector and its derivatives.  
% A closed loop step test is carried out on the reactor cooling water temp.
% Inputs which are requested are the PI parameters and final simulation time.

function reactorcw
%global yp y idv valve xmeas u index
idv=zeros(20,1);
load steady;% steady state for TE process
clc;
disp('********************************************************************')
disp('*               REACTION CURVE TESTS FOR CASCADE LOOPS             *')
disp('*                      THE TENNESSEE EASTMAN PLANT                 *')
disp('*                    ****************************                  *')
disp('*                       REACTOR COOLING WATER                      *')
disp('*              ========================================            *')
disp('*                                                                  *')
disp('********************************************************************')
disp(' ')


% Input the PI parameters and the final time
disp(' ')
	KC=input ('INPUT THE CONTROLLER GAIN FOR THE REACTOR CW TEMP: ');
	while (isempty(KC)==1)
		disp('PLEASE ENTER KC AGAIN. ')
		KC=input ('INPUT THE CONTROLLER GAIN FOR THE REACTOR CW TEMP: ');
	end
disp(' ')
	TR=input ('INPUT THE RESET TIME (min) FOR THE REACTOR CW TEMP: ');
	while (isempty(TR)==1)
		disp('PLEASE ENTER TR AGAIN. ')
		TR=input ('INPUT THE RESET TIME FOR THE REACTOR CW TEMP: ');
	end	
disp(' ')
	tfin=input ('INPUT FINAL SIMULATION TIME (sec) FOR THE REACTOR CW TEMP: ');
	while (isempty(tfin)==1)
		disp('PLEASE ENTER FINAL TIME AGAIN. ')
		tfin=input ('INPUT FINAL SIMULATION TIME (sec) FOR THE REACTOR CW TEMP: ');
	end	
disp(' ')
disp(' ')	
disp('A +1 deg C STEP IN THE REACTOR CW TEMPERATURE WILL BE MADE')
disp(' ')

% numerical integration begins here
deltat=1/3600;% time step = 1 second, (1/3600 hours)

[xmeas,yp,y]=tefunc(0.,idv,u,y,1);% initialize simulation (iflag=1)
load steady;% load in steady state values 
rcwt=xmeas(21);	

% main integration begins here. Use Euler method and integrate for tfin secs
% delta t is 1 second
time=0;
zz=0;
rcwt=rcwt+1;% 1 deg step in reactor CW temp
tfin=tfin+1;
for t=1:tfin+1
	[xmeas,yp,y]=tefunc(time,idv,u,y,0);
	err=rcwt-xmeas(21);
	dxmv=KC*(err-zz+err*deltat*60/TR);% deltat is in hours
	zz=err;
	u(10)=u(10)+dxmv;
	zout(t)=u(10);
	time=time+deltat;
	yout(t)=xmeas(21);
	z1out(t)=xmeas(8);% Reactor level
	z2out(t)=xmeas(7);% Reactor pressure
	z3out(t)=xmeas(6);% Reactor feed
	y=y+yp*deltat;
end
subplot(2,2,1),plot([0:tfin],yout)
title('Reactor CW Temperature')
xlabel('Time sec')
subplot(2,2,2),plot([0:tfin],z1out)
title('Reactor Level')
xlabel('Time sec')
subplot(2,2,3),plot([0:tfin],z2out)
title('Reactor Pressure')
xlabel('Time sec')
subplot(2,2,4),plot([0:tfin],z3out)
title('Reactor Feed')
xlabel('Time sec')



