% Copyright 1997 Thomas Mc Avoy
% MATLAB routine to set a plantwide control system for the TE process. User 
% enters tuning parameters for inner cascades, the level pairings, the
% PI parameters for the level controllers, and the outer loop pairings and
% PI tuning parameters.  The outer loops should be tested in closed loop.
% Variables to be plotted (up to 6) can be selected. The plantwide system
% is tested using IDV upsets, and several setpoint changes
%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 since recycle and agitator are not used
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(3600,6);
conlbl1=zeros(1,16);
conlbl=zeros(1,20);
con=[0 0];
kkk=[2 7 8];
KCP=0;
TRP=1.e6;
clf;
disp('********************************************************************')
disp('*                  PLANTWIDE CONTROL SYSTEM TEST FOR               *')
disp('*                      THE TENNESSEE EASTMAN PLANT                 *')
disp('*                    ****************************                  *')
disp('*                       COMPLETE PLANT TESTING                     *')
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
idv=zeros(20,1);% reset disturbaces here
	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 LOOP ARCHITECTURE & PI PARAMETERS')
	disp('    4.  INPUT OUTER LOOP ARCHITECTURE & PI PARAMETERS')
	disp('    5.  INPUT WHICH VARIABLE TO CHANGE')
    disp('    6.  SAVE YOUR PI PARAMETERS')
    disp('    7.  LOAD YOUR PI PARAMETERS')
	k=input(' ');
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(kkk,KCL,TRL); 
% kkk gives the LC pairings, KCL & TRL the PI parameters
elseif k==4
	[KCP,TRP,con,conlbl,conlbl1]=archit(kkk,KCP,TRP,con,conlbl,conlbl1);
elseif k==5
	[kk,d,titl1]=pltst(con,conlbl,conlbl1);
% kk variables to be changed, d bump size if applicable
	if isempty(d)
	else
	break
	end
elseif k==6
        % 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==7
% RETRIEVES PREVIOUSLY SAVED DATA    added by Yunbing 97-04-20
prevcase=input('Enter the name of the file to be retrieved as name.mat : ','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./3600;% time step = 1 second, (1/3600 hours)

load steady;% load in steady state values and fix set points
[xmeas,yp,y]=tefunc(0.,idv,u,y,1);% initialize simulation (iflag=1)
load steady;
% Calculate ratios
xmeas(42)=xmeas(40)/xmeas(41);
xmeas(43)=xmeas(23)/xmeas(25);
xmeas(44)=xmeas(26)/xmeas(27); 
setpt=[xmeas(ncas(:,2));50];
setptL=[xmeas(8);xmeas(12);xmeas(15)];% level measurements
[n,m]=size(con);
setptP=zeros(n,1);
setptP=xmeas(con(:,2));% set points for primary loops

% kk=variable to change
	if kk<=15
% kk<=15 means IDV(kk) to be tested
		idv(kk)=1;
	elseif kk==16;% production change of -15%
		disp(' ')
		disp('THE PRODUCTION RATE WILL BE RAMPED DOWN IN 1 HOUR')
		disp(' ')
		for i=1:n
			if con(i,1)==4;% C feed used for production control
				ramp=xmeas(17)*ones(1,3600)-[1:3600]*(9.6393e-4);
				setptP(i)=ramp(1);
				ik=i;
			end
		end
	elseif kk==17;% 50G/50H change to 40G/60H
		for i=1:n
			if con(i,1)==1;% D/E ratio feed used for product mix control
				setptP(i)=(40/60)*setptP(i);
			end
		end
	elseif kk==18;% reactor pressure change of -60kPa
		for i=1:n
			if con(i,2)==7;% find reactor pressure
				setptP(i)=setptP(i)-60;
			end
		end
	elseif kk==19;% 2% change in B composition in purge
		for i=1:n
			if con(i,2)==30;% find B composition
				setptP(i)=setptP(i)+2;
			end
		end
	else
		kP=kk-19;%find which loop to bump
		setptP(kP)=setptP(kP)*(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
%find G/H loop
for i=1:n
	if con(i,2)==42
		kD=i;
		setptDE=xmeas(2)/xmeas(3);
	end
end
% If E not used for reactor level control ratio it to C
kE=0;% don't ratio E-Feed to C-Feed if kE=0
for i=1:n
	if kkk(1)~=2;% check for E-Feed Reactor level pairing
		setptEE=xmeas(3)/xmeas(4);
		kE=1;
	end
end
time=0;
zz=zeros(12,1);% used in velocity form of PI controller
zzL=zeros(3,1);% used in velocity form of Level Controllers
zzP=zeros(n,1);% used in velocity form of Primary controllers
zzDE=0.;% used for D/E ratio
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 and outer controllers first to get set points for inner loops
	errL=setptL-[xmeas(8);xmeas(12);xmeas(15)];
	errP=setptP-xmeas(con(:,2));
	errDE=setptP(kD)-xmeas(con(kD,2));
	dxmvL=KCL'.*(errL-zzL+errL*deltat*60./TRL');
	dxmvP=KCP'.*(errP-zzP+errP*deltat*60./TRP');
	dxmvDE=KCP(kD)*(errDE-zzDE+errDE*deltat*60/TRP(kD));
	zzL=errL;
	zzP=errP;
	zzDE=errDE;
	if kk==16 & t>3600;% stop ramp after 3600 secs
		setptP(ik)=ramp(3600);
	elseif kk==16 & t<=3600
		setptP(ik)=ramp(t);
	end
	setpt(kkk)=setpt(kkk)+dxmvL;
		if kkk(1)==5
		dxrec=dxmvL(1);
		u(5)=u(5)+dxrec;
	end
% 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(1:4)=u(1:4)+dxmv(1:4);
	u(6:12)=u(6:12)+dxmv(6:12);
	setpt(con(:,1))=setpt(con(:,1))+dxmvP;
	setptDE=setptDE+dxmvDE;
	setpt(1)=setptDE*xmeas(3);% multiply D/E set point by E
	if kE==1;% ratio E-Feed to C
		setpt(2)=setptEE*xmeas(4);
	end
% now do inner cascades
	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;
% check for valve saturation	
while any(u>100);[ixx,jxx]=max(u);u(jxx)=100;end;
while any(u<0);[ixx,jxx]=min(u);u(jxx)=0;end;
	time=time+deltat;

%  if system shut down
%if xmeas(7)>=3000 | xmeas(8)>= 100 | xmeas(9)>=175 % | xmeas(12)>=12
%	yout(kt+1:tfin/60,:)=ones(tfin/60-kt-1,1)*xmeas(plt)';
%	tfin = t;
%	break
%end;

		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
orient tall
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)
% pause
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')
		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
