Task1: Matlab code for Sizing of PV and BESS

% Clear the workspace to avoid conflicts
clear;
clear hours;

% Load data
load('data.mat');

% Extract data from the tables
pv_yield = table_PV(:); % PV power for a 1kWp panel in watts (W)
grainmill_demand = table_grainmill(:); % (W)
household_demand = table_household(:); % (W)
pump_demand = table_pump(:); % (W)

% Convert watts to kilowatts (kW)
pv_yield = pv_yield / 1000; % Convert W to kW for hourly period
grainmill_demand = grainmill_demand / 1000; % Convert W to kW
household_demand = household_demand / 1000; % Convert W to kW
pump_demand = pump_demand / 1000; % Convert W to kW

% Total demand is the sum of all individual demands
total_demand = grainmill_demand + household_demand + pump_demand;

% Create a datetime array for the time variable
start_time = datetime(2017,1,1,0,0,0);
num_hours = length(total_demand);
time = start_time + hours(0:num_hours-1);

% Create a timetable for easier handling
data = timetable(time', total_demand, pv_yield);

% Plot the overall data for a quick overview
figure;
subplot(2, 1, 1);
plot(data.Time, data.total_demand);
title('Electricity Demand');
xlabel('Time');
ylabel('Demand (kW)');

subplot(2, 1, 2);
plot(data.Time, data.pv_yield);
title('PV Power Yield');
xlabel('Time');
ylabel('PV Yield (kW)');

% Normalize the data for percentage plotting and finding the critical period
total_demand_percentage = (data.total_demand / max(data.total_demand)) * 100; % Percentage of maximum demand
pv_yield_percentage = (data.pv_yield / max(data.pv_yield)) * 100; % Percentage of maximum PV yield

% Define the period length (in days) for analysis
period_days = 10;
period_hours = period_days * 24;

% Calculate the difference between demand and PV generation in percentage
net_demand_percentage = total_demand_percentage - pv_yield_percentage;

% Initialize variables to store the critical period
max_net_demand_percentage = -inf;
critical_start_index = 1;

% Loop through the data to find the period with the highest net demand in percentage
for i = 1:(height(data) - period_hours + 1)
    current_net_demand_percentage = sum(net_demand_percentage(i:i + period_hours - 1));
    if current_net_demand_percentage > max_net_demand_percentage
        max_net_demand_percentage = current_net_demand_percentage;
        critical_start_index = i;
    end
end

% Extract the critical period data
critical_period_data = data(critical_start_index:(critical_start_index + period_hours - 1), :);
critical_period_demand_percentage = total_demand_percentage(critical_start_index:(critical_start_index + period_hours - 1));
critical_period_pv_yield_percentage = pv_yield_percentage(critical_start_index:(critical_start_index + period_hours - 1));



% Plot the critical period for visualization in percentage
figure;
plot(critical_period_data.Time, critical_period_demand_percentage, 'b', 'DisplayName', 'Demand (percentage)');
hold on;
plot(critical_period_data.Time, critical_period_pv_yield_percentage, 'r', 'DisplayName', 'PV Yield (percentage)');
title(['Electricity Demand and PV Power Yield (percentage) during the Critical Period (', num2str(period_days), ' days)']);
xlabel('Time');
ylabel('Percentage Value');
legend;
hold off;

% Use the start and end times from critical_period_data for sizing the PV and BESS
critical_start_time = critical_period_data.Time(1);
critical_end_time = critical_period_data.Time(end);

% Extract the data for the critical period
critical_period_data = data(data.Time >= critical_start_time & data.Time <= critical_end_time, :);

% Calculate system size for the critical period
pv_size_kWp = calculate_pv_size(critical_period_data);

% Calculate daily energy balance using scaled PV yield
daily_data = retime(critical_period_data, 'daily', 'sum');
daily_net_energy = daily_data.pv_yield * pv_size_kWp - daily_data.total_demand;

% Determine the maximum daily energy deficit
max_daily_deficit = -min(daily_net_energy); % Negative sign to get deficit

% Calculate BESS size based on max_daily_deficit
[bess_size, cost_pv, cost_bess] = calculate_bess_size(pv_size_kWp, max_daily_deficit);
fprintf('PV System Size: %.2f kWp\n', pv_size_kWp);
fprintf('BESS Size: %.2f kWh\n', bess_size);
fprintf('Cost of PV System: £%.2f\n', cost_pv);
fprintf('Cost of BESS: £%.2f\n', cost_bess);
total = cost_pv + cost_bess;
fprintf('Total cost: £%.2f\n', total);
% Function to calculate the PV size

function pv_size_kWp = calculate_pv_size(selected_data)
    % Calculate total demand and PV yield over the selected period
    total_demand = sum(selected_data.total_demand);
    total_pv_yield = sum(selected_data.pv_yield);

    % PV system size (kWp)
    pv_size_kWp = total_demand / total_pv_yield;
end

% Function to calculate the BESS size
function [bess_size_kWh, cost_pv, cost_bess] = calculate_bess_size(pv_size_kWp, max_daily_deficit)
    % BESS size (kWh)
    bess_size_kWh = max_daily_deficit;

    % Costs
    cost_pv = pv_size_kWp * 1000 * 2.73;
    cost_bess = bess_size_kWh * 86;
end