clc;clear;close all
%% Plotting options
set(0,'DefaultTextFontName','Times',...
     'DefaultTextFontSize',14,...
     'DefaultAxesFontName','Times',...
     'DefaultAxesFontSize',14,...
     'DefaultLineLineWidth',4,...
     'DefaultLineMarkerSize',4,...
     'DefaultStairLineWidth',4)
orange  = [1 0.5 0];
green   = [0 0.5 0];
greencg = [83,188,102]/256;

teal    = [67 183 176]/256;
maroon  = [192 0 0]/256;
gold    = [255 192 0]/256;
skyblue = [39 141 251]/256;
indigo   = [0.2930 0 0.5078];
lilac    = [0.7383    0.5898    0.7930];
hotpink  = [0.8945    0.4102    0.5898];
cerulean = [0.1094    0.6641    0.8359];
plum     = [194, 147, 141]/256;
purple   = [176 27 129]/256;

pos     = [100 100 450 300]; % figure position on screen
poslgd  = [0.5888    0.6268    0.377    0.1963];
Vmax    = 5000; % [m3] maximum tank volume
Vmin    = 0;    % [m3] minimum tank volume
qmax    = 1000; % [m3/h] maximum flowrate
qmin    = 0;    % [m3/h] minimum flowrate
tend    = 15;   % [h] plot end time

%% Plot open loop. step responses for qS
load('ResultsOL_qS') % step in hot water supply
Vh=yOL(:,1); p=yOL(:,2);
qS=dOL(:,1); q=dOL(:,2);
qSP=uOL(:,1); qTP=uOL(:,2); qEP=uOL(:,3);
for i =1:length(Vh)
     if   Vh(i)  >= Vmax
          qD(i)  = qS(i)-qSP(i); % dumped
          qST(i) = 0;      % charge
     elseif Vh(i)< Vmax
          qD(i)  = 0;       % dump
          qST(i) = qS(i)-qSP(i)  % charge
     end
end

figure(); % MVs
plot(tOL,qS,'Color',green)
hold on
plot(tOL,q,'Color',skyblue)                                      % demand q
plot(tOL,qSP,'--','Color',gold)                  % qSP from boiler to plant
plot(tOL,qD,'Color','k')                                          % qD dump
plot(tOL,qST,'--','Color',orange );                %qST from boiler to tank
plot(tOL,qTP,'-','Color',teal)                     % qTP from tank to plant
plot(tOL,qEP,'--','Color',purple)       % qEP from electric boiler to plant
lgd=legend('q_S','q','q_{SP}','q_D','q_{ST}','q_{TP}','q_{EP}','Location','east'); 
legend boxoff
lgd.NumColumns=4; %lgd.Position=poslgd;
xlabel('Time [h]');ylabel('Flows [m^3/h]')
xlim([0 tend]); ylim([-50; 800;])
set(gcf,'Position',pos);print('-depsc2','-r600','flows_qSdist')


figure(); % hot water tank volume
plot(tOL,Vh,'red')
hold on
h = plot([tOL(1);tOL(end)],[Vmax;Vmax],'k--');
h.Annotation.LegendInformation.IconDisplayStyle = 'off';
h = plot([tOL(1);tOL(end)],[Vmin;Vmin],'k--');
xlabel('Time [h]'); ylabel('Hot water volume [m^3]')
xlim([0 tend]); ylim([-100; 5100])
legend('V_h','Location','East');  legend boxoff
set(gcf,'Position',pos); print('-depsc2','-r600','V_qSdist')

figure(); % network pressure
plot(tOL,p,'Color',indigo)
legend('p'); legend boxoff
xlabel('Time [h]');ylabel('Network pressure [bar]')
xlim([0 tend]); ylim([4.5; 6;])
set(gcf,'Position',pos);print('-depsc2','-r600','p_qSdist')




