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;
bluecg  = [44,144,195]/256; 
movcg   = [198,102,177]/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 400 300]; % figure position on screen
poslgd  = [0.5263    0.6032    0.3770    0.1963];[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    = 35;   % [h] plot end time
lw      = 2; % thickness limits
%% Plot SRC
load('ResultsSRC')
figure(1); % supply. same for SRC and 3C
plot(tSRC,qS,'Color',green)
xlabel('Time [h]');ylabel({'Hot water from';'waste heat boiler [m^3/h]'})
xlim([0 tend]); ylim([200; 1100;])
legend('qS');legend boxoff
set(gcf,'Position',pos);print('-depsc2','-r600','qS.eps')

figure(2); % hot water demand
plot(tSRC,qref,'Color',skyblue)
xlabel('Time [h]');ylabel('Hot water demand [m^3/h]')
xlim([0 tend]); ylim([200; 1100;])
legend('q','Location','west');legend boxoff
set(gcf,'Position',pos); print('-depsc2','-r600','q.eps')

figure(3); % hot water tank volume
plot(tSRC,ySRC(:,1),'red')
hold on
h = plot([tSRC(1);tSRC(end)],[Vmax;Vmax],'k--','LineWidth',lw);
h.Annotation.LegendInformation.IconDisplayStyle = 'off';
h = plot([tSRC(1);tSRC(end)],[Vmin;Vmin],'k--','Linewidth',lw);
xlabel('Time [h]'); ylabel('Hot water volume [m^3]')
xlim([0 tend]); ylim([-200; 5200])
legend('V_h','Location','East');  legend boxoff
set(gcf,'Position',pos); print('-depsc2','-r600','VSRC.eps')

figure(4); % network pressure
plot(tSRC,ySRC(:,2),'Color',indigo)
hold on
plot(tSRC,pref,'--','Color',greencg)
legend('p','p^s','Location','east'); legend boxoff
xlabel('Time [h]');ylabel('Network pressure [bar]')
xlim([0 tend]); ylim([0; 7.5;])
set(gcf,'Position',pos);print('-depsc2','-r600','pSRC.eps')

figure(5); %MVs: qSP, qD, qTP, qEB
plot(tSRC,qSPSRC,'Color',gold)
hold on
plot(tSRC,qDSRC,'Color','black')
plot(tSRC,qTPSRC,'Color',teal)
plot(tSRC,qEPSRC,'Color',purple)
h = plot([tSRC(1);tSRC(end)],[qmax;qmax],'k--','Linewidth',lw);
h.Annotation.LegendInformation.IconDisplayStyle = 'off';
h = plot([tSRC(1);tSRC(end)],[qmin;qmin],'k--','Linewidth',lw);
h.Annotation.LegendInformation.IconDisplayStyle = 'off';
lgd=legend('q_{SP}','q_D','q_{TP}','q_{EP}','Location','east'); legend boxoff
lgd.NumColumns=2; lgd.Position=poslgd;
xlabel('Time [h]');ylabel('Flows (MVs) [m^3/h]')
xlim([0 tend]); ylim([-100; 1100;])
set(gcf,'Position',pos);print('-depsc2','-r600','MVSRC.eps')

%% Plot parallel control
filename = {'Results3C_05','Results3C_025'};
figname  = {'05','025'};

for i =1
load(filename{i})
figure(10*i);                                       % hot water tank volume
plot(t3C,y3C(:,1),'red')
hold on
h = plot([t3C(1);t3C(end)],[Vmax;Vmax],'k--','LineWidth',lw);
h.Annotation.LegendInformation.IconDisplayStyle = 'off';
h = plot([t3C(1);t3C(end)],[Vmin;Vmin],'k--','LineWidth',lw);
xlabel('Time [h]'); ylabel('Hot water volume [m^3]')
xlim([0 tend]); ylim([-100; 5100;])
legend('V_h','Location','East'); legend boxoff
fgn = strcat('V3C',figname{i});
set(gcf,'Position',pos); print('-depsc2','-r600',fgn)

figure(20*i);                                            % network pressure
plot(t3C,y3C(:,2),'Color',indigo)
hold on
plot(tSRC,pref,'--','Color',greencg)
plot(t3C,pTPref,'--','Color',cerulean)
plot(t3C,pEPref,'--','Color',orange)
lgd=legend('p','p_{SP}^s','p_{TP}^s','p_{EP}^s'); legend boxoff
lgd.NumColumns=2; lgd.Position=[0.1707    0.1796    0.3770    0.2283];
xlabel('Time [h]');ylabel('Network pressure [bar]')
xlim([0 tend]); ylim([0; 7.5])
fgn = strcat('p3C',figname{i});
set(gcf,'Position',pos);print('-depsc2','-r600',fgn)

figure(30*i);                                       %MVs: qSP, qD, qTP, qEB
plot(t3C,qSP3C,'Color',gold)
hold on
plot(t3C,qD3C,'Color','black')
plot(t3C,qTP3C,'Color',teal)
plot(t3C,qEP3C,'Color',purple)
h = plot([t3C(1);t3C(end)],[qmax;qmax],'k--','LineWidth',lw);
h.Annotation.LegendInformation.IconDisplayStyle = 'off';
h = plot([t3C(1);t3C(end)],[qmin;qmin],'k--','LineWidth',lw);
h.Annotation.LegendInformation.IconDisplayStyle = 'off';
lgd = legend('q_{SP}','q_D','q_{TP}','q_{EP}','Location','northeast'); legend boxoff
lgd.NumColumns=2; lgd.Position=poslgd;
xlabel('Time [h]');ylabel('FLows (MVs) [m^3/h]')
xlim([0 tend]); ylim([-100; 1100;])
set(gcf,'Position',pos)
fgn = strcat('MV3C',figname{i});
print('-depsc2','-r600',fgn)
end
%% Plot MPC

load ResultsMPC_Ph10

tgrid = 0:st:time; % each sampling time

figure(300); % hot water tank volume
plot(tgrid,[2500;x_vec(:,1)],'red')
hold on
h = plot([tgrid(1);tgrid(end)],[Vmax;Vmax],'k--','Linewidth',lw);
h.Annotation.LegendInformation.IconDisplayStyle = 'off';
h = plot([tgrid(1);tgrid(end)],[Vmin;Vmin],'k--','Linewidth',lw);
xlabel('Time [h]'); ylabel('Hot water volume [m^3]')
xlim([0 tend]); ylim([-200; 5200;])
legend('V_h','Location','East'); legend boxoff
set(gcf,'Position',pos); print('-depsc2','-r600','VMPC.eps')

figure(400); % network pressure
plot(tgrid,[5;x_vec(:,2)],'Color',indigo)
hold on
plot(tSRC,pref,'--','Color',greencg)
legend('p','p_{SP}^s','Location','east'); legend boxoff
xlabel('Time [h]');ylabel('Network pressure [bar]')
xlim([0 tend]); ylim([0; 7.5;])
set(gcf,'Position',pos);print('-depsc2','-r600','pMPC.eps')

figure(500); %MVs: qSP, qD, qTP, qEB
plot(tgrid,[500;u_vec(:,1)],'Color',gold)
hold on
plot(tgrid,[0;u_vec(:,2)],'Color','black')
plot(tgrid,[0;u_vec(:,3)],'Color',teal)
plot(tgrid,[0;u_vec(:,4)],'Color',purple)
h = plot([t3C(1);t3C(end)],[qmax;qmax],'k--');
h.Annotation.LegendInformation.IconDisplayStyle = 'off';
h = plot([t3C(1);t3C(end)],[qmin;qmin],'k--');
h.Annotation.LegendInformation.IconDisplayStyle = 'off';
lgd = legend('q_{SP}','q_D','q_{TP}','q_{EP}','Location','northeast'); legend boxoff
lgd.NumColumns=2; lgd.Position=poslgd;
xlabel('Time [h]');ylabel('Flows (MVs) [m^3/h]')
xlim([0 tend]); ylim([-100; 1100;])
set(gcf,'Position',pos)
print('-depsc2','-r600','MVMPC.eps')

%% Compare network pressure for SRC and MPC. short time
figure(); % network pressure
plot(tSRC,ySRC(:,2),'Color',bluecg)
hold on
plot(tgrid,[5;x_vec(:,2)],'Color',indigo)
plot(tSRC,pref,'--','Color',greencg)
legend('p^{SRC}','p^{MPC}','p^s','Location','east'); legend boxoff
xlabel('Time [h]');ylabel('Network pressure [bar]')
xlim([11.99 12.1]); %ylim([3.4; 5.2;])
set(gcf,'Position',pos);print('-depsc2','-r600','pSRCvsMPC.eps')

%% Compare qEP for SRC, 3C and MPC 
figure(); % teal, purple, indigo
plot(tSRC,qEPSRC,'Color',indigo)
hold on
plot(t3C,qEP3C,'Color',greencg)
plot(tgrid,[0;u_vec(:,4)],'Color',bluecg)
legend('SRC','3C','MPC','Location','east'); legend boxoff
xlabel('Time [h]');ylabel('q_{EP} [m^3/h]')
xlim([24 26]); ylim([-100; 600;])
set(gcf,'Position',pos);print('-depsc2','-r600','qEPSRC-3C-MPC.eps')

%% Compare qTP for SRC, 3C and MPC 
figure(); % teal, purple, indigo
plot(tSRC,qTPSRC,'Color',indigo)
hold on
plot(t3C,qTP3C,'Color',greencg)
plot(tgrid,[0;u_vec(:,3)],'Color',bluecg)
legend('SRC','3C','MPC','Location','east'); legend boxoff
xlabel('Time [h]');ylabel('q_{TP} [m^3/h]')
xlim([24 26]); ylim([-100; 600;])
set(gcf,'Position',pos);print('-depsc2','-r600','qTPSRC-3C-MPC.eps')
