tesfunc.m

Tennessee eastman 的Matalb仿真程序

    ADISTVAL = S5;
    BDISTVAL = SP5;
    CDISTVAL = (3 * (S1 - S5) - H5 * (S1P + 2 * SP5)) / H5^2;
    DDISTVAL = (2 * (S5 - S1) + H5 * (S1P + SP5)) / H5^3;
    TNEXTVAL = TLASTVAL + H5;
	ADIST(I,1) = ADISTVAL; 
	BDIST(I,1) = BDISTVAL;
	CDIST(I,1) = CDISTVAL;
	DDIST(I,1) = DDISTVAL;
	TNEXT(I,1) = TNEXTVAL;  
	%
	% End of TESUB5
	%
	end
end
for I = 10:12
    if TIME >= TNEXT(I,1)
    HWLK = TNEXT(I,1) - TLAST(I,1);
    SWLK = ADIST(I,1) + HWLK * (BDIST(I,1) + HWLK * (CDIST(I,1) + HWLK * DDIST(I,1)));
    SPWLK = BDIST(I,1) + HWLK * (2 * CDIST(I,1) + 3 * HWLK * DDIST(I,1));
    TLAST(I,1) = TNEXT(I,1);
    	if SWLK > 0.1
      	ADIST(I,1) = SWLK;
      	BDIST(I,1) = SPWLK;
      	CDIST(I,1) = -(3 * SWLK + 0.2 * SPWLK) / 0.01;
      	DDIST(I,1) = (2 * SWLK + 0.1 * SPWLK) / 0.001;
      	TNEXT(I,1) = TLAST(I,1) + 0.1;
      	else
      	ISD = -1;
      	HWLK = HSPAN(I,1) * TESUB7(ISD) + HZERO(I,1);
     	ADIST(I,1) = 0;
      	BDIST(I,1) = 0;
      	CDIST(I,1) = (IDVWLK(I,1)) / HWLK^2;
      	DDIST(I,1) = 0;
      	TNEXT(I,1) = TLAST(I,1) + HWLK;
  		end
	end
end
if TIME == 0
	for I = 1:12
    ADIST(I,1)=SZERO(I,1);
    BDIST(I,1)=0;
    CDIST(I,1)=0;
    DDIST(I,1)=0;
    TLAST(I,1)=0.0;
    TNEXT(I,1)=0.1;
  	end
end

XST(1,4) = TESUB8(1,TIME) - IDV(1,1) * 0.03 - IDV(2,1) * 2.43719 * 10 ^ (-3);
XST(2,4) = TESUB8(2,TIME) + IDV(2,1) * 0.005;
XST(3,4) = 1 - XST(1,4) - XST(2,4);
TST(1,1) = TESUB8(3,TIME) + IDV(3,1) * 5;
TST(4,1) = TESUB8(4,TIME);
TCWR = TESUB8(5,TIME) + IDV(4,1) * 5;
TCWS = TESUB8(6,TIME) + IDV(5,1) * 5;
R1F = TESUB8(7,TIME);
R2F = TESUB8(8,TIME);
for I = 1:3
      UCVR(I,1)=YY(I,1);
      UCVS(I,1)=YY(I+9,1);
      UCLR(I,1)=0.0;
      UCLS(I,1)=0.0;
end
for I = 4:8
      UCLR(I,1) = YY(I,1);
      UCLS(I,1) = YY(I+9,1);
end
  
for I = 1:8
      UCLC(I,1) = YY(I+18,1);
      UCVV(I,1) = YY(I+27,1);
end
  
      ETR = YY(9,1);
      ETS = YY(18,1);
      ETC = YY(27,1);
      ETV = YY(36,1);
      TWR = YY(37,1);
      TWS = YY(38,1);
	  
%	  TWS shows an error here.  For the initialization,
%     it's value is  7.729698353000001e+01, when it should be
%     77.296983530000
	  
for I = 1:12
      VPOS(I,1) = YY(I+38,1);
end
      UTLR = 0.0;
      UTLS = 0.0;
      UTLC = 0.0;
      UTVV = 0.0;
	  
for I = 1:8
      UTLR = UTLR + UCLR(I,1);
      UTLS = UTLS + UCLS(I,1);
      UTLC = UTLC + UCLC(I,1);
      UTVV = UTVV + UCVV(I,1);
