📄 tesfunc.m
字号:
FCM(I,2)=XST(I,2)*FTM(2,1);
FCM(I,3)=XST(I,3)*FTM(3,1);
FCM(I,4)=XST(I,4)*FTM(4,1);
FCM(I,6)=XST(I,6)*FTM(6,1);
FCM(I,8)=XST(I,8)*FTM(8,1);
FCM(I,9)=XST(I,9)*FTM(9,1);
FCM(I,10)=XST(I,10)*FTM(10,1);
FCM(I,11)=XST(I,11)*FTM(11,1);
FCM(I,13)=XST(I,13)*FTM(13,1);
end
if FTM(11,1) > 0.1
if TCC > 170
TMPFAC=TCC-120.262;
elseif TCC < 5.292
TMPFAC=0.1;
else
TMPFAC=363.744/(177.-TCC)-2.22579488;
end
VOVRL=FTM(4,1)/FTM(11,1)*TMPFAC;
SFR(4,1)=8.5010*VOVRL/(1.0+8.5010*VOVRL);
SFR(5,1)=11.402*VOVRL/(1.0+11.402*VOVRL);
SFR(6,1)=11.795*VOVRL/(1.0+11.795*VOVRL);
SFR(7,1)=0.0480*VOVRL/(1.0+0.0480*VOVRL);
SFR(8,1)=0.0242*VOVRL/(1.0+0.0242*VOVRL);
else
SFR(4,1)=0.9999;
SFR(5,1)=0.999;
SFR(6,1)=0.999;
SFR(7,1)=0.99;
SFR(8,1)=0.98;
end
for I=1:8
FIN(I,1)=0.0;
FIN(I,1)=FIN(I,1)+FCM(I,4);
FIN(I,1)=FIN(I,1)+FCM(I,11);
end
FTM(5,1)=0.0;
FTM(12,1)=0.0;
for I=1:8
FCM(I,5)=SFR(I,1)*FIN(I,1);
FCM(I,12)=FIN(I,1)-FCM(I,5);
FTM(5,1)=FTM(5,1)+FCM(I,5);
FTM(12,1)=FTM(12,1)+FCM(I,12);
end
for I=1:8
XST(I,5)=FCM(I,5)/FTM(5,1);
XST(I,12)=FCM(I,12)/FTM(12,1);
end
TST(5,1)=TCC;
TST(12,1)=TCC;
% TESUB1(XST(1,5),TST(5,1),HST(5,1),1);
% TESUB1(Z,T,H,ITY)
% Substitution of TESUB1 for call statement. MWB
Z(1:8,1)=XST(1:8,5);
T=TST(5,1);
H=HST(5,1);
ITY=1;
if ITY == 0
H = 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;
H = H + Z(I,1) * XMW(I,1) * HI;
end
else
H = 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);
H = H + Z(I,1) * XMW(I,1) * HI;
end
end
if ITY == 2
R1 = 3.57696 / 1000000;
H = H - R1 * (T + 273.15);
end
HST(5,1) = H;
% end of TESUB1
% TESUB1(XST(1,12),TST(12,1),HST(12,1),0);
% TESUB1(Z,T,H,ITY)
% Substitution of TESUB1 for call statement. MWB
Z(1:8,1)=XST(1:8,12);
T=TST(12,1);
H=HST(12,1);
ITY=0;
if ITY == 0
H = 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;
H = H + Z(I,1) * XMW(I,1) * HI;
end
else
H = 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);
H = H + Z(I,1) * XMW(I,1) * HI;
end
end
if ITY == 2
R1 = 3.57696 / 1000000;
H = H - R1 * (T + 273.15);
end
HST(12,1) = H;
% end of TESUB1
FTM(7,1)=FTM(6,1);
HST(7,1)=HST(6,1);
TST(7,1)=TST(6,1);
for I=1:8
XST(I,7)=XST(I,6);
FCM(I,7)=FCM(I,6);
end
if VLR/7.8 > 50.0
UARLEV=1.0;
elseif VLR/7.8 < 10.0
UARLEV=0.0;
else
UARLEV=0.025*VLR/7.8-0.25;
end
UAR=UARLEV*(-0.5*AGSP^2+2.75*AGSP-2.5)*855490.E-6;
QUR=UAR*(TWR-TCR)*(1-0.35*TESUB8(10,TIME));
UAS=0.404655*(1.0-1.0/(1.0+(FTM(8,1)/3528.73)^4));
QUS=UAS*(TWS-TST(8,1))*(1-0.25*TESUB8(11,TIME));
QUC=0;
if TCC < 100
QUC=UAC*(100.0-TCC);
end
XMEAS(1,1)=FTM(3,1)*0.359/35.3145;
XMEAS(2,1)=FTM(1,1)*XMWS(1,1)*0.454;
XMEAS(3,1)=FTM(2,1)*XMWS(2,1)*0.454;
XMEAS(4,1)=FTM(4,1)*0.359/35.3145;
XMEAS(5,1)=FTM(9,1)*0.359/35.3145;
XMEAS(6,1)=FTM(6,1)*0.359/35.3145;
XMEAS(7,1)=(PTR-760.0)/760.0*101.325;
XMEAS(8,1)=(VLR-84.6)/666.7*100.0;
XMEAS(9,1)=TCR;
XMEAS(10,1)=FTM(10,1)*0.359/35.3145;
XMEAS(11,1)=TCS;
XMEAS(12,1)=(VLS-27.5)/290.0*100.0;
XMEAS(13,1)=(PTS-760.0)/760.0*101.325;
XMEAS(14,1)=FTM(11,1)/DLS/35.3145;
XMEAS(15,1)=(VLC-78.25)/VTC*100.0;
XMEAS(16,1)=(PTV-760.0)/760.0*101.325;
XMEAS(17,1)=FTM(13,1)/DLC/35.3145;
XMEAS(18,1)=TCC;
XMEAS(19,1)=QUC*1.04E3*0.454;
XMEAS(20,1)=CPDH*0.0003927E6;
XMEAS(20,1)=CPDH*0.29307E3;
XMEAS(21,1)=TWR;
XMEAS(22,1)=TWS;
ISD=0;
if XMEAS(7,1) > 3000.0
ISD=1;
end
if VLR/35.3145 > 24.0
ISD=1;
end
if VLR/35.3145 < 2.0
ISD=1;
end
if XMEAS(9,1) > 175.0
ISD=1;
end
if VLS/35.3145 > 12.0
ISD=1;
end
if VLS/35.3145 < 1.0
ISD=1;
end
if VLC/35.3145 > 8.0
ISD=1;
end
if VLC/35.3145 < 1.0
ISD=1;
end
if TIME > 0.0 & ISD == 0
for I=1:22
% TESUB6(XNS(I,1),XMNS)
% TESUB6(STD,X)
% Substitution of TESUB6 for call statement. MWB
STD = XNS(I,1);
X = XMNS;
X = 0;
for Iter =1:12
X = X + TESUB7(Iter);
end
X = (X - 6) * STD;
XMNS = X;
% End of TESUB6.
XMEAS(I,1)=XMEAS(I,1)+XMNS;
end
end
XCMP(23,1)=XST(1,7)*100.0;
XCMP(24,1)=XST(2,7)*100.0;
XCMP(25,1)=XST(3,7)*100.0;
XCMP(26,1)=XST(4,7)*100.0;
XCMP(27,1)=XST(5,7)*100.0;
XCMP(28,1)=XST(6,7)*100.0;
XCMP(29,1)=XST(1,10)*100.0;
XCMP(30,1)=XST(2,10)*100.0;
XCMP(31,1)=XST(3,10)*100.0;
XCMP(32,1)=XST(4,10)*100.0;
XCMP(33,1)=XST(5,10)*100.0;
XCMP(34,1)=XST(6,10)*100.0;
XCMP(35,1)=XST(7,10)*100.0;
XCMP(36,1)=XST(8,10)*100.0;
XCMP(37,1)=XST(4,13)*100.0;
XCMP(38,1)=XST(5,13)*100.0;
XCMP(39,1)=XST(6,13)*100.0;
XCMP(40,1)=XST(7,13)*100.0;
XCMP(41,1)=XST(8,13)*100.0;
if TIME == 0
for I=23:41
XDEL(I,1)=XCMP(I,1);
XMEAS(I,1)=XCMP(I,1);
end
TGAS=0.1;
TPROD=0.25;
end
if TIME >= TGAS
for I=23:36
XMEAS(I,1)=XDEL(I,1);
% TESUB6(XNS(I,1),XMNS);
% TESUB6(STD,X)
% Substitution of TESUB6 for call statement. MWB
STD = XNS(I,1);
X = XMNS;
X = 0;
for Iter=1:12
X = X + TESUB7(Iter);
end
X = (X - 6) * STD;
XMNS = X;
% End of TESUB6
XMEAS(I,1)=XMEAS(I,1)+XMNS;
XDEL(I,1)=XCMP(I,1);
end
TGAS=TGAS+0.1;
end
if TIME >= TPROD
for I=37:41
XMEAS(I,1)=XDEL(I,1);
% TESUB6(XNS(I,1),XMNS)
% TESUB6(STD,X)
% Substitution of TESUB6 for call statement. MWB
STD = XNS(I,1);
X = XMNS;
X = 0;
for Iter=1:12
X = X + TESUB7(Iter);
end
X = (X - 6) * STD;
XMNS = X;
% End of TESUB6
XMEAS(I,1)=XMEAS(I,1)+XMNS;
XDEL(I,1)=XCMP(I,1);
end
TPROD=TPROD+0.25;
end
%
for iter = 1:8
YP(iter,1) = FCM(iter,7) - FCM(iter,8) + CRXR(iter,1);
YP(iter+9,1) = FCM(iter,8) - FCM(iter,9) - FCM(iter,10) - FCM(iter,11);
YP(iter+18,1) = FCM(iter,12) - FCM(iter,13);
YP(iter+27,1) = FCM(iter,1) + FCM(iter,2) + FCM(iter,3) + FCM(iter,5) + FCM(iter,9) - FCM(iter,6);
end
% Here is the "correct" version of the separator energy balance:
%
% YP(18)=HST(8)*FTM(8)-(HST(9)*FTM(9)-cpdh)-HST(10)*FTM(10)-HST(11)*FTM(11)+QUS
%
% Here is the original version
YP(9,1)=HST(7,1) * FTM(7,1) - HST(8) * FTM(8) + RH + QUR;
YP(18,1)=HST(8,1)*FTM(8,1)-HST(9,1)*FTM(9,1)-HST(10,1)*FTM(10,1)-HST(11,1)*FTM(11,1)+QUS;
YP(27,1)=HST(4,1)*FTM(4,1)+HST(11,1)*FTM(11,1)-HST(5,1)*FTM(5,1)-HST(13,1)*FTM(13,1)+QUC;
YP(36,1)=HST(1,1)*FTM(1,1)+HST(2,1)*FTM(2,1)+HST(3,1)*FTM(3,1)+HST(5,1)*FTM(5,1)+HST(9,1)*FTM(9,1)-HST(6,1)*FTM(6,1);
YP(37,1)=(FWR*500.53*(TCWR-TWR)-QUR*1000000/1.8)/HWR;
YP(38,1)=(FWS*500.53*(TCWS-TWS)-QUS*1000000/1.8)/HWS;
IVST(10,1)=IDV(14,1);
IVST(11,1)=IDV(15,1);
IVST(5,1)=IDV(19,1);
IVST(7,1)=IDV(19,1);
IVST(8,1)=IDV(19,1);
IVST(9,1)=IDV(19,1);
for I=1:12
if TIME == 0 | abs(VCV(I,1)-XMV(I,1)) > VST(I,1)*IVST(I,1)
VCV(I,1)=XMV(I,1);
end
if VCV(I,1) < 0.0
VCV(I,1)=0.0;
end
if VCV(I,1) > 100.0
VCV(I,1)=100.0;
end
YP(I+38,1)=(VCV(I,1)-VPOS(I,1))/VTAU(I,1);
end
if TIME > 0.0 & ISD ~= 0
for I=1:NN
YP(I,1)=0.0;
end
end
% This is the final output of the subroutine:
% Everything up to this point (except what has been
% documented in the code) runs very similar to the FORTRAN 77
% version. MWB
% 7/12/98
% Matlab: FORTRAN77:
% YY
% 1.040491389000000e+01 1.040491389000000e+01
% 4.363996017000000e+00 4.363996017000000e+00
% 7.570059737000000e+00 7.570059737000000e+00
% 4.230042431000000e-01 4.230042431000000e-01
% 2.415513437000000e+01 2.415513437000000e+01
% 2.942597645000000e+00 2.942597645000000e+00
% 1.543770655000000e+02 1.543770655000000e+02
% 1.591865960000000e+02 1.591865960000000e+02
% 2.808522723000000e+00 2.808522723000000e+00
% 6.375581199000000e+01 6.375581199000000e+01
% 2.674026066000000e+01 2.674026066000000e+01
% 4.638532432000000e+01 4.638532432000000e+01
% 2.464521543000000e-01 2.464521543000000e-01
% 1.520484404000000e+01 1.520484404000000e+01
% 1.852266172000000e+00 1.852266172000000e+00
% 5.244639459000000e+01 5.244639459000000e+01
% 4.120394008000000e+01 4.120394008000000e+01
% 5.699317760000000e-01 5.699317760000000e-01
% 4.306056376000000e-01 4.306056376000000e-01
% 7.990620078300001e-03 7.990620078300001e-03
% 9.056036089000000e-01 9.056036089000000e-01
% 1.605425821600000e-02 1.605425821600000e-02
% 7.509759687000001e-01 7.509759687000001e-01
% 8.858285595500000e-02 8.858285595500000e-02
% 4.827726193000000e+01 4.827726193000000e+01
% 3.938459028000000e+01 3.938459028000000e+01
% 3.755297257000000e-01 3.755297257000000e-01
% 1.077562698000000e+02 1.077562698000000e+02
% 2.977250546000000e+01 2.977250546000000e+01
% 8.832481135000000e+01 8.832481135000000e+01
% 2.303929507000000e+01 2.303929507000000e+01
% 6.285848794000000e+01 6.285848794000000e+01
% 5.546318688000000e+00 5.546318688000000e+00
% 1.192244772000000e+01 1.192244772000000e+01
% 5.555448243000000e+00 5.555448243000000e+00
% 9.218489761999999e-01 9.218489761999999e-01
% 9.459927549000000e+01 9.459927549000000e+01
% 7.729698353000001e+01 7.729698353000001e+01
% 6.305263039000000e+01 6.305263039000000e+01
% 5.397970677000000e+01 5.397970677000000e+01
% 2.464355755000000e+01 2.464355755000000e+01
% 6.130192144000000e+01 6.130192144000000e+01
% 2.221000000000000e+01 2.221000000000000e+01
% 4.006374673000000e+01 4.006374673000000e+01
% 3.810034370000000e+01 3.810034370000000e+01
% 4.653415582000000e+01 4.653415582000000e+01
% 4.744573456000000e+01 4.744573456000000e+01
% 4.110581288000000e+01 4.110581288000000e+01
% 1.811349055000000e+01 1.811349055000000e+01
% 5.000000000000000e+01 5.000000000000000e+01
% YP:
% -9.311351132055279e-04 -2.245521557142600e-06 <- Sig Dif
% -3.680818849716161e-07 -3.680818849716200e-07
% -9.313903218526320e-04 -1.341411575595000e-06 <- Sig Dif
% -6.614342966599907e-04 -4.106221069832800e-07 <- Sig Dif
% -2.689011108429895e-04 -1.116128657940900e-06 <- Sig Dif
% -1.312917689233473e-06 -9.964246872051801e-08 <- Sig Dif
% 6.609202557683602e-04 -1.843536381329600e-07 <- Sig Dif
% 2.688889922808357e-04 -5.530856128643800e-08 <- Sig Dif
% 5.904527482059052e-05 -1.406638006074000e-08 <- Sig Dif
% -3.216924339355387e-06 -3.216923774473900e-06
% -1.320602298626739e-06 -1.320602072141200e-06
% -2.292647105051060e-06 -2.292646652080100e-06
% -1.172087313872083e-07 -1.172087100709300e-07
% -1.564706110457337e-06 -1.564705769396800e-06
% -2.046840421598972e-07 -2.046840066327600e-07
% -3.064993734369637e-07 -3.064992597501300e-07
% -1.538568028536247e-07 -1.538566607450800e-07
% -3.970321582613678e-08 -3.970320694435300e-08
% -2.569517931760856e-10 -2.569517931760900e-10
% -2.359501483084614e-13 -2.359501483084600e-13
% -1.614086642121038e-10 -1.614086642121000e-10
% -1.656258463711424e-11 -1.656258463711400e-11
% -1.189860210359939e-09 -1.189860210359900e-09
% -1.897565438113702e-11 -1.897565438113700e-11
% -4.671457531912893e-08 -4.671457531912900e-08
% -1.934833449013240e-08 -1.934833449013200e-08
% -2.352361638813250e-10 -2.352361638813200e-10
% 5.300687917042524e-06 5.300687462295200e-06
% 1.688940415078832e-06 1.688940187705200e-06
% 3.474454388197046e-06 3.474453933449700e-06
% 4.445469699021487e-07 4.445469699021500e-07
% 2.613836613818421e-06 2.613836272757900e-06
% 3.049074734917667e-07 3.049074308592000e-07
% 5.932505473538185e-07 5.932504336669800e-07
% 3.048364050073360e-07 3.048362486879300e-07
% 5.976254513484491e-08 5.976253802941801e-08
% -4.380748794538104e-07 -4.380756709460800e-07
% 3.318385615567584e-07 3.318384779400800e-07
% 0 0
% 0 0
% 0 0
% 0 0
% 0 0
% 0 0
% 0 0
% 0 0
% 0 0
% 0 0
% 0 0
% 0 0
% End of TEFUNC
%disp ('initialized')
x0 = YY;
YP0 = YP;
%
% str is always an empty matrix
%
str = [];
%
% initialize the array of sample times
%
ts = [0 0];
% end mdlInitializeSizes
%
%=============================================================================
% mdlDerivatives
% Return the derivatives for the continuous states.
%=============================================================================
%
function sys=mdlDerivatives(t,x,u)
global YP0 YY YP
TIME = t;
if t == 0
sys = YP0;
else
global XMEAS XMV
global NN
IDV = u(1:20,1);
XMV = u(21:32,1);
% TEFUNC(NN,TIME,YY,YP);
% TEFUNCVAL = TEFUNC(NN,TIME,YY,YP)
%
% Substitution of TEFUNC for call statement. MWB
%
% Tennessee Eastman Process Control Test Problem
%
% James J. Downs and Ernest F. Vogel
%
% Process and Control Systems Engineering
% Tennessee Eastman Company
% P.O. Box 511
% ⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -