% Test 2, Standard Model, Usoro, MIT, 1977
% "A drum boiler-turbine power plant under emergency state control"
% eqm 15feb'19
%------------------------------------------------------------------
% read variables from file
clear
close all
format short e
global dxdt
stdinit775       % modificar de acuerdo a la carga inicial
% Miscellanous initializations
t = -Ts; xdot=zeros(27,1);
%-----------------------------------------------------------------
%-------------------         >>>   MAIN LOOP   <<<           -----
%-----------------------------------------------------------------
for i=1:samples,
  t = t + Ts;
  if(t>=10 && t<=120)           % t<=100 /120
      ldc=3.875-0.0125*(t-10);    % ldc=3.875+0.0125*(t-10);
  end
%jc=jc+1;
%----------------falta---------------
% set controlled input values
card=limchk(card);
cfld=limchk(cfld);
cgrd=limchk(cgrd);
cwryd=limchk(cwryd);
cxggd=limchk(cxggd);
csyd=limchk(csyd);
c12tr=check(c12tr,kn5,kn0);
acv=xducer(kn0,kn5,kavl,kavu,c12tr);
war=xducer(kcl,kcu,kwarl,kwaru,card);
wfl=xducer(kcl,kcu,kwfll,kwflu,cfld);
kwgr=kc2gr*war;
kcwgr=xducer(kwgrl,kwgru,kcl,kcu,kwgr);
if(cgrd>kcwgr); cgrd=kcwgr; end
wgr=xducer(kcl,kcu,kwgrl,kwgru,cgrd);
cry=cwryd*kc1ry;
wry=xducer(kcryl,kcryu,kwryl,kwryu,cry);
xgg=xducer(kcl,kcu,kxggl,kxggu,cxggd);
zxgg=xgg*57.3;
csy=csyd*kc4sy;
wsy=xducer(kcsyl,kcsyu,kwsyl,kwsyu,csy);
%-----------------------------------------------------------------
% state properties
[rdrw,hdrs,hdrw,pdrs,tdrs]=drstat(rdrs);
[psso,tsso,ssso]=shstat(rsso,hsso);
[prho,trho]=rhstat(rrho,hrho);
%-----------------------------------------------------------------
% steam flow
kd1=1;
pshd=pdrs-psso;
if(pshd<kn0); kd1=-1; end
wdrs=kd1*sqrt(abs(pshd)*rdrs/kfsh);
wsso=kcv*acv*psso;
whpo=kwhp*wsso;
wrho=kwrh*prho;
%------------------------------------------------------------------
% turbine and reheater
ehp=0.589+2.317e-4*wsso;
kelp=0.84;
phpo=prho+kfrh*wrho*wrho/rrho;
[hhpo,thpo]=hpstat(ssso,phpo,ehp,hsso);
mwhp=wsso*(hsso-hhpo)*kj;
mwlp=wrho*(hrho-khlpo)*kj*kelp;
mwtro=mwhp+mwlp;
mwo=kmwx*mwtro;
mwtrpu=mwtro/kmwr;
mwgnpu=kn2*sin(delta);
mwgn=mwgnpu*kmwr;
%------------------------------------------------------------------
% heat transfer
wwwg=war+wfl+wgr;
kuwwgm=73.801+1.4173*wfl;
tfn1=kt0+(wfl*khfl*kefl+wfl*ksfl*(ktfl-kt0)+war*ksar*(ktar-kt0)...
    +wgr*ksgr*(ktecgo-kt0))/(wwwg*ksfng);
uxgg=k1xgg+k2xgg*sin(xgg)/cos(xgg);
uwwgm=kuwwgm*uxgg;
zfn1=wwwg*ksfng+knp5*uwwgm;
twwge=(wwwg*ksfng*tfn1+knp5*uwwgm*twwm)/zfn1;
qwwgm=uwwgm*(twwge-twwm);
twwgo=twwge-knp5*qwwgm/(wwwg*ksfng);
wwwg2=wwwg*wwwg;
qgp=wwwg*ksgpe*(twwgo-ktecgo);
zwwwg=wfl+war;
zwwwg2=zwwwg*zwwwg;
zqrh=0.2710-7.803e-6*zwwwg-1.337e-8*zwwwg2;
zqec=0.6912-8.852e-4*zwwwg+3.971e-7*zwwwg2;
qec=zqec*qgp;
qrh=zqrh*qgp;
qsh=qgp-qec-qrh;
qwwmw=kuwwmw*(twwm-tdrs);
%-----------------------------------------------------------------
% drum
zdr1=kdr2+kn2*kdr3*rdrs;
zdr2=kvdr/kvdrw-kn1+zdr1;
zdr3=kvdrs*(hdrs+rdrs*(kdr5+kn2*kdr6*rdrs))...
    +kvdrw*(rdrw*(kdr8+kn2*kdr9*rdrs)+hdrw*zdr1);
