% Copyright 1997 Thomas Mc Avoy
% MATLAB routine to check tuning of primary loops in TE process.  User 
% enters tuning parameters for inner cascades which should then be tested
% in closed loop.  Next the level loops can be tuned and tested in closed loop.
% variables to be plotted (up to 6) can be selected.
%clear
%global yp y idv valve xmeas u index
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.
idv=zeros(20,1);
load steady;
KC=zeros(12,1);% cascade KC (5,12) dummy variables
TR=ones(12,1)*1.e6;% cascade reset times (5,12) dummy variables
TRL=ones(1,3)*1.e6;% level reset time
KCL=zeros(1,3);% level KC
yout=zeros(1080,6);
clf;
disp('********************************************************************')
disp('*               REACTION CURVE TESTS FOR CASCADE LOOPS             *')
disp('*                      THE TENNESSEE EASTMAN PLANT                 *')
disp('*                    ****************************                  *')
disp('*                     PRIMARY LOOP BUMP TESTS                      *')
disp('*              ========================================            *')
disp('*                                                                  *')
disp('********************************************************************')
disp(' ')
kkk=[2 7 8];% initialize with E-react level, Sep exit-sep level, prod-strip level
[plt,titl,pltno]=varplo1;% initialize plotting
while 1
while 1
	disp(' ')
	disp('SELECT FROM THE FOLLOWING MENU')
	disp('    1.  INPUT PI PARAMETERS FOR INNER CASCADE LOOPS')
	disp('    2.  INPUT VARIABLES TO BE PLOTTED')
	disp('    3.  INPUT LEVEL CONTROLLER PARAMETERS')
	disp('    4.  INPUT WHICH VARIABLE TO BUMP')
    disp('    5.  SAVE YOUR PI PARAMETERS')
    disp('    6.  LOAD YOUR PI PARAMETERS')
%	disp('    0.  Return to Matlab prompt')
	k=input(' ');
%if k==0
%	return;
%else
if k==1
	[KC,TR]=cascpi(KC,TR);
% KC & TR give the PI parameters for inner cascades
elseif k==2
	[plt,titl,pltno]=varplo1;
% plt gives variables to be plotted, titl plot titles, and pltno # plots
elseif k==3
	[kkk,KCL,TRL]=levcon([2 7 8],KCL,TRL); 
% kkk gives the LC pairings, KCL & TRL the PI parameters
elseif k==4
	[kk,d,titl1]=pribmp;
% kk variables to be bumped, d bump size
	break
elseif k==5
	% SAVE DATA TO FILE  added by Yunbing 97-04-20
	savecase=input('Please enter the file name : ','s'); 
	if isempty(savecase)~=1
	    disp('                          Saving Data To A File')
	    eval(['save ',savecase])            % SAVES DATA
	    pause(1);
	end
elseif k==6
% RETRIEVES PREVIOUSLY SAVED DATA    added by Yunbing 97-04-20
prevcase=input('Enter the name of the file to be retrieved : ','s');
 if (exist(prevcase)>=2) & (isempty(prevcase)~=1)
	disp('                      Loading Data From File');
	eval(['load ',prevcase])        % LOADS IN FILE
	pause(1);
 elseif (isempty(prevcase)~=1)
    disp('The file cannot be found, please try again.')
    pause(2);
 end
end
end

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

load steady;% load in steady state values and fix set points
setpt=[xmeas(ncas(:,2));50];
setptL=[xmeas(8);xmeas(12);xmeas(15)];% level measurements
[xmeas,yp,y]=tefunc(0.,idv,u,y,1);% initialize simulation (iflag=1)
load steady; 
DEflg=0;
% bump variable kk
	if kk<13&kk~=5&kk~=12
% kk<13 means cascade, recycle, agitator to be bumped
		setpt(kk)=setpt(kk)*(1+d);
	elseif kk==5|kk==12
% kk=5 recycle, kk=12 agitator
		u(kk)=u(kk)*(1+d);
	elseif kk==16
% kk=16 means D/E to be bumped
		DEflg=1;
		setptD=(setpt(1)/setpt(2))*(1+d);
	elseif kk>12&kk~=16
