Matlab code for Sustainable Electrical Systems Coursework 2

                                                                                                                                                                                                                                                                                                                                                                                                 Sustainable Electrical Systems Coursework 2
Authors: Thiraphat Likitkumjorn
Date: 18 March 2024
Task 1
Let's now acquire the wind speed data first.
% Acquire the path of the current folder
currentFolder = pwd;
csvFileName = 'E4.50 Coursework_No2_2024_Data_F 2.xls';
csvFilePath = fullfile(currentFolder, csvFileName);

% read the csv file and acquire wind speed data
data = readtable(csvFilePath);
% Convert date to datetime
data.date = datetime(data.date, 'InputFormat', 'dd-MMM-yy');

% Plot Wind Output Profile
figure;
plot(data.date + hours(data.period/2), data.('RealWindOutput'), 'LineWidth', 1.5);
title('Wind Output Profile');
xlabel('Half hourly period');
ylabel('Real Wind Output (kWh)');
grid on;
% Update y-axis tick labels to show real values
ax = gca;
ax.YAxis.Exponent = 0;

% Plot SBP and SSP Profiles
figure;
plot(data.date + hours(data.period/2), data.SBP, 'LineWidth', 1.5, 'DisplayName', 'SBP (£/MWh)', 'Color', 'green');
hold on;
plot(data.date + hours(data.period/2), data.SSP, 'LineWidth', 1.5, 'DisplayName', 'SSP (£/MWh)', 'Color', 'red');
title('SBP and SSP Profiles');
xlabel('Half hourly period');
ylabel('Price (£/MWh)');
legend('Location', 'Best');
grid on;

hold off;
Task 2
% 'Real Wind Output' is the column representing the actual wind output
actual_wind_output = data.('RealWindOutput');

% Persistence forecast for 3.5 hours
forecast_horizon = 3.5 * 2; % 2 data points per hour
predicted_wind_output = [actual_wind_output(1:forecast_horizon); actual_wind_output(1:end-forecast_horizon)];

% Calculate forecast error
forecast_error = actual_wind_output - predicted_wind_output;
% Calculate Mean Absolute Error (MAE)
mae = mean(abs(forecast_error(9:721, 1)));

% Display MAE
disp(['Mean Absolute Error (MAE) for a 3.5-hour persistence forecast: ', num2str(mae)]);
% Persistence forecast for 1 hour
forecast_horizon2 = 1 * 2; % 2 data points per hour
predicted_wind_output2 = [actual_wind_output(1:forecast_horizon2); actual_wind_output(1:end-forecast_horizon2)];

% Calculate forecast error
forecast_error2 = actual_wind_output - predicted_wind_output2;
% Calculate Mean Absolute Error (MAE)
mae2 = mean(abs(forecast_error2(4:721, 1)));

% Display MAE
disp(['Mean Absolute Error (MAE) for a 1-hour persistence forecast: ', num2str(mae2)]);
% Plot the forecast error for 3.5 hours gate closure
figure;
plot(data.date + hours(data.period/2), forecast_error, 'LineWidth', 1.5, 'DisplayName', '3.5 Hours Gate Closure');
hold on;

% Plot the forecast error for 1 hour gate closure
plot(data.date + hours(data.period/2), forecast_error2, 'LineWidth', 1.5, 'DisplayName', '1 Hour Gate Closure');

title('Forecast Error using Persistence Forecast Technique');
xlabel('Date');
ylabel('Forecast Error (kWh)');
legend('Location', 'Best');
% Update y-axis tick labels to show real values
ax = gca;
ax.YAxis.Exponent = 0;
ytickValues = -14000:2000:14000;  % Adjust range as needed
yticks(ytickValues);
% Set the y-axis limits
ylim([-13000, 13000]);
grid on;
hold off;
Task 3 and 4
ssp = data.SSP;
sbp = data.SBP;

% Assume that the start calculation is at the 8th period for fair comparison between
% both gate closure
% Calculate Imbalance Income/Charges for a 3.5-hour persistence forecast

imbalance_income_charges = zeros(size(forecast_error) - 8); % Adjust size for starting from 9th period