end
for I = 1:8
      XLR(I,1) = UCLR(I,1)/UTLR;
      XLS(I,1) = UCLS(I,1)/UTLS;
      XLC(I,1) = UCLC(I,1)/UTLC;
      XVV(I,1) = UCVV(I,1)/UTVV;
end
  
      ESR = ETR/UTLR;
      ESS = ETS/UTLS;
      ESC = ETC/UTLC;
      ESV = ETV/UTVV;
% Additional precision errors show up in XLR XLS XLC XVV ESR ESS ESC ESV 

%     TESUB2(XLR,TCR,ESR,0);
%	  TESUB2(Z,T,H,ITY)
% Substitution of TESUB2 for call statement. MWB
% TESUB2 appears to update T and only T.
Z = XLR;
T = TCR;
H = ESR;
ITY = 0;
	  TIN=T;
for J = 1:100
	%      TESUB1(Z,T,HTEST,ITY);
	%	   TESUB1VAL = TESUB1(Z,T,H,ITY)	
	% Substitution of TESUB1 for call statement MWB
	%
	if ITY == 0
    	HTEST = 0;
		for I = 1:8
      		HI = T * (AH(I,1) + BH(I,1) * T / 2 + CH(I,1) * T^2 / 3);
      		HI = 1.8 * HI;
      		HTEST = HTEST + Z(I,1) * XMW(I,1) * HI;
		end
	else
      	HTEST = 0;
		for I = 1:8 
      		HI = T * (AG(I,1) + BG(I,1) * T / 2 + CG(I,1) * T^2 / 3);
      		HI = 1.8 * HI;
      		HI = HI + AV(I,1);
      		HTEST = HTEST + Z(I,1) * XMW(I,1) * HI;
		end
	end
	if ITY == 2
      	R1 = 3.57696 / 1000000;
     	HTEST = HTEST - R1 * (T + 273.15);
  	end
	%
	% End of TESUB1 sub
	%
	
    ERR = HTEST - H;
	
	% 		TESUB3VAL = TESUB3(Z,T,DH,ITY)
	%		TESUB3VAL = TESUB3(Z,T,DH,ITY)
	% Substitution of TESUB3 for call statement MWB
	%
	if ITY == 0
      	DH = 0.0;
    	for I = 1:8
      		DHI = AH(I,1) + BH(I,1) * T + CH(I,1) * T^2;
      		DHI = 1.8 * DHI;
      		DH = DH + Z(I,1) * XMW(I,1) * DHI;
		end
	else
      	DH=0.0;
    	for I = 1:8
      		DHI = AG(I,1) + BG(I,1) * T + CG(I,1) * T^2;
      		DHI = 1.8 * DHI;
      		DH = DH + Z(I,1) * XMW(I,1) * DHI;
		end
	end
	if ITY == 2
      	R3 = 3.57696/1000000;
     	DH = DH - R3;
 	end
	%
	% End of TESUB3 substitution.
	%
    DT = -ERR / DH;
    T = T + DT;
  	if abs(DT) < 1*10^(-12)
    	break
  	end
  	if J == 100
  		T = TIN;
  	end
end
TCR = T;
% End of TESUB2.

      TKR = TCR+273.15;

