clear;

data = xlsread('Hourly_electricity_data_2016.xlsx','H2:N8785');

consumption = data(:, 7); % 7th column of data

% Power data _______________________________________________________________________________________________________________________________

figure(1);
area(data(:, 1:6), 'LineStyle', 'none');
colormap prism(6)
xlabel('Hours (h)');
ylabel('Power (MW)');
title('Hourly power data'); 
legend('Nuclear power', 'Hydropower', 'Wind power', 'CHP', 'Separate thermal production', 'Imports');
grid on;
pause;
print(figure(1), '-dpng', '-r400', 'demo_hourly_data');

% Duration curve ___________________________________________________________________________________________________________________________

sortedConsumption = sort(consumption, 'descend');

figure(2);
plot(sortedConsumption,'LineWidth', 1.5);
colormap jet
xlabel('Hours (h)');
ylabel('Power (MW)');
title('Power duration curve'); 
legend('Consumption');
axis([0 9000 0 15500]);
grid on;
pause;
print(figure(2), '-dpng', '-r400', 'demo_load_duration_curve');

% Peak and base loads ______________________________________________________________________________________________________________________

baseLoad = zeros(size(consumption));
peakLoad = zeros(size(consumption));
border   = 0.5*max(consumption); % Dividing the consumption at 50 % of the maximum

for j = 1:length(consumption)
   baseLoad(j) = min(border, sortedConsumption(j));
   peakLoad(j) = sortedConsumption(j)-baseLoad(j);
end

base_share = sum(baseLoad)/sum(sortedConsumption)*100;
peak_share = sum(peakLoad)/sum(sortedConsumption)*100;

figure(3);
area([baseLoad,peakLoad]);
colormap jet(2)
xlabel('Hours (h)');
ylabel('Power (MW)');
title('Power duration curve'); 
legend('Base load','Peak load');
text(1000,4000,strcat(num2str(base_share),' %'));
text(1000,9000,strcat(num2str(peak_share),' %'));
axis([0 9000 0 15500]);
grid on;
pause;
print(figure(3), '-dpng', '-r400', 'demo_base_peak');

% Variability histogram ____________________________________________________________________________________________________________________

timestep = 1; % 1 hour
ramps = zeros(length(consumption) - 1, 1);

for i = 1:length(ramps)
    ramps(i) = (consumption(i + 1) - consumption(i)) / timestep;     % delta P / delta t
end

minBin      = 100 * floor(min(min(ramps)) / 100);
maxBin      = 100 * ceil (max(max(ramps)) / 100);
rampBins    = minBin : 20 : maxBin;
frequencies = hist(ramps, rampBins);

figure(4);
plot(rampBins,frequencies,'LineWidth',1.5);
xlabel('Ramp-up/down (MW/h)');
ylabel('Frequency (number of 1 h timeslots)');
title('Variability analysis');
legend('Consumption');
grid on;
pause;
print(figure(4), '-dpng', '-r400', 'demo_variability');