zdr4=zdr3-kvdrw*zdr2*heco;
frdrs=(qwwmw+wdrs*(heco-hdrs))/zdr4;
wfw=kvdrw*zdr2*frdrs+wdrs;
%
hcpo=hcnw+khcpd;
hlhd=113.12+28.3e-3*wsso;
hded=43.62+14.817e-3*wsso;
hfpo=hdew+khfpd;
hhhd=136.66+56.2e-3*wsso;
%-----------------------------------------------------------------
% controls - boiler master demand
kcpsso=4.204;
kcntr=4.351;
cntr=xducer(kntrl,kntru,kcl,kcu,ntr);
cpsso=xducer(kpssol,kpssou,kcl,kcu,psso);
c1md=kcpsso-cpsso;
c2md=c1md*kc1md;
c3md=limchk(c3md);
c4md=c2md+c3md;
c4md=limchk(c4md);
c5md=kc2md*(kcntr-cntr);
cbmd=c4md+c5md;
cbmd=limchk(cbmd);
%----------------------------------------------------------------
% air flow control
cwar=xducer(kwarl,kwaru,kcl,kcu,war);
c2ar=cbmd;
if(c2ar<kc2arl); c2ar=kc2arl; end
c3ar=c2ar-cwar;
c4ar=c3ar*kc1ar;
c5ar=limchk(c5ar);
c6ar=c5ar+c4ar;
c6ar=limchk(c6ar);
%----------------------------------------------------------------
% fuel flow control
cwfl=cfld;
c2fl=cbmd;
if(cwar<cbmd); c2fl=cwar; end
c3fl=c2fl-cfld;
c4fl=c3fl*kc1fl;
c5fl=limchk(c5fl);
c6fl=c4fl+c5fl;
c6fl=limchk(c6fl);
%----------------------------------------------------------------
% gas recirculation control
kcxgg=3;
kcgr=0;
if(abs(xgg)>knp087); kcgr=1; end
c1gr=cxggd-kcxgg;
c2gr=limchk(c2gr);
%----------------------------------------------------------------
% reheat temperature control
ktrh=k1trh+k2trh*wsso;
ktrh=check(ktrh,k4trh,k3trh);
kctrh=xducer(ktrhol,ktrhou,kcl,kcu,ktrh);
ctrho=xducer(ktrhol,ktrhou,kcl,kcu,trho);
c3rh=kctrh-ctrho;
c4rh=c3rh*kc1rh;
c5rh=limchk(c5rh);
c6rh=c4rh+c5rh;
c6rh=limchk(c6rh);
cwry=c6rh;
cxgg=c6rh;
%----------------------------------------------------------------
% superheat temperature control
ktss=k1tss+k2tss*wsso;
ktss=check(ktss,k4tss,k3tss);
kctsso=xducer(ktssol,ktssou,kcl,kcu,ktss);
ctsso=xducer(ktssol,ktssou,kcl,kcu,tsso);
c3sy=kctsso-ctsso;
c4sy=c3sy*kc1sy;
c5sy=limchk(c5sy);
c6sy=c4sy+c5sy;
c6sy=limchk(c6sy);
%----------------------------------------------------------------
% turbine control
c1tr=kc2tr*(kcntr-cntr);
cmwtro=xducer(kmwtrl,kmwtru,kn0,kn6,mwtro);
c2tr=kc1tr*(ldc-cmwtro);
c4tr=c1tr/kcvreg;
c6tr=check(c6tr,kn1,km1);
c7tr=c6tr+c2tr;
c8tr=c7tr+ldc;
c8tr=check(c8tr,kn5,kn0);
c9tr=check(c9tr,kn5,kn0);
c10tr=c4tr+c9tr;
%----------------------------------------------------------------
% dynamic equations
% equivalence
% x0=[heco, rdrs, twwm, rsso, hsso, rrho, hrho,...
%     ntr, hlho, hdew, hhho, c3md, c5ar, c5fl,...
%     c2gr, c5rh, c5sy, card, cfld, cgrd, cwryd,...
%     cxggd, csyd, c6tr, c9tr, c12tr, delta];
xdot(1)=(wfw*(hhho-heco)+qec)/kmece;
xdot(2)=frdrs;
xdot(3)=(qwwgm-qwwmw)/(kmwwme*kswwm);
xdot(4)=(wdrs+wsy-wsso)/kvsh;
xdot(5)=(wdrs*hdrs+wsy*hfpo-wsso*hsso+qsh)/kmshe;
xdot(6)=(whpo+wry-wrho)/kvrh;
xdot(7)=(whpo*hhpo+wry*hfpo-wrho*hrho+qrh)/kmrhe;
xdot(8)=(mwtro-mwgn)/(ntr*kjtre);
xdot(9)=(hcpo+hlhd-hlho)/ktclh;
xdot(10)=(hlho+hded-hdew)/ktcde;
xdot(11)=(hfpo+hhhd-hhho)/ktchh;
xdot(12)=c2md/ktc1md;
xdot(13)=c4ar/ktc1ar;
xdot(14)=c4fl/ktc1fl;
xdot(15)=c1gr*kc1gr*kcgr;
xdot(16)=c4rh/ktc1rh;
xdot(17)=c4sy/ktc1sy;
xdot(18)=(c6ar-card)/ktc2ar;
xdot(19)=(c6fl-cfld)/ktc2fl;
xdot(20)=(c2gr-cgrd)/ktc1gr;
xdot(21)=(cwry-cwryd)/ktc2rh;
xdot(22)=(cxgg-cxggd)/ktc3rh;
xdot(23)=(c6sy-csyd)/ktc2sy;
xdot(24)=c2tr/ktc1tr;
xdot(25)=(c8tr-c9tr)/ktc2tr;
xdot(26)=(c10tr-c12tr)/ktc3tr;
xdot(27)=kc1gn*(ntr-kntr);
dxdt=xdot;
%----------------------------------------------------------------
%   elseif strcmp(simul,'matlab')
%     ugl = u_old(1);
%     [time,x] = ode45(mat_model,[t-Ts t],x0);
%     x0 = x(length(time),:)';
%     eval(['y = ' model_out '(x0);']);
%options=odeset('Refine',1);
if(t>=tprint)
xx2=[xx2;t,ntr,mwo,psso,wsso,c3md,c10tr,c6fl,c6ar,...
    hsso,hrho,csyd,cxggd,tsso,trho,twwm,cgrd];
    tprint=tprint+1; % graficas de resultados
