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

69 views (last 30 days)
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

Answers (2)

Bagas Bimantoro
Bagas Bimantoro on 7 Sep 2022
Hi mr. Joao
would you please share the paper or journal that you use as a reference sir
sincerely
Bagas

Hector Aguirre
Hector Aguirre on 26 Jul 2021
Hi, Joao.
I know it is a long time since you asked this, but any chance you solved the problem and are willing to share how?
Regardsm
Hector
  5 Comments
Johanes Kasilima
Johanes Kasilima on 20 May 2023
Hi Joao
Please send me through my email niwakasilima@gmail.com to see how you managed to solve it.
Rizk Masoud
Rizk Masoud on 2 Jan 2024
Hi Joao
Please send me through my email: rizk_masoud@yahoo.com to see how you managed to solve it.

Sign in to comment.

Products


Release

R2019a

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!