% kk>12&kk~=16 means level set point to be bumped
		setptL(kk-12)=setptL(kk-12)*(1+d);
	end
	
% main integration begins here. Use Euler method and integrate for tfin sec
disp(' ')
	tfin=input ('INPUT FINAL SIMULATION TIME (min): ');
	while (isempty(tfin)==1)
		disp('PLEASE ENTER FINAL TIME AGAIN. ')
		d=input ('INPUT FINAL SIMULATION TIME (min): ');
	end	
disp(' ')
tfin=tfin*60;% routines use 1 second integration time. convert to secs

time=0;
zz=zeros(12,1);% used in velocity form of PI controller
zzL=zeros(3,1);% used in velocity form of Level Controllers
minit=60;% plot first point when minutes are used. For tfin < 3hours 
% plots are in seconds
kt=1;% used for plotting in minutes

for t=1:tfin+1
	[xmeas,yp,y]=tefunc(time,idv,u,y,0);
	xmeas(42)=xmeas(40)/xmeas(41);
	xmeas(43)=xmeas(23)/xmeas(25);
	xmeas(44)=xmeas(26)/xmeas(27);
% do level controllers first to get set points for inner loops
	errL=setptL-[xmeas(8);xmeas(12);xmeas(15)];
	dxmvL=KCL'.*(errL-zzL+errL*deltat*60./TRL');
	zzL=errL;
	setpt(kkk)=setpt(kkk)+dxmvL;
% now do inner cascades
	if DEflg==1
		setpt(1)=setptD*xmeas(3);
	end;
	err=setpt-[xmeas(ncas(:,2));50];
	dxmv=KC.*(err-zz+err*deltat*60./TR);% deltat is in hours, TR in minutes
	zz=err;
	u=u+dxmv;
	time=time+deltat;
		if tfin>3600&minit==60
			yout(kt,:)=xmeas(plt)';
			kt=kt+1;
			minit=0;
		elseif tfin<=3600
			yout(t,:)=xmeas(plt)';
		else
		end
	y=y+yp*deltat;% integrate state equations
	minit=minit+1;
end
% Integration finished

xlbl='Time sec';
if tfin>3600
	tfin=tfin/60;
	xlbl='Time min';
end

subplot(3,2,1),plot([0:tfin-1],yout(1:tfin,1))
ylabel(titl(1,:))
title(titl1)
subplot(3,2,2),plot([0:tfin-1],yout(1:tfin,2))
ylabel(titl(2,:))

subplot(3,2,3),plot([0:tfin-1],yout(1:tfin,3))
ylabel(titl(3,:))

subplot(3,2,4),plot([0:tfin-1],yout(1:tfin,4))
ylabel(titl(4,:))

subplot(3,2,5),plot([0:tfin-1],yout(1:tfin,5))
ylabel(titl(5,:))
xlabel(xlbl)
subplot(3,2,6),plot([0:tfin-1],yout(1:tfin,6))
ylabel(titl(6,:))
xlabel(xlbl)
loop=0;
while loop==0
disp(' ')
disp('WOULD YOU LIKE TO PLOT A SINGLE VARIABLE?')
disp('INPUT 1 FOR A SINGLE PLOT OR ANYTHING ELSE TO CONTINUE ON')
disp(' ')
	temp=input(' ');
	if temp==1
		disp('THE FOLLOWING VARIABLES CAN BE PLOTTED. PLEASE SELECT ONE.')
		disp(' ')
		plos=['	1	';'	2	';'	3	';'	4	';'	5	';'	6	'];
		plos=[plos,titl];
		disp(plos)
		disp(' ')
		while 1
			tem=input(' ');
			while (isempty(tem)==1)|tem<1|tem>6
				disp('PLEASE INPUT A VALUE BETWEEN 1 AND 6')
				tem=input(' ')
			end;
			clf
			plot([0:tfin-1],yout(1:tfin,tem))
			ylabel(titl(tem,:))
			xlabel(xlbl)
			title(titl1)
			break
		end
	else
		loop=1;
		disp(' ')
		disp('HIT RETURN TO CONTINUE')
		pause
	end
end
	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