reward = 0;
penalty = 0;

% Loop through forecast_error starting from 9th period
for i = 9:length(forecast_error)
    % Check for NaN in forecast_error
    if ~isnan(forecast_error(i))
        if forecast_error(i) >= 0
            imbalance_income_charges(i-8) = forecast_error(i) * ssp(i) / 1000;
            % Calculate Energy Sales Income
            reward = reward + imbalance_income_charges(i-8);
        else
            imbalance_income_charges(i-8) = forecast_error(i) * sbp(i) / 1000;
            % Calculate Total Penalties
            penalty = penalty - imbalance_income_charges(i-8);
        end
    end
end

% Calculate Imbalance Income/Charges for a 1-hour persistence forecast

imbalance_income_charges2 = zeros(size(forecast_error2) - 8); % Adjust size for starting from 9th period

reward2 = 0;
penalty2 = 0;

% Loop through forecast_error2 starting from 9th period
for i = 9:length(forecast_error2)
    % Check for NaN in forecast_error2
    if ~isnan(forecast_error2(i))
        if forecast_error2(i) >= 0
            imbalance_income_charges2(i-8) = forecast_error2(i) * ssp(i) / 1000;
            % Calculate Energy Sales Income
            reward2 = reward2 + imbalance_income_charges2(i-8);
        else
            imbalance_income_charges2(i-8) = forecast_error2(i) * sbp(i) / 1000;
            % Calculate Total Penalties
            penalty2 = penalty2 - imbalance_income_charges2(i-8);
        end
    end
end

% Plot the results
figure;
plot(data.date(9:end) + hours(data.period(9:end)/2), imbalance_income_charges, 'LineWidth', 1.5);
hold on;
plot(data.date(9:end) + hours(data.period(9:end)/2), imbalance_income_charges2, 'LineWidth', 1.5);
title('Half-Hourly Imbalance Income/Charges');
xlabel('Date');
ylabel('Imbalance Income/Charges (£)')
legend('3.5-hours Gate Closure', '1-hour Gate Closure');
ytickValues = -800:100:400;  % Adjust range as needed
yticks(ytickValues);
ylim([-800, 250]);
grid on;
hold off;

% Calculate total penalties for a 3.5-hour persistence forecast
total_penalties = penaltie;
% Calculate Energy Sales Income
% Assume that the market price is the average of SBP and SSP
market_price= (ssp+sbp)./2;
energy_sales_income = sum(predicted_wind_output(9:721, 1) .* market_price(9:721, 1)./1000)+reward;
% Calculate Overall Income for a 3.5-hour persistence forecast
overall_income = energy_sales_income - total_penalties;
fprintf('Energy Sales Income (3.5-hours persistence forecast): £%.2f\n', energy_sales_income);
fprintf('Total Penalties (3.5-hours persistence forecast): £%.2f\n', total_penalties);
fprintf('Overall Income (3.5-hours persistence forecast): £%.2f\n', overall_income);
% Calculate total penalties for a 1-hour persistence forecast
total_penalties2 = penaltie2;
% Calculate Energy Sales Income for a 1-hour persistence forecast
energy_sales_income2 = sum(predicted_wind_output2(9:721, 1) .* market_price(9:721, 1)./1000)+reward2;
% Calculate Overall Income for a 1-hour persistence forecast
overall_income2 = energy_sales_income2 - total_penalties2;
fprintf('Energy Sales Income (1-hour persistence forecast): £%.2f\n', energy_sales_income2);
fprintf('Total Penalties (1-hour persistence forecast): £%.2f\n', total_penalties2);
fprintf('Overall Income (1-hour persistence forecast): £%.2f\n', overall_income2);

%Output for task4:
%Energy Sales Income (3.5-hours persistence forecast): £199661.28
%Total Penalties (3.5-hours persistence forecast): £23408.66
%Overall Income (3.5-hours persistence forecast): £176252.62
%Energy Sales Income (1-hour persistence forecast): £191066.94
%Total Penalties (1-hour persistence forecast): £9222.81
%Overall Income (1-hour persistence forecast): £181844.13