MATLAB Answers

Optimal Size of a Wind, Solar, and Battery to supply a load

9 views (last 30 days)
Joao Gomes
Joao Gomes on 22 Jun 2020 at 9:13
Hi everyone,
I’m developing a program that finds the perfect capacity size of a solar PV power plant, a wind farm, and an electrical battery. The objective is to supply the load demand of an isolated system.
First I started the code by calculating the solar and wind generation profiles, then I developed the battery storage procedure. In this part I don’t have any problem.
Now I want to find the cheapest configuration of the solar, wind and battery that will guarantee that the load demand is always supplied. I tried to use fmincon, linear prog and some other optimization functions without success.
By analyzing to my code do you guys have some suggestions on how to solve this problem?
%Load Management of Renewable Suppplied Microgrid
% Our Objective is to create a 100 % renewable microgrid in an island We want to find the minimal value of wind and solar to be installed to achieve 100% of the load %demand is supplied by renewable energy. To match the electricity supply and demand we are using an electrical battery.
Capacity_PV = 6; %% should be the decision variable
Capacity_Wind = 2; %% should be the decision variable
Capacity_Battery = 2000; %% should be the decision variable
%General Parameters:
T=24; %hours of the day
%% Solar Parameters
SystemLoss=0.14; % System losses from Solar PV Generation
%% Wind Parameters
Cp=0.3; %Reference Power Coefficient for Wind Speed
Air_D=1.15; %Air Density
Blade_R=15; %Radius of Turbine Blade
TurbinePower=300; %Initial power of the turbine - should be a variable later
%% Transmission Loss Parameters
%TransmissionWind=0.02; %Loss for Wind Generation Transmission until Battery
%TransmissionSolar=0.01; %Loss for Solar Generation Transmission until Battery
%% Operating Costs Parameters
Cost_Wind= 600000 ; %€ per 300 kW
Cost_Maintenance_Wind=0.02; %€ per kWh
Cost_Solar= 900 ; %€/kW
Cost_Maintenance_Solar=0.4; %€/kWh
Cost_Battery=2500; %€/kWh
%% Battery Parameters
SoCmin = 0.1; %Battery stored minimum value (%)
inv=0.98; %Inverter efficiency
dischargf=0.95; %Battery Discharching Efficiency
chargf=0.95; %Battery Charching efficiency
%% Ideas for Later - put a restriction of the solar and wind capacity based on the area available in the island
Available_Area=5100; %Available area in the island m2
PV_Occupation=0.054; %kW/m2
Wind_Occupation=78; %turbine/m2
Battery_Occupation = 0;
%% Reading of the excel files
[Load] = xlsread('Load2019','C:C');
[Temperature] = xlsread('IrradiationNew','D:D'); %Temperature Values
[Irradiation] = xlsread('IrradiationNew','B:B'); %Irradiation Values
[WindSpeed] = xlsread('Wind','H:H'); %wind speed (m2/s)
Hours=[1:1:8760];
Day=[1:1:24];
%% Calculation of the Solar Generation throught the year
Profile_PV = (1-SystemLoss)*Irradiation.*(1-0.005*(Temperature-25)); %Calculation of the Solar PV Generation kWh
Gen_PV=Capacity_PV * Profile_PV;
% %% Profile of Solar Generation throught the year
% figure
% plot(Hours,Profile_PV,'r') %Plot the PV Generation Profile
% title('Solar Generation in Corvo')
% xlabel('Hours')
% ylabel('kWh/kWp')
% xlim([1,8760])
% %% Average Solar Generation Day %%%%%%%
% Solarx=reshape(Profile_PV,24,[]);
% SolD = sum(Solarx,2)/365;
% figure
% plot(Day,SolD)
% title('Typical Solar Generation Day')
% xlabel('Hours')
% ylabel('kWh/kWp')
% xlim([1,24])
%% Calculation of the Wind Generation throught the year
RotorArea=pi*Blade_R^2; %Rotor area (m2)
Profile_Wind = (WindSpeed.^3*0.5*RotorArea*Cp*Air_D)/1000; %Calculation of the Wind Generation kWh
Gen_Wind=Capacity_Wind * Profile_Wind;
% %% Profile of Wind Generation throught the year
% figure
% plot(Hours,Profile_Wind,'r') %Plot the Wind Generation Profile
% title('Wind Generation in Corvo for E-30 Turbine')
% xlabel('Hours')
% ylabel('kWh')
% xlim([1,8760])
% %% Average Wind Generation Day %%%%%%%
% Windx=reshape(Profile_Wind,24,[]);
% WindD = sum(Windx,2)/365;
% figure
% plot(Day,WindD)
% title('Typical Wind Generation Day for E-30 Turbine')
% xlabel('Hours')
% ylabel('kW')
% xlim([1,24])
%% Battery Storage Calculation%%
for i=2:8760
SoCeop = 1; %% initial stored electricity, assume it is full
Electricity_discharged(1)=0;
Battery_Stored(1)=Capacity_Battery*SoCeop; %Initial electricity stored inside the battery in the beggining of the year
Battery_AvailabletoStore(i)=Gen_Wind(i)+Gen_PV(i)-Load(i); %Electricity that can be stored every hour
Battery_Stored(i)=Battery_Stored(i-1)+Battery_AvailabletoStore(i).*chargf; %Storage behavior
if Battery_Stored(i)/Capacity_Battery>1 %% this means it was fully charged
SoCeop = 1;
Battery_Stored(i)= Capacity_Battery * SoCeop; %%% here it means the battery is fully charged, so Battery_stored(i) = Capacity_Battery
Electricity_discharged(i)=0;
Electricity_charged(i) = (Battery_Stored(i)-Battery_Stored(i-1))/chargf;
Electricity_curtailed(i) = Battery_AvailabletoStore(i)-Electricity_charged(i);
elseif Battery_Stored(i)/Capacity_Battery<SoCmin %%%here it means that the battery discharged into the lower limit
SoCeop = SoCmin;
Battery_Stored(i)= Capacity_Battery * SoCeop; %%% Battery.Stored(i)=Capacity_Battery*SoCmin
Electricity_discharged(i)= (Battery_Stored(i-1)-Battery_Stored(i))*dischargf ; %%% Electricity_discharged(i) = (Battery_Stored(i-1)-Battery.Stored(i))*dischargf
Electricity_charged(i)=0;
Electricity_curtailed(i)=0;
else
Electricity_discharged(i) = max(-Battery_AvailabletoStore(i),0);
Electricity_charged(i) = max(Battery_AvailabletoStore(i),0);
Electricity_curtailed(i)=0;
end
end
%Energy_Lack(Capacity_Wind,Capacity_PV)=Load-Gen_Wind-Gen_PV
% %% Optimization part %%%%%%%%%%%%%%%%%%%
% % x = [Capacity_PV, Capacity_Wind, Capacity_Battery]; %% decision variables: capacity of each system
lb = [0,0,0]; ub = [1000,1000,1000];
f = [Cost_PV, Cost_Wind , Cost_Battery]; %% Cost per unit of supplied energy
A = [PV_Occupation, Wind_Occupation, Battery_Occupation]; %% Occupation area of each system
b = Available_Area; %%Island total area
Aeq = [Profile_PV, Profile_Wind,Electricity_discharged']; %% Aeq.x equals to the total energy output from the system
beq = Load + Electricity_curtailed' + Electricity_charged' ; %% total energy balance: Aeq.x = beq
[x,fval] = linprog(f,A,b,Aeq,beq,lb,ub);
Thank you for your time.
Best regards,
João

  0 Comments

Sign in to comment.

Answers (0)

Products


Release

R2019a