%     TESUB2(XLS,TCS,ESS,0);
%	  TESUB2(Z,T,H,ITY)
% Substitution of TESUB2 for call statement. MWB
% TESUB2 appears to update T and only T.
Z = XLS;
T = TCS;
H = ESS;
ITY = 0;
TIN=T;
for J = 1:100
	%      TESUB1(Z,T,HTEST,ITY);
	%	   TESUB1VAL = TESUB1(Z,T,H,ITY)	
	% Substitution of TESUB1 for call statement MWB
	%	
	if ITY == 0
		HTEST = 0;
		for I = 1:8
			HI = T * (AH(I,1) + BH(I,1) * T / 2 + CH(I,1) * T^2 / 3);
			HI = 1.8 * HI;
			HTEST = HTEST + Z(I,1) * XMW(I,1) * HI;
		end
	else
		HTEST = 0;
		for I = 1:8 
			HI = T * (AG(I,1) + BG(I,1) * T / 2 + CG(I,1) * T^2 / 3);
			HI = 1.8 * HI;
			HI = HI + AV(I,1);
			HTEST = HTEST + Z(I,1) * XMW(I,1) * HI;
		end
	end
	if ITY == 2
		R1 = 3.57696 / 1000000;
		HTEST = HTEST - R1 * (T + 273.15);
	end
	%
	% End of TESUB1 sub
	%
 	ERR = HTEST - H;
	% 		TESUB3VAL = TESUB3(Z,T,DH,ITY)
	%		TESUB3VAL = TESUB3(Z,T,DH,ITY)
	% Substitution of TESUB3 for call statement MWB
	%
	if ITY == 0
   		DH = 0.0;
   		for I = 1:8
   			DHI = AH(I,1) + BH(I,1) * T + CH(I,1) * T^2;
   			DHI = 1.8 * DHI;
   			DH = DH + Z(I,1) * XMW(I,1) * DHI;
		end
	else
   		DH=0.0;
   		for I = 1:8
   			DHI = AG(I,1) + BG(I,1) * T + CG(I,1) * T^2;
   			DHI = 1.8 * DHI;
   			DH = DH + Z(I,1) * XMW(I,1) * DHI;
		end
	end
	if ITY == 2
    	R3 = 3.57696/1000000;
    	DH = DH - R3;
  	end
	%
	% End of TESUB3 substitution.
	%
    DT = -ERR / DH;
    T = T + DT;
	if abs(DT) < 1*10^(-12)
    	break
  	end
	if J == 100
  		T = TIN;
  	end
end
TCS = T;
% End of TESUB2.

      TKS = TCS+273.15;
      
%     TESUB2(XLC,TCC,ESC,0);
%	  TESUB2(Z,T,H,ITY)
% Substitution of TESUB2 for call statement. MWB
% TESUB2 appears to update T and only T.
Z = XLC;
T = TCC;
H = ESC;
ITY = 0;
TIN=T;
for J = 1:100
	%      TESUB1(Z,T,HTEST,ITY);
	%	   TESUB1VAL = TESUB1(Z,T,H,ITY)	
	% Substitution of TESUB1 for call statement MWB
	%
	if ITY == 0
		HTEST = 0;
		for I = 1:8
			HI = T * (AH(I,1) + BH(I,1) * T / 2 + CH(I,1) * T^2 / 3);
			HI = 1.8 * HI;
			HTEST = HTEST + Z(I,1) * XMW(I,1) * HI;
		end
	else
		HTEST = 0;
		for I = 1:8 
			HI = T * (AG(I,1) + BG(I,1) * T / 2 + CG(I,1) * T^2 / 3);
			HI = 1.8 * HI;
			HI = HI + AV(I,1);
			HTEST = HTEST + Z(I,1) * XMW(I,1) * HI;
		end
	end
	if ITY == 2
		R1 = 3.57696 / 1000000;
		HTEST = HTEST - R1 * (T + 273.15);
	end
	%
	% End of TESUB1 sub
	%
	ERR = HTEST - H;
	% 		TESUB3VAL = TESUB3(Z,T,DH,ITY)
	%		TESUB3VAL = TESUB3(Z,T,DH,ITY)
	% Substitution of TESUB3 for call statement MWB
	%
	if ITY == 0
		DH = 0.0;
    	for I = 1:8
			DHI = AH(I,1) + BH(I,1) * T + CH(I,1) * T^2;
			DHI = 1.8 * DHI;
			DH = DH + Z(I,1) * XMW(I,1) * DHI;
		end
	else
		DH=0.0;
		for I = 1:8
			DHI = AG(I,1) + BG(I,1) * T + CG(I,1) * T^2;
			DHI = 1.8 * DHI;
			DH = DH + Z(I,1) * XMW(I,1) * DHI;
		end
	end
	if ITY == 2
		R3 = 3.57696/1000000;
		DH = DH - R3;
	end
	%
	% End of TESUB3 substitution.
	%
	DT = -ERR / DH;
	T = T + DT;
	if abs(DT) < 1*10^(-12)
		break
	end
	if J == 100
		T = TIN;
	end
end
TCC = T;
% End of TESUB2.
      
