% Script to plot the saved data from a dynamic simulation of a H2P cycle
% 16.10.2019.
% Figure(1):power; figure(2):valve; figure(3):fuel
% Disturbance setpoint change

clc;clear;close all
%% Set default options for plotting. You can change it to your preferences
set(0,'DefaultTextFontName','Times','DefaultTextFontSize',18,...
'DefaultAxesFontName','Times','DefaultAxesFontSize',18,...
'DefaultLineLineWidth',4,'DefaultLineMarkerSize',7.75,...
'DefaultTextInterpreter','latex')
%% Define new colors

green    = [0       0.5      0];
gold     = [0.9805, 0.7773,  0.25];
hotpink  = [0.8945, 0.4102,  0.5898];
lilac    = [0.7383, 0.5898,  0.7930];
indigo   = [0.293 0 0.5078]; 
cerulean = [0.1094, 0.6641,  0.8359];

R        = 8.314462618*1e-5;            % [m3*bar/K/mol] ideal gas constant
Mwt      = 18*1e-3;                             % [kg/mol] water molar mass
MSh0     = 100;                                   % [kg] superheater holdup
pSh0     = 8.805290778575733e+01;            % [bar] superheater pressure
pShsp    = pSh0;
TSh0     = 868;
VSh      = MSh0/Mwt*R*TSh0/pSh0;  
zV0      = 0.9;
tsim     = 300;
%% Load data for floating pressure
filename = 'FloatingPressure' ;    
load(filename)
% Inputs: umG - fuel used
Pm   = v_save(:,end);                        %[kW] - mechanical power
MSh  = v_save(:,5);                                  % superheater holdup
TSh  = v_save(:,3);                             % superheater temperature
pSh  = MSh./Mwt.*R.*TSh./VSh;                % [bar] superheater pressure

t    = t_save;                                            % simulation time  
Pmsp = ones(length(t),1)*Pmsp;                         % create step vector

mG0      = ctrlPm(1); KcG=ctrlPm(2); tauIPm=ctrlPm(3);   % power controller
errP     = Pm-Pmsp;                             % error from power setpoint
interrPm = v_save(:,11);                        % integrate error for power
mG       = mG0-KcG*errP-KcG/tauIPm*interrPm;    % [kg/s] flue gass flowrate

Pm0 = Pm0/1e3;                                  % [kW] nominal power output
Pm  = Pm/1e3;                                           % [kW] power output
%% Figure with all simulation
figure(1)% power
plot([-50;0;t],[Pm0;Pm0;Pm],'Color',green)
hold on

figure(2)%pressure
plot([-50;0;t],[pSh0;pSh0;pSh],'Color',green)
hold on

figure(3)% valve 
plot([-50;0;t],[1;1;ones(length(t),1)]*zV0,'Color',green) 
hold on

figure(4)% fuel
plot([-50;0;t],[mG0;mG0;mG,],'Color',green) 
hold on



%% Load data for boiler driven
clear t mG zV Pm pSh KcPm tauIPm errPm interPm Pmsp
filename = 'BoilerDriven' ;    
load(filename)
% Inputs: mG - fuel used; zV valve position
Pm     = v_save(:,end);                             %[kW] - mchanical power
MSh    = v_save(:,5);                                  % superheater holdup
TSh    = v_save(:,3);                             % superheater temperature
pSh    = MSh./Mwt.*R.*TSh./VSh;                % [bar] superheater pressure

t    = t_save;                                            % simulation time  
Pmsp = ones(length(t),1)*Pmsp; % create step vector

mG0      = ctrlPm(1); KcG=ctrlPm(2); tauIPm=ctrlPm(3);    % power controller
errPm    = Pm-Pmsp;                  % error from mechanical power setpoint
interrPm = v_save(:,12);                           % integrated error power
mG       = mG0-KcG*errPm-KcG/tauIPm*interrPm; % [kg/s] flue gass flowrate

Kcz      = ctrlpsh(2); tauIz=ctrlpsh(3); tauTzV = tauIz; % pressure control
errp     = pSh-pShsp;
interrp  = v_save(:,11);                        % integrated error pressure
interrzV = v_save(:,16);        % integrated error valve position. tracking
zV       = zV0-Kcz*errp-Kcz/tauIz*interrp-1/tauTzV*interrzV; 
zVp      = min(1,max(zV,0));                         % valve is opened [0,1]
 