%% Plot open loop. Step on qSP
load('ResultsOL_qSP')         % step in hot water supply. sent to consumers
clear qD qST
Vh=yOL(:,1); p=yOL(:,2);
qS=dOL(:,1); q=dOL(:,2);
qSP=uOL(:,1); qTP=uOL(:,2); qEP=uOL(:,3);
for i =1:length(Vh)
     if   Vh(i)  >= Vmax
          qD(i)  = qS(i)-qSP(i); % dumped
          qST(i) = 0;      % charge
     elseif Vh(i)< Vmax
          qD(i)  = 0;       % dump
          qST(i) = qS(i)-qSP(i);  % charge
     end
end

tOL=tOL*3600; tend=6;
figure(); % MVs
plot(tOL,qS,'Color',green)
hold on
plot(tOL,q,'Color',skyblue)                                      % demand q
plot(tOL,qSP,'--','Color',gold)                  % qSP from boiler to plant
plot(tOL,qD,'Color','k')                                          % qD dump
plot(tOL,qST,'--','Color',orange );                %qST from boiler to tank
plot(tOL,qTP,'-','Color',teal)                     % qTP from tank to plant
plot(tOL,qEP,'--','Color',purple)       % qEP from electric boiler to plant
lgd=legend('q_S','q','q_{SP}','q_D','q_{ST}','q_{TP}','q_{EP}','Location','east'); 
legend boxoff
lgd.NumColumns=4; %lgd.Position=poslgd;
xlabel('Time [s]');ylabel('Flows [m^3/h]')
xlim([0 tend]); ylim([-50; 800;])
set(gcf,'Position',pos);print('-depsc2','-r600','flows_qSPdist')

figure(); % hot water tank volume
plot(tOL,Vh,'red')
hold on
h = plot([tOL(1);tOL(end)],[Vmax;Vmax],'k--');
h.Annotation.LegendInformation.IconDisplayStyle = 'off';
h = plot([tOL(1);tOL(end)],[Vmin;Vmin],'k--');
xlabel('Time [s]'); ylabel('Hot water volume [m^3]')
xlim([0 tend]); ylim([-100; 5100])
legend('V_h','Location','NorthEast');  legend boxoff
set(gcf,'Position',pos); print('-depsc2','-r600','V_qSPdist')

figure(); % network pressure
plot(tOL,p,'Color',indigo)
legend('p','Location','East'); legend boxoff
xlabel('Time [s]');ylabel('Network pressure [bar]')
xlim([0 tend]); ylim([4.5; 6;])
set(gcf,'Position',pos);print('-depsc2','-r600','p_qSPdist')


%% Plot open loop. Step on q
load('ResultsOL_q')         % step in hot water supply. sent to consumers
clear qD qST
Vh=yOL(:,1); p=yOL(:,2);
qS=dOL(:,1); q=dOL(:,2);
qSP=uOL(:,1); qTP=uOL(:,2); qEP=uOL(:,3);
for i =1:length(Vh)
     if   Vh(i)  >= Vmax
          qD(i)  = qS(i)-qSP(i); % dumped
          qST(i) = 0;      % charge
     elseif Vh(i)< Vmax
          qD(i)  = 0;       % dump
          qST(i) = qS(i)-qSP(i);  % charge
     end
end

tOL=tOL*3600; tend=6;
figure(); % MVs
plot(tOL,qS,'Color',green)
hold on
plot(tOL,q,'Color',skyblue)                                      % demand q
plot(tOL,qSP,'--','Color',gold)                  % qSP from boiler to plant
plot(tOL,qD,'Color','k')                                          % qD dump
plot(tOL,qST,'--','Color',orange );                %qST from boiler to tank
plot(tOL,qTP,'-','Color',teal)                     % qTP from tank to plant
plot(tOL,qEP,'--','Color',purple)       % qEP from electric boiler to plant
lgd=legend('q_S','q','q_{SP}','q_D','q_{ST}','q_{TP}','q_{EP}','Location','east'); 
legend boxoff
lgd.NumColumns=4; %lgd.Position=poslgd;
xlabel('Time [s]');ylabel('Flows [m^3/h]')
xlim([0 tend]); ylim([-50; 800;])
set(gcf,'Position',pos);print('-depsc2','-r600','flows_qdist')

