% Copyright 1997 Thomas Mc Avoy
% Routine to tune the inner cascade PI controllers for the Tennessee Eastman 
% process.  A reaction curve is generated for the selected controller, from
% which PI tuning constants can be determined. This routine calls the TE
% 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.  
% k gives the valve to be perturbed for bump test

%global yp y idv valve xmeas u index

clf
idv=zeros(20,1);
load steady;
clc;
disp('********************************************************************')
disp('*               REACTION CURVE TESTS FOR CASCADE LOOPS             *')
disp('*                      THE TENNESSEE EASTMAN PLANT                 *')
disp('*                    ****************************                  *')
disp('*                ALL CASCADES EXCEPT REACTOR CW TEMP               *')
disp('*              ========================================            *')
disp('*                                                                  *')
disp('********************************************************************')
disp(' ')
disp(' ')
while 1
while 1
	disp(' ')
	disp('SELECT FROM THE FOLLOWING MENU')
	disp('    1.  INPUT PI PARAMETERS FOR REACTOR CW TEMPERTURE')
	disp('    2.  CARRY OUT BUMP TEST FOR INNER CASCADE')
	k=input(' ');
	if k==2
		break
	elseif k==1
		disp(' ')
		disp('INPUT THE CONTROLLER GAIN FOR REACTOR CW TEMP');
		temp=input(' ');
		while (isempty(temp)==1)
			disp('PLEASE ENTER KC AGAIN')
			temp=input(' ');
		end
		KC=temp;
		disp(' ')
		disp('INPUT THE RESET TIME (min) FOR REACTOR CW TEMP');
		temp=input(' ');
		while (isempty(temp)==1)
			disp('PLEASE ENTER TR (min) AGAIN')
			temp=input(' ');
		end
		TR=temp;
	else
	end
end

[d,k,titl,titl1]=whicas;

% numerical integration begins here
deltat=1.000001/3600;% time step = 1 second, (1/3600 hours)
% A value slightly greater than 1 second is needed to synchronize with
% the FORTRAN routines for sampled variables


ncas=[1,2;2,3;3,1;4,4;5,5;6,10;7,14;8,17;9,19;10,21;11,22];
% ncas gives cascade measurement and its associated manipulated valve.

load steady;% load in steady state values
[xmeas,yp,y]=tefunc(0.,idv,u,y,1);% initialize simulation (iflag=1)
load steady; 
rcwt=xmeas(21);
uout=u(k)*d;
you0=xmeas(ncas(k,2));
% bump valve k
u(k)=u(k)*(1+d);

% main integration begins here. Use Euler method and integrate for tfin sec
disp(' ')
	tfin=input ('INPUT FINAL SIMULATION TIME (sec): ');
	while (isempty(tfin)==1)
		disp('PLEASE ENTER FINAL TIME AGAIN. ')
		d=input ('INPUT FINAL SIMULATION TIME (sec): ');
	end	
disp(' ')
% delta t is 1 second
time=0;
zz=0;% used in velocity form of PI controller
for t=1:tfin+1
	[xmeas,yp,y]=tefunc(time,idv,u,y,0);
	err=rcwt-xmeas(ncas(10,2));
	dxmv=KC*(err-zz+err*deltat*60/TR);% deltat is in hours, TR in minutes
	zz=err;
	u(10)=u(10)+dxmv;
	u10(t)=u(10);
	zout(t)=u(10);
	time=time+deltat;
	yout(t)=xmeas(ncas(k,2));
	z1out(t)=xmeas(21);
	y=y+yp*deltat;
end


plot([0:tfin],[yout(1:tfin+1)])
title(titl1)
xlabel('Time sec')
ylabel(titl)

disp('HIT RETURN TO CONTINUE')
pause

	disp(' ')
	disp('WOULD YOU LIKE TO TEST ANOTHER LOOP?')
	an=input ('INPUT 1 FOR YES AND 0 TO EXIT ');
	while (isempty(an)==1)&an~=0&an~=1
		disp('PLEASE ENTER A 1 OR 0')
		an=input ('WOULD YOU LIKE TO TEST ANOTHER LOOP?');
	end	
	disp(' ')

	if an==0
		break
	end
end


