%task        : Model of a steam cycle with 1 pressure level
%author      : Cristina Zotica
%contact     : cristina.f.zotica@ntnu.no
%organization: Department of Chemical Engineering, NTNU, Norway
%since       : Jan 21, 2019
%change      : 
%version     : v1.1
%todo        : add drum.
%requires    : 
%note        : Data from simple steam cycle in Skogestad, Prosessteknikk
%components  : Pump(P);Bypass(B) Drum(D); Economizer(E); Condenser(C);
%              Mixer before drum of economizer and its bypass (MB); 
%              Superheater(S);Turbine valve(V), Turbine(T); 
%              generator(G), attemperator(A).

clc;clear;close all

load DesignParam.mat%var: UAE, UAD, UASh, cvD, Kv, phid, m, mB, mP,mA, mG W

%% Physical constants
Tref  = 0;                        % [K] reference tempereature for enthalpy
TrefB = 576; % [K] reference boiling temperature = drum nominal temperature
dHvap = 1382;                    % [kJ/kg] water vaporization heat at 576 K
cpg   = 1.25;                    % [kJ/kg/K] flue gas average heat capacity
cpW   = 4.18;                               % [kJ/kg/K] water specific heat
cpSt  = 3;                                  % [kJ/kg/K] steam specific heat
R     = 8.314462618*1e-5;               % [m3*bar/K/mol] ideal gas constant
Mwt   = 18*1e-3;                                % [kg/mol] water molar mass

%% Initial condition states

TE0   = 576;                            % [K] economizer outlet temperature 
TD0   = TE0;                                  % [K] drum outlet temperature 
TSh0  = 868;    % [K] supeheated steam nominal temperature to calculate VSh
TA0   = 529+273;                       % [K] temperature after attemperator
MSh0  = 100;                              % [kg] superheater nominal holdup
MD0   = 3000;                                    % [kg] drum nominal holdup
M0    = 7.56;         % [pu MW/pu rad /s] angular momentum per machine base. from Wood&Wollenberg 2014. example 10.A
D0    = 2;              % [pu MW/rad/s] electrical damping per machine base. from Wood&Wollenberg 2014. example 10.A
w0    = 50;                       % [Hz] nominal grid and turbine frequency
W0    = W;                                              % [kW] turbine work 
Pm0   = W0;                    % [kW] base load - nominal power production

%% Parameters. process variables kept constant. can be changed to states
T0    = 45+273;                                     % [K] inlet temperature
TGE0  = 150+273;                          % [K] flue gas outlet temperature
ME    = 100;                                       % [kg] economizer holdup
pC    = 0.096;                 % [bar] pressure at turbine outlet. constant
VSh   = MSh0/Mwt*R*TSh0/pSh;                      % [m3] superheater volume

%% Input steps
mP0    = mA+m;                            % [kg/s] pump nominal outlet flow
mB0    = 0;                                      % [kg/s] economizer bypass
mG0    = mG;                                    % [kg/s] flue gas mass flow
mA0    = mA;       % [kg/s] attemperation flowrate MV for superheater steam
zV0    = 0.9;                                    % [-] steam  valve opening
PL     = Pm0;  
%% Controller parameters
%Drum level control
MDsp   = MD0; % drum level setpoint
kD     = 0.97;                 % gain from feedawter flowrate to drum level
thetaD = 0;                                                     % [s] delay
taucD  = 10;                         %[s] desired closed loop time constant
KcD    = 1/kD*1/(taucD+thetaD);   % proportional gain drum level controller
tauID  = 4*(taucD+thetaD);                                  % integral time

%Flue gas temperature controller
TGsp    = 423;                                   % [K] flue gas temperature
kTG     = 9.767;
tauTG   = 10;
thetaTG = 0;
taucTG  = 20;
KcTG    = 1/kTG*tauTG/(taucTG+thetaTG);%proportional gain drum level controller
tauITG   = min(tauTG,4*(taucTG+thetaTG));                               % integral time