Pm0 = Pm0/1e3;                                  % [kW] nominal power output
Pm  = Pm/1e3;                                           % [kW] power output
%% Figure with all simulation
figure(1)% power
plot([-50;0;t],[Pm0;Pm0;Pm],'Color',hotpink)
hold on
xlabel('Time, [s]')
ylabel('Power, [kW]')
xlim([-50,max(t)])

figure(2)%pressure
plot([-50;0;t],[pSh0;pSh0;pSh],'Color',hotpink)
hold on

figure(3)% valve 
plot([-50;0;t],[zV0;zV0;zVp],'Color',hotpink) 
hold on


figure(4)% fuel
plot([-50;0;t],[mG0;mG0;mG],'Color',hotpink) 
hold on


%% Load data for turbine driven
clear t umG uzV Pm pSh
filename = 'TurbineDriven' ;    
load(filename)
% Inputs: mG - fuel used; zV - steam valve
Pm     = v_save(:,end);                             %[kW] - mchanical power
MSh    = v_save(:,5);                                  % superheater holdup
TSh    = v_save(:,3);                             % superheater temperature
pSh    = MSh./Mwt.*R.*TSh./VSh;                % [bar] superheater pressure

t    = t_save;                                            % simulation time  
Pmsp = ones(length(t),1)*Pmsp; % create step vector

mG0      = ctrlG(1); KcG=ctrlG(2); tauIG=ctrlG(3);    % power controller
errp     = pSh-pShsp;
interrp  = v_save(:,11);                        % integrated error pressure
mG       = mG0-KcG*errp-KcG/tauIG*interrp; % [kg/s] flue gass flowrate

Kcz      = ctrlz(2); tauIz=ctrlz(3); tauTzV = tauIz; % pressure control
errPm    = Pm-Pmsp;                  % error from mechanical power setpoint
interrPm = v_save(:,12);                           % integrated error power
interrzV = v_save(:,16);        % integrated error valve position. tracking
zV       = zV0-Kcz*errPm-Kcz/tauIz*interrPm-1/tauTzV*interrzV; 
zVp      = min(1,max(zV,0));                         % valve is opened [0,1]

Pm0 = Pm0/1e3;                                  % [kW] nominal power output
Pm  = Pm/1e3;                                           % [kW] power output
%% Figure with all simulation
figure(1)% power
plot([-50;0;t],[Pm0;Pm0;Pm],'Color',indigo)
hold on

figure(2)%pressure
plot([-50;0;t],[pSh0;pSh0;pSh],'Color',indigo)
hold on

figure(3)% valve 
plot([-50;0;t],[zV0;zV0;zVp],'Color',indigo) 
hold on

figure(4)% fuel
plot([-50;0;t],[mG0;mG0;mG],'Color',indigo) 
hold on


%% Load data for valve position controller
clear t umG uzV Pm pSh
filename = 'VPC' ;    
load(filename)
% Inputs: mG - fuel used; zV-steam valve
Pm   = v_save(:,end);                               %[kW] - mchanical power
MSh  = v_save(:,5);                                    % superheater holdup
TSh  = v_save(:,3);                               % superheater temperature
pSh  = MSh./Mwt.*R.*TSh./VSh;                  % [bar] superheater pressure

t    = t_save;                                            % simulation time  
Pmsp = ones(length(t),1)*Pmsp; % create step vector

KczPm     = ctrlzPm(2); tauIzPm = ctrlzPm(3);tauTzV = tauIzPm; %power ctrl
errP      = Pm-Pmsp;                        % deviation from power setpoint
interrPm  = v_save(:,11);
interruzV = v_save(:,16); 
zV        = zV0-KczPm*errP-KczPm/tauIzPm*interrPm-1/tauTzV*interruzV; 
zVp       = min(1,max(zV,0));                       % valve is opened [0,1]

KcGz      = ctrlGz(2); tauIGz = ctrlGz(3);                            % VPC
erryzV    = zV-zV0;                        % error from setpoint - for VPC 
interryzV = v_save(:,12);
mG        = mG0-KcGz*erryzV-KcGz/tauIGz*interryzV;%[kg/s] flue gass flowrate

Pm0 = Pm0/1e3;                                  % [MW] nominal power output
Pm  = Pm/1e3;                                           % [MW] power output
%% Figure with all simulation
figure(1)% power
plot([-50;0;t],[Pm0;Pm0;Pm],'Color',gold)
hold on

figure(2)%pressure
plot([-50;0;t],[pSh0;pSh0;pSh],'Color',gold)
hold on

figure(3)% valve 
plot([-50;0;t],[zV0;zV0;zVp],'Color',gold) 
hold on

