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

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.

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_PV=600000;
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 = 10;%%% changing

%% 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;
% load
% figure
% plot(Hours,Load)

% %% 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])
 %energy calculation
 energysolar=SolD*Day;
 % figure
 %plot(Day,energysolar)
 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=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,2000];
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,exitflag,output] = linprog(f,A,b,Aeq,beq,lb,ub);

Answer 1

[Jasvin]

Hi Hammad,

While I cannot explicitly solve your problem, I have some general tips when it comes to Optimization Problems with MATLAB and a possible bug in your code to list out.

The following lines should appear outside the for loop,

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

Currently this code is run on every single iteration of the loop, while this doesn't really change anything about the output of the code, I just wanted to point it out.

And when it comes to Optimization problems in MATLAB, it is often better to use the Problem-Based approach for solving the problem instead of using solver functions directly.

For example, in your use-case it becomes quite difficult to read through all the logic being used for your code because the MATLAB code and the optimization problem are so coupled together that it becomes difficult to separate out and focus on the actual problem being optimized. Whereas the Problem-Based approach produces code that is quite similar to the equations that you would write on a sheet of paper as part of the problem statement.

The following page gives a starter guide on how to use the Problem-Based approach for solving optimization problems:

https://www.mathworks.com/help/optim/ug/mixed-integer-linear-programming-basics-problem-based.html

If you have computed the parameters and input values correctly as per the provided code and your problem statement and you are still not getting any feasible points, then there really is no other option but to increase the size of the potential feasibility region by loosening the restrictions on the bounds of the decision variables like increasing the battery capacity's upper bound and so on.

I have another question as well, how is the initial guess for the decision variable being used to compute the constraints for optimizing the decision variable itself?

For example:

The "Electricity_discharged" variable is used in "Aeq" and the "Electricity_charged" variable is used in "Beq" to formulate the equality constraints for the problem. However these two variables were computed indirectly with the help of the "Capacity_Battery" variable which is initialized to 2000 at the start of the code and is itself a decision variable. From a code point of view all this works, but logically speaking from the optimization problem's point of view, it does not make sense to do this. Maybe a reformulation of the problem itself is in order?

Hope I was able to provide some insights and thinking points to take this further!

[Rizk Masoud]

any suggestion