%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       : Oct, 2019
%change      : 
%version     : v1.3. with antiwindup and kpi
%todo        : 
%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
%TG0 = TG0+10;
%% 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
mP    = mA+m;                             % [kg/s] pump nominal outlet flow
mB    = 0;                                       % [kg/s] economizer bypass
mG    = mG;                                     % [kg/s] flue gas mass flow
mA    = mA;        % [kg/s] attemperation flowrate MV for superheater steam
zV    = 0.9;                                     % [-] steam  valve opening
PL    = Pm0;                                       % [kW]; load disturbance
%% 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  = 300;                         %[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
tauIG   = 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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Boiler drive: fuel-load (slow) & pressure-turbine valve (fast-tune first)
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
Pmsp    = Pm0*1.1;
kPm     = 432.33;
tauPm   = 40; 
thetaPm = 3;
taucPm  = 10*thetaPm;
KcPm    = 1/kPm*tauPm/(taucPm+thetaPm);%proportional gain SH temp controller
tauIPm  = min(tauPm,4*(taucPm+thetaPm));                    % integral time

pShsp    = pSh;
kpsh     = -5; % process gain (-0.5/1)
taupsh   = 37; 
thetapsh = 0;
taucpsh  = 5;
Kcz      = 1/kpsh*taupsh/(taucpsh+thetapsh);%proportional gain SH temp controller
tauIz    = min(taupsh,4*(taucpsh+thetapsh));                % integral time

%%%%%%%%%%%%%%%%%% Ctrl vector %%%%%%%%%%%%%%%%%%
ctrlMD  = [mP; KcD; tauID];                          % controller for level
ctrlTA  = [mA; KcA; tauIA];  % controller for superheated steam temperature
ctrlTG  = [mB; KcTG; tauIG];          % controller for flue gas temperature
ctrlPm  = [mG; KcPm; tauIPm];
ctrlpsh = [zV; Kcz; tauIz];
ctrl    = [ctrlMD; ctrlTA; ctrlTG; ctrlPm; ctrlpsh];
%% Simulation
% initial condition- differential states
x0 = [TE0;TD0;TSh0;MD0; MSh0;w0;TA0;0;0;0;0;0;0;0;0;0]; 
nx = length(x0);                            % number of differential states
z0 = [TGE0;TGD;TGSh;pT; Pm0;mG;zV];     % initial condition- algebraic states
nz = length(z0);                               % number of algebraic states
v0 = [x0;z0];                                          % initial conditions
u  = [Pmsp; MDsp;TGsp; TAsp;pShsp];                                % inputs
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

Mass  = blkdiag(eye(nx),zeros(nz));    % mass matrix. 1-diff, 0- alg states
opts  = odeset('Mass',Mass);%,'MaxStep',0.01);         % simulation options

% step changes in power setpoint. +-10 % from nominal with "step" % ramp.
step   = 0.1;                                                  % step size
ymin   = -0.1;                              % min step change. not used now
ymax   = 0;                                             % max step change 
nmin   = -ymin/step;
n      = [-step:-step:ymin, (ymin+step):step:ymax];          % step changes 
%n      = 0;                           % no step change for disturbance test
t_save = [];                                % empty vector to store results
v_save = [];                                % empty vector to store results
tsim   = 300;                                      % [sec] integration time
tspan  = [0;tsim];                                  % [sec] simulation time
for i = 1:length(n)
    Pmsp(i)= Pm0*(1+n(i));                            % [kW] power setpoint
    u      = [Pmsp(i); MDsp;TGsp; TAsp;pShsp];                     % inputs
    [t,v]=ode15s(@(t,v)H2P_Boiler_driven(t,v,u,param,ctrl),tspan,v0,opts);
    t_save = [t_save;t];                                 % save time values
    v_save = [v_save;v];                                % save state values
    v0     = v(end,:);       % new intial conditions = last variable values
    tspan  = t(end)+[0 tsim];     % extend simulation time from last time
    
end

figure(1)
plot(t_save,v_save(:,end)); % plot last state
hold on
%% 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 for setpoint changes
folder   = '\\home.ansatt.ntnu.no\cristinz\Documents\phd\Power Cycles\Models\Simple Model\Dynamic\Antiwindup\Data\Setpoint\';                            % select folder to save data

% folder for disturbance in flue gas inlet temperature
%folder   = '\\home.ansatt.ntnu.no\cristinz\Documents\phd\Power Cycles\Models\Simple Model\Dynamic\Antiwindup\Data\TG0plus10K\';                            % select folder to save data

save(strcat(folder,filename),'t_save','v_save','ctrlPm','ctrlpsh','Pmsp','Pm0') % save the date in the given filename