figure(); % hot water tank volume
plot(tOL,Vh,'red')
hold on
h = plot([tOL(1);tOL(end)],[Vmax;Vmax],'k--');
h.Annotation.LegendInformation.IconDisplayStyle = 'off';
h = plot([tOL(1);tOL(end)],[Vmin;Vmin],'k--');
xlabel('Time [s]'); ylabel('Hot water volume [m^3]')
xlim([0 tend]); ylim([-100; 5100])
legend('V_h','Location','NorthEast');  legend boxoff
set(gcf,'Position',pos); print('-depsc2','-r600','V_qdist')

figure(); % network pressure
plot(tOL,p,'Color',indigo)
legend('p','Location','East'); legend boxoff
xlabel('Time [s]');ylabel('Network pressure [bar]')
xlim([0 tend]); ylim([4.5; 6;])
set(gcf,'Position',pos);print('-depsc2','-r600','p_qdist')

%% Plot open loop. Step on qTP
load('ResultsOL_qTP')         % step in hot water supply. sent to consumers
clear qD qST
Vh=yOL(:,1); p=yOL(:,2);
qS=dOL(:,1); q=dOL(:,2);
qSP=uOL(:,1); qTP=uOL(:,2); qEP=uOL(:,3);
for i =1:length(Vh)
     if   Vh(i)  >= Vmax
          qD(i)  = qS(i)-qSP(i); % dumped
          qST(i) = 0;      % charge
     elseif Vh(i)< Vmax
          qD(i)  = 0;       % dump
          qST(i) = qS(i)-qSP(i);  % charge
     end
end

tOL=tOL*3600; tend=6;
figure(); % MVs
plot(tOL,qS,'Color',green)
hold on
plot(tOL,q,'Color',skyblue)                                      % demand q
plot(tOL,qSP,'--','Color',gold)                  % qSP from boiler to plant
plot(tOL,qD,'Color','k')                                          % qD dump
plot(tOL,qST,'--','Color',orange );                %qST from boiler to tank
plot(tOL,qTP,'-','Color',teal)                     % qTP from tank to plant
plot(tOL,qEP,'--','Color',purple)       % qEP from electric boiler to plant
lgd=legend('q_S','q','q_{SP}','q_D','q_{ST}','q_{TP}','q_{EP}','Location','east'); 
legend boxoff
lgd.NumColumns=4; %lgd.Position=poslgd;
xlabel('Time [s]');ylabel('Flows [m^3/h]')
xlim([0 tend]); ylim([-50; 800;])
set(gcf,'Position',pos);print('-depsc2','-r600','flows_qTPdist')

figure(); % hot water tank volume
plot(tOL,Vh,'red')
hold on
%h = plot([tOL(1);tOL(end)],[Vmax;Vmax],'k--');
%h.Annotation.LegendInformation.IconDisplayStyle = 'off';
%h = plot([tOL(1);tOL(end)],[Vmin;Vmin],'k--');
xlabel('Time [s]'); ylabel('Hot water volume [m^3]')
xlim([0 tend]); ylim([2499; 2501])
legend('V_h','Location','NorthEast');  legend boxoff
set(gcf,'Position',pos); print('-depsc2','-r600','V_qTPdist')

figure(); % network pressure
plot(tOL,p,'Color',indigo)
legend('p','Location','East'); legend boxoff
xlabel('Time [s]');ylabel('Network pressure [bar]')
xlim([0 tend]); ylim([4.5; 6;])
set(gcf,'Position',pos);print('-depsc2','-r600','p_qTPdist')