%Superheated steam temperature controller
TAsp    = 802;                                               % [K] setpoint
kTA     = -86.5;
tauTA   = 1;
thetaTA = 0;
taucTA  = 15;
KcA     = 1/kTA*tauTA/(taucTA+thetaTA);%proportional gain SH temp controller
tauIA   = min(tauTA,4*(taucTA+thetaTA));                    % integral time

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%    VPC for the steam valve, to use the valve dynamically, and then take it to it's nominal point. 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
zVsp    = 0.9; % steam valve nominal opening. s
kGz     = -0.8; % use the fuel valve to bring the steam valve to nominal
tauGz   = 115; 
thetaGz = 1;
taucGz  = 100;
KcGz    = 1/kGz*tauGz/(taucGz+thetaGz);%proportional gain SH temp controller
tauIGz  = min(tauGz,4*(taucGz+thetaGz));                    % integral time

Pmsp     = Pm0*0.9; % use the steam valve to control the power.
kzPm     = 500; % initial gain from fuel to steam valve. same as TD&2Ctrl
tauzPm   = 1; 
thetazPm = 0;
tauczPm  = 20;
KczPm    = 1/kzPm*tauzPm/(tauczPm+thetazPm);%proportional gain SH temp controller
tauIzPm  = min(tauzPm,4*(tauczPm+thetazPm));                               % integral time

%% Simulation
x0 = [TE0;TD0;TSh0;MD0; MSh0;w0;TA0]; % initial condition- differential states
nx = length(x0);                            % number of differential states
z0 = [TGE0;TGD;TGSh;pT; Pm0];     % initial condition- algebraic states
nz = length(z0);                               % number of algebraic states
v0 = [x0;z0]; 
nv = length(v0);   
u0 = [mG0;zV0;mP0;mB0;mA0];
param = [Tref;TrefB;dHvap;cpg;cpW;cpSt;R;Mwt;...
         T0;TG0;ME;pC;...
         UAE; UAD; UASh; CvD; Kv; VSh; phid;...
         Pm0; w0; M0; D0; PL; nx];                   % vector of parameters
tsim = 10000;                                        % [sec] simulation time

np = length(param);
nu = length(u0);

Kdyn = [eye(nx),zeros(nx,nz)];% selection matrix for differential equations
Kalg = [zeros(nz,nx),eye(nz)];   % selection matrix for algebraic equations

KMD  = zeros(nv,1); KMD(4) = 1;              % select level from all states
KTG  = zeros(nv,1); KTG(8) = 1;      % select flue gas temp from all states
KTA  = zeros(nv,1); KTA(7) = 1;         % select steam temp from all states
KPm  = zeros(nv,1); KPm(12)= 1;              % select power from all states


sim('simH2P')


MSh    = v(:,5); % superheater holdup
TSh    = v(:,3); % superheater temperature
pSh    = MSh./Mwt.*R.*TSh./VSh;                % [bar] superheater pressure
pSh0   = MSh0./Mwt.*R.*TSh0./VSh;              % [bar] superheater pressure
pT     = v(:,end-1);
Pm     = v(:,end); % mechanical power

figure
plot([-200;0],[pSh0;pSh0],'b')
hold on
plot([0;0],[pSh0;pSh(1)],'b')
plot(t,pSh,'b')


figure
plot([-200;0],[Pm0;Pm0],'r')
hold on
plot([0;0],[Pm0;Pm(1)],'r')
plot(t,Pm,'r')%Power. invers response. very small steady-state gain

%% Calculate the input mG and zV

Pm       = v(:,end);
errPm    = (Pm-Pmsp);
interrPm = v(:,11);
uzV      = zV-KczPm*errPm-KczPm\tauIzPm*interrPm;

interrzV = v(:,12);
errzV    = uzV-zVsp;
umG      = mG-KcGz*errzV-KcGz/tauIGz*interrzV;

figure
plot(t,uzV,'r')
hold on
plot([-50,0,0],[zV,zV,uzV(1)],'r')
figure
plot(t,umG,'r')
hold on
plot([-50,0,0],[mG,mG,umG(1)],'r')

%% Save data
promt    = 'Select a file name to save the date:';           % ask the user
filename = input(promt,'s');               % receive the name from the user
folder   = '\\home.ansatt.ntnu.no\cristinz\Documents\phd\Power Cycles\Models\Simple Model\Dynamic\Data/';                            % select folder to save data
save(strcat(folder,filename),'t','v','uzV','umG','Pmsp','ctrl') % save the date in the given filename