end
[time,x]=ode45(mat_model,[t-Ts t],x0);
x0 = x(length(time),:)';
%x=x0+xdot*Ts;
xx=[xx;x0'];
% variables de estado
heco=x0(1); rdrs=x0(2); twwm=x0(3); rsso=x0(4); hsso=x0(5);
rrho=x0(6); hrho=x0(7); ntr=x0(8); hlho=x0(9); hdew=x0(10);
hhho=x0(11); c3md=x0(12); c5ar=x0(13); c5fl=x0(14); c2gr=x0(15);
c5rh=x0(16); c5sy=x0(17); card=x0(18); cfld=x0(19); cgrd=x0(20);
cwryd=x0(21); cxggd=x0(22); csyd=x0(23); c6tr=x0(24); c9tr=x0(25);
c12tr=x0(26); delta=x0(27);
%x0(8)=x00(8); % pba para fijar los controles
end
xx3=[1:12;mwo,ntr,wsso,psso,tsso-459.7,trho-459.7,...
    prho,pdrs,hlho,hhho,war,wfl]'; % steady state values (Table V)
xx4=[(1:27)',x0]; % state variables
%----------------------------------------------------------------
%----------------         >>>   END OF MAIN LOOP   <<<        ---
%----------------------------------------------------------------
figure;
subplot(221); plot(xx2(:,1),xx2(:,2));
legend('turbine speed, rad/s'); legend('boxoff');
axis([0 t 373 381]);
subplot(222); plot(xx2(:,1),xx2(:,3));
legend('power output, MW'); legend('boxoff');
%axis([0 t 460 600]); axis([0 t 280 480]); axis([0 t 300 480]);
axis([0 t 280 480]); %axis([0 t 460 620]);
subplot(223); plot(xx2(:,1),xx2(:,4));
legend('throttle pressure, psia'); legend('boxoff');
%axis([0 t 2380 2470]);axis([0 t 2410 2445]);axis([0 t 2370 2440]);
axis([0 t 2410 2445]); %axis([0 t 2360 2450]);
subplot(224); plot(xx2(:,1),xx2(:,5));
legend('main steam flow, lb/s'); legend('boxoff');
%axis([0 t 800 1150]); axis([0 t 500 850]); axis([0 t 500 850]);
axis([0 t 500 850]); %axis([0 t 800 1150]);
figure;
subplot(221); plot(xx2(:,1),xx2(:,6));
legend('boiler master demand, V'); legend('boxoff');
%axis([0 t 3.7 4.6]); axis([0 t 2.7 3.8]); axis([0 t 2.7 3.8]);
axis([0 t 2.7 3.8]); %axis([0 t 3.6 5.0]);
subplot(222); plot(xx2(:,1),xx2(:,7));
legend('governor valve control, V'); legend('boxoff');
%axis([0 t 3.6 5]); axis([0 t 2.2 3.8]); axis([0 t 2.4 3.8]);
axis([0 t 2.2 3.8]); %axis([0 t 3.6 5.0]);
subplot(223); plot(xx2(:,1),xx2(:,8));
legend('fuel flow control, V'); legend('boxoff');
%axis([0 t 3.6 4.6]); axis([0 t 2.4 3.8]); axis([0 t 2.6 4.0]);
axis([0 t 2.4 3.8]); %axis([0 t 3.6 4.8]);
subplot(224); plot(xx2(:,1),xx2(:,9));
legend('air flow control, V'); legend('boxoff');
%axis([0 t 3.6 4.6]); axis([0 t 2.4 3.8]); axis([0 t 2.6 4.2]);
axis([0 t 2.4 3.8]); %axis([0 t 3.6 5.0]);
figure;
subplot(221); plot(xx2(:,1),xx2(:,10));
legend('superheat enthalpy, BTU/lb'); legend('boxoff');
%axis([0 t 1440 1490]);axis([0 t 1435 1460]);axis([0 t 1445 1480]);
axis([0 t 1435 1460]); %axis([0 t 1435 1485]);
subplot(222); plot(xx2(:,1),xx2(:,11));% 'Location','northoutside'
legend('reheat enthalpy, BTU/lb'); legend('boxoff');  
%axis([0 t 1510 1530]);axis([0 t 1500 1525]);axis([0 t 1500 1540]);
axis([0 t 1500 1525]); %axis([0 t 1514 1532]);
subplot(223); plot(xx2(:,1),xx2(:,12));
legend('superheat temp. control, V'); legend('boxoff');
%axis([0 t 1.0 4.5]); axis([0 t 2.5 5.0]); axis([0 t 1.0 5.0]);
axis([0 t 2.5 5.0]); %axis([0 t 1.5 4.5]);
subplot(224); plot(xx2(:,1),xx2(:,13));
legend('reheat temp. control, V'); legend('boxoff');
%axis([0 t 2.5 5.0]); axis([0 t 2.6 4.4]); axis([0 t 1.5 4.5]);
axis([0 t 2.6 4.4]); %axis([0 t 1.4 2.8]);
figure;
subplot(221); plot(xx2(:,1),xx2(:,14));
legend('superheat temperature, °R'); legend('boxoff');
%axis([0 t 1455 1480]);axis([0 t 1430 1460]);axis([0 t 1440 1490]);
axis([0 t 1430 1460]); %axis([0 t 1448 1470]);
subplot(222); plot(xx2(:,1),xx2(:,15));
legend('reheat temperature, °R'); legend('boxoff');
%axis([0 t 1435 1470]);axis([0 t 1410 1460]);axis([0 t 1410 1490]);
axis([0 t 1410 1460]); %axis([0 t 1450 1480]);
subplot(223); plot(xx2(:,1),xx2(:,16));
legend('waterwall metal temp., °R'); legend('boxoff');
%axis([0 t 1142 1158]);axis([0 t 1120 1160]);axis([0 t 1125 1160]);
axis([0 t 1120 1160]); % axis([0 t 1144 1158]);
subplot(224); plot(xx2(:,1),xx2(:,17));
legend('gas recirc. control, V'); legend('boxoff');
%axis([0 t 2.4 4.0]); axis([0 t 3.65 3.90]); axis([0 t 3.2 4.0]);
axis([0 t 3.65 3.9]); %axis([0 t 2.4 3.8]);