clear;

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

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

% Power duration curve

sorted_consumption = sort(consumption,'descend');

% Calculating the share of peak and base load

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

for j = 1:length(consumption)
   base_load(j) = min(border, sorted_consumption(j));
   peak_load(j) = sorted_consumption(j)-base_load(j);
end

base_share = sum(base_load)/sum(sorted_consumption)*100;
peak_share = sum(peak_load)/sum(sorted_consumption)*100;
  
% Ramp-up/down

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

rbins = 100*floor(min(ramps)/100):20:100*ceil(max(ramps)/100); % Ramp-up bins
freq = hist(ramps,rbins); % Frequencies

% Plots

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;

figure(2);
plot(sorted_consumption,'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;

figure(3);
area([base_load,peak_load]);
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;

figure(4);
plot(rbins,freq,'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(1), '-dpng', '-r400', 'demo_hourly_data');
print(figure(2), '-dpng', '-r400', 'demo_load_duration_curve');
print(figure(3), '-dpng', '-r400', 'demo_base_peak');
print(figure(4), '-dpng', '-r400', 'demo_variability');