%     TESUB2(XVV,TCV,ESV,2);
%	  TESUB2(Z,T,H,ITY)
% Substitution of TESUB2 for call statement. MWB
% TESUB2 appears to update T and only T.
Z = XVV;
T = TCV;
H = ESV;
ITY = 2;
TIN=T;
for J = 1:100
	%      TESUB1(Z,T,HTEST,ITY);
	%	   TESUB1VAL = TESUB1(Z,T,H,ITY)	
	% Substitution of TESUB1 for call statement MWB
	%
	if ITY == 0
		HTEST = 0;
		for I = 1:8
			HI = T * (AH(I,1) + BH(I,1) * T / 2 + CH(I,1) * T^2 / 3);
			HI = 1.8 * HI;
			HTEST = HTEST + Z(I,1) * XMW(I,1) * HI;
		end
	else
		HTEST = 0;
		for I = 1:8 
			HI = T * (AG(I,1) + BG(I,1) * T / 2 + CG(I,1) * T^2 / 3);
			HI = 1.8 * HI;
			HI = HI + AV(I,1);
			HTEST = HTEST + Z(I,1) * XMW(I,1) * HI;
		end
	end
	if ITY == 2
		R1 = 3.57696 / 1000000;
		HTEST = HTEST - R1 * (T + 273.15);
	end
	%
	% End of TESUB1 sub
	%
	ERR = HTEST - H;
	% 		TESUB3VAL = TESUB3(Z,T,DH,ITY)
	%		TESUB3VAL = TESUB3(Z,T,DH,ITY)
	% Substitution of TESUB3 for call statement MWB
	%
	if ITY == 0
		DH = 0.0;
    	for I = 1:8
			DHI = AH(I,1) + BH(I,1) * T + CH(I,1) * T^2;
			DHI = 1.8 * DHI;
			DH = DH + Z(I,1) * XMW(I,1) * DHI;
		end
	else
		DH=0.0;
    	for I = 1:8
			DHI = AG(I,1) + BG(I,1) * T + CG(I,1) * T^2;
			DHI = 1.8 * DHI;
			DH = DH + Z(I,1) * XMW(I,1) * DHI;
		end
	end
	if ITY == 2
		R3 = 3.57696/1000000;
		DH = DH - R3;
	end
	%
	% End of TESUB3 substitution.
	%
	DT = -ERR / DH;
	T = T + DT;
	if abs(DT) < 1*10^(-12)
		break
	end
	if J == 100
		T = TIN;
	end
end
TCV = T;
% End of TESUB2.
	  
TKV = TCV+273.15;
      
%	  TESUB4(XLR,TCR,DLR);
%	  TESUB4(X,T,R)
% Substitution of TESUB4 for call statement. MWB
X = XLR;
T4 = TCR;
R4 = DLR;
V = 0.0;
for I = 1:8
	V = V + X(I,1) * XMW(I,1) / (AD(I,1) + (BD(I,1) + CD(I,1) * T4) * T4);
end
R4 = 1.0 / V;
DLR = R4;
% End of TESUB4.

      
%	  TESUB4(XLS,TCS,DLS);
%	  TESUB4(X,T,R)
% Substitution of TESUB4 for call statement. MWB
X = XLS;
T4 = TCS;
R4 = DLS;
V = 0.0;
for I = 1:8
	V = V + X(I,1) * XMW(I,1) / (AD(I,1) + (BD(I,1) + CD(I,1) * T4) * T4);
end
R4 = 1.0 / V;
DLS = R4;
% End of TESUB4.
      
%	  TESUB4(XLC,TCC,DLC);
%	  TESUB4(X,T,R)
% Substitution of TESUB4 for call statement. MWB
X = XLC;
T4 = TCC;
R4 = DLC;
V = 0.0;
for I = 1:8
      V = V + X(I,1) * XMW(I,1) / (AD(I,1) + (BD(I,1) + CD(I,1) * T4) * T4);
end
R4 = 1.0 / V;
DLC = R4;
% End of TESUB4.
VLR = UTLR/DLR;
VLS = UTLS/DLS;
VLC = UTLC/DLC;
VVR = VTR-VLR;
VVS = VTS-VLS;
RG = 998.9;
PTR = 0.0;
PTS = 0.0;
for I = 1:3
      PPR(I,1) = UCVR(I,1)*RG*TKR/VVR;
      PTR = PTR+PPR(I,1);
      PPS(I,1) = UCVS(I,1)*RG*TKS/VVS;
      PTS = PTS+PPS(I,1);
  end
% Some precision error with PTS