figure(4)% fuel
plot([-50;0;t],[mG0;mG0;mG],'Color',gold) 
hold on

%% Load data for 2 controllers
clear t umG uzV Pm pSh
filename = '2Controllers' ;    
load(filename)
% Inputs: mG - fuel used; zV - steam valve position
Pm   = v_save(:,end);                               %[kW] - mchanical power
MSh  = v_save(:,5);                                    % superheater holdup
TSh  = v_save(:,3);                               % superheater temperature
pSh  = MSh./Mwt.*R.*TSh./VSh;                  % [bar] superheater pressure

t    = t_save;                                            % simulation time  
Pmsp = ones(length(t),1)*Pmsp; % create step vector

KcGPm    = ctrlG(2); tauIGPm = ctrlG(3);
errPm    = Pm-Pmsp;                                           % power error
interrPm = v_save(:,11);                           % integrated error power
mG       = mG0-KcGPm*errPm-KcGPm/tauIGPm*interrPm;% [kg/s] flue gass flowrate
 
KczPm   = ctrlz(2);
zV      = zV0-KczPm*errPm;       % P-control for steam valve. no antiwindup
zVp     = min(1,max(zV,0));

Pm0  = Pm0/1e3;                                 % [MW] nominal power output
Pm   = Pm/1e3;                                          % [MW] power output                                  
Pmsp = Pmsp/1e3;                                      % [MW] power setpoint
%% Figure with all simulation
figure(1)% power
plot([-50;0;t],[Pm0;Pm0;Pm],'Color',cerulean)
plot([-50;0;t],[Pm0;Pm0;Pmsp],'k--') % plot setpoint
xlabel('Time, [s]')
ylabel('CV7 (Power), [MW]')
xlim([-50,300])
%ylim([14.5 17])
textlegend = {'Floating Pressure','Boiler Driven','Turbine Driven',...
       'VPC','2 Controlles','Setpoint'};
legend(textlegend,'Location','northeast','NumColumns',1)
legend boxoff
set(gcf,'Position',[100 100 750 500])                 % set the figure size                 
print('-depsc2','-r600','Power_dist.eps')                      % save eps 

figure(2)%pressure
plot([-50;0;t],[pSh0;pSh0;pSh],'Color',cerulean)
plot([-50;0;t],[pSh0;pSh0;ones(length(t),1)*pSh0],'k--') % plot setpoint
xlim([-50,300])
%ylim([78 92])
xlabel('Time, [s]')
ylabel('CV8 (Pressure), [bar]')
textlegend = {'Floating Pressure','Boiler Driven','Turbine Driven',...
       'VPC','2 Controlles','Setpoint'};
legend(textlegend,'Location','northeast','NumColumns',1)
legend boxoff
set(gcf,'Position',[100 100 750 500])                 % set the figure size                 
print('-depsc2','-r600','Pressure_dist.eps')                      % save eps 


figure(3)% valve 
plot([-50;0;t],[zV0;zV0;zVp],'Color',cerulean) 
hold on
xlabel('Time, [s]')
ylabel('MV2 input usage, [-]')
xlim([-50,300])
ylim([0.75 1])
textu = {'Floating Pressure','Boiler Driven','Turbine Driven','VPC','2 Controlles'};
legend(textu,'Location','northeast','NumColumns',1)
legend boxoff
set(gcf,'Position',[100 100 750 500])                 % set the figure size                 
print('-depsc2','-r600','Valve_dist.eps')                      % save eps 


figure(4)% fuel
plot([-50;0;t],[mG0;mG0;mG],'Color',cerulean) 
hold on
xlabel('Time, [s]')
ylabel('MV1 input usage, [kg s$^{-1}$]')
xlim([-50,300])
%ylim([26 38])
legend(textu,'Location','northeast','NumColumns',1)
set(gcf,'Position',[100 100 750 500]) 
legend boxoff% set the figure size                 
print('-depsc2','-r600','Fuel_dist.eps')                      % save eps 

%% Plot flue gas temperature
TG0 = ones(length(t+2),1)*(1000+273);
TG1 = TG0+10;

figure
plot([-50;0;t],[TG0(1);TG0(1);TG1],'red')
ylabel('Hot gas temperature [K]')
xlabel('Time, [s]')
xlim([-50,300])
ylim([1270 1285])
set(gcf,'Position',[100 100 750 500]) 
print('-depsc2','-r600','FlueGasT.eps')                      % save eps 



