Clear Filters
Clear Filters

optimising hybrid energy design using genetic algorithm

91 views (last 30 days)
i'm working on optimising a design of a hybrid PV/Wind energy system (with battery) using Genetic Algorithms ,and based on a research paper i have been able to code the following :
% Solar PV Generator script
function solargen = solar(Areapv)
% Input:
% Areapv = total area of all PV modules, one of the variables to be optimized
% x = solar irradiance (data read from an Excel file)
x = xlsread('éclairement_tlemcen.xlsx','A1:A8784');
effpv = 0.227; % efficiency of the PV module specified by manufacturer
% Output:
% solargen = power produced by the PV generator (kW)
% solargen = Areapv * effpv * Irrad * 0.001;
w = effpv * Areapv * 0.001;
solargen = w .* x;
end
wind turbine script
% Wind Generator script
function windgen = windpower(Areawt)
% Inputs:
% Areawt = total rotor swept area of all wind turbines used (m2)
% Areawt is one of the variables to be optimized
coeff_wt = 45/100; % coefficient of power of wind turbine
rho = 1.1839; % density of air (kg/m3) at 1 atm pres and 25 degC temp
vci = 2.01168; % wind turbine cut-in/min wind speed (m/s)
vrat = 13.8582; % wind turbine rated wind speed (m/s)
vco = 17.8816; % wind turbine cut-out/max wind speed (m/s)
% vact = actual wind speed (data read from an Excel file)
vact = xlsread('wind_speed_tlemcen.xlsx','A1:A8784');
y = length(vact); % number of wind speed data points
% windgen = power produced by the wind generator (kW)
b = zeros(1, y); % preallocation of memory outside loop
windgen = b.';
for a = 1:y
if ((vact(a) < vci)||(vact(a) > vco)) % wind speed less than cut-in or greater than cut-out
z = 0;
elseif ((vact(a) >= vci)&&(vact(a) < vrat)) % wind speed between cut-in and rated
z = 0.5 * coeff_wt * rho * Areawt * (vrat.^3) * ((vact(a).^3 - vci.^3)/(vrat.^3 - vci.^3)) * 0.001;
else % wind speed between rated and cut-out
z = 0.5 * coeff_wt * rho * Areawt * (vrat.^3) * 0.001;
end
windgen(a) = z;
end
end
battery script
% Battery System script
function [battpower,Pload] = battery(Areapv,Areawt,Capbatt)
% Inputs:
% Areapv = total area of all PV modules used (m2)
% Areawt = total rotor swept area of all wind turbines used (m2)
% Capbatt = total energy capacity of all batteries used (kWh)
% The three inputs are the variables to be optimized
% Pload = load demand (data read from an Excel file)
Pl = xlsread('charge_tlemcen.xlsx','A1:A8784');
Pload = 0.001 * Pl; % load power in kW
eff_inv = 88/100; % efficiency of inverter
eff_batt_cha = 92/100; % efficiency of battery charging
eff_batt_disch = 100/100; % efficiency of battery discharging
Ppv = solar(Areapv); % generation from solar
Pwind = windpower(Areawt); % generation from wind
soc_max = 98/100; % maximum state of charge of battery
dod_max = 80/100; % maximum depth of discharge
y = length(Pload);
% preallocation of memory outside loop for some variables
c = zeros(1, y);
Pgen = c.';
d = zeros(1, y);
soc = d.';
e = zeros(1, y);
dod = e.';
f = zeros(1, y);
Pdump = f.';
g = zeros(1, y);
battpower = g.';
h = zeros(1, y);
Pdef = h.';
soc(1) = soc_max; % initial state of charge of battery
dod(1) = 1 - soc(1); % initial state of discharge of battery
for b = 1:y
Pgen(b) = Ppv(b) + Pwind(b);
soc(b) = soc(b) + (battpower(b)/(1000*Capbatt));
if Pgen(b) > (Pload(b)/eff_inv) % generation > load
if soc(b) < soc_max % battery not charged fully
battpower(b) = (Pgen(b) - (Pload(b)/eff_inv)) *eff_batt_cha; % battery charges, battpower is +ve
else
soc(b) >= soc_max; % battery charged to maximum
battpower(b) = 0; % no more charging
Pdump(b) = Pgen(b) - (Pload(b)/eff_inv);
% surplus power is dumped after battery charges to maximum
end
elseif Pgen(b) < (Pload(b)/eff_inv) % generation < load
if (dod(b)) < dod_max % battery not dicharged to maximum
battpower(b) = -((Pload(b)/eff_inv) - Pgen(b)) * eff_batt_disch; % battery discharges, battpower is -ve
else
dod(b) >= dod_max; % battery discharged to maximum
battpower(b) = 0; % no more discharging
Pdef(b) = Pload(b) - (Pgen(b) + ((1000 * Capbatt)* ((soc(b)-(1-dod_max))))* eff_inv);
% deficit power persists after battery discharged fully
end
else % Pgen(b) = (Pload(b)/eff_inv) i.e. generation = load
battpower(b) = 0; % No charging or discharging
end
end
end
LPSP script which is the non linear constraint script
function [LPSP_value,LPSP_eq] = LPSP(Areapv,Areawt,Capbatt)
eff_inv = 88/100; % efficiency of inverter
s = solar(Areapv); % call solar function
w = windpower(Areawt); % call wind function
Pg = s + w; % generation = solar + wind
[b,Pl] = battery(Areapv,Areawt,Capbatt); % call battery function
Pgen = sum(Pg); % total generation
Pload = sum(Pl); % total load
battp = sum(b); % total battery power
LPS = sum(Pdef);
LPSP_value = (LPS / Pload) - 0.05; % LPSP <= 0.05
LPSP_eq = [];
end
the cost script which is the objective function
% Total System Cost script
% The total system cost is the sum of the following costs:
% Capital/investment costs
% Operation and maintenance costs
% Replacement costs
% Salvage revenue (negative cost)
% The present value of the above cost components are found for each of the
% three main system components: solar pv generator, wind generator and
% battery system; all costs are then added to get the total system cost
% The total system cost is the objective function to be optimized by a
% genetic algorithm. The constraints of this function are the loss of power
% supply probability (defined in a separate function) and the input
% variables.
function system_cost = cost(Areapv,Areawt,Capbatt)
%-------------------------------------------------------------------------%
% Constraints
Areapv_max = 20 * 1.63; % 20 PV modules maximum (32.6 m2 max area)
Areapv(Areapv>Areapv_max) = Areapv_max;
Areapv_min = 10 * 1.63; % 10 PV modules minimum (16.3 m2 min area)
Areapv(Areapv<Areapv_min) = Areapv_min;
Areawt_max = 10 * pi * (0.85344.^2); % 10 wind turbines maximum
% (22.8821 m2 max area)
Areawt(Areawt>Areawt_max) = Areawt_max;
Areawt_min = 3 * pi * (0.85344.^2); % 3 wind turbines minimum
% (6.8646 m2 min area)
Areawt(Areawt<Areawt_min) = Areawt_min;
Capbatt_max = 20 * 1.68; % 20 battery units maximum (33.6 kWh max)
Capbatt(Capbatt>Capbatt_max) = Capbatt_max;
Capbatt_min = 10 * 1.68; % 10 battery units minimum (16.8 kWh min)
Capbatt(Capbatt<Capbatt_min) = Capbatt_min;
%-------------------------------------------------------------------------%
% Project lifetime
%proj_life = 25; % years
%-------------------------------------------------------------------------%
% Rates applicable
%int = 5/100; % interest rate that affects all costs
%infl = 3/100; % inflation rate that affects salvage costs
%inc = 4/100; % non-inflation rate at which non-salvage costs increase
%-------------------------------------------------------------------------%
% Solar specifications
cap_pv = 300/1.63; % capital cost of PV module (184.0491 UK pounds/m2)
oandm_pv = 7.5/1.63; % o & m cost of PV module (4.6012 UK pounds/m2/yr)
sal_pv = 60/1.63; % salvage revenue of PV module (36.8098 UK pounds/m2)
%life_pv = 25; % lifetime of PV module (years)
%rep_pv = (proj_life / life_pv) - 1; % number of replacements in project (0)
%-------------------------------------------------------------------------%
% Wind specifications
cap_wt = 1125/(pi*(0.85344.^2)); % capital cost of turbine
% (491.6507 UK pounds/m2)
oandm_wt = 168.75/(pi*(0.85344.^2)); % o & m cost of turbine
% (73.7476 UK pounds/m2/yr)
sal_wt = 225/(pi*(0.85344.^2)); % salvage revenue of turbine
% (98.3301 UK pounds/m2)
%life_wt = 12.5; % lifetime of turbine (years)
%rep_wt = (proj_life / life_wt) - 1; % number of replacements in project (1)
%-------------------------------------------------------------------------%
% Battery specifications
cap_batt = 364/1.68; % capital cost of battery (216.6667 UK pounds/kWh)
oandm_batt = 3.64/1.68; % o & m cost of battery (2.1667 UK pounds/kWh/yr)
sal_batt = 36.4/1.68; % salvage revenue of battery (21.6667 UK pounds/kWh)
%life_batt = 2.5; % lifetime of battery (years)
%rep_batt = (proj_life / life_batt) - 1; % number of replacements in project;
% (9)
%-------------------------------------------------------------------------%
% Useful factors for net present value
%fac1 = (1+inc)/(1+int); % 0.9905
%fac2 = (1+infl)/(1+int); % 0.9810
%factor1a = symsum(fac1.^((k-1)*life_wt),k,1,rep_wt);
factor1a = 1; % summation of fac1^(k-1)*life_wt) for turbine replacements
%factor1b = symsum(fac1.^((k-1)*life_batt),k,1,rep_batt);
factor1b = 8.1943; % summation of fac1.^((k-1)*life_batt) for battery
% replacements
%factor2 = symsum(fac1.^k,k,1,proj_life);
factor2 = 22.1282; % summation of fac1^k for project life
% factor2 = fac1 + (fac1.^2) + (fac1.^3) + ... + (fac1.^25)
%factor3a = ((1+infl).^proj_life)/((1+int).^proj_life);
factor3a = 0.6183;
%factor3b = symsum(fac2.^(x*life_wt),x,1,rep_wt);
factor3b = 0.7863; % summation of fac2^(x*life_wt) for turbine life
%factor3c = symsum(fac2.^(x*life_batt),x,1,rep_batt);
factor3c = 7.1315; % summation of fac2^(x*life_batt) for battery life
%-------------------------------------------------------------------------%
% Capital costs and replacement costs
% Solar
pv_caprep = cap_pv * Areapv;
pv_caprep_npv = pv_caprep;
const_pv1 = pv_caprep_npv / Areapv;
% Wind
windg_caprep = cap_wt * Areawt;
windg_caprep_npv = windg_caprep * factor1a;
const_wt1 = windg_caprep_npv / Areawt;
% Battery
batt_caprep = cap_batt * Capbatt;
batt_caprep_npv = batt_caprep * factor1b;
const_batt1 = batt_caprep_npv / Capbatt;
%-------------------------------------------------------------------------%
% Operation and maintenance costs
% Solar
pv_oandm = oandm_pv * Areapv;
pv_oandm_npv = pv_oandm * factor2;
const_pv2 = pv_oandm_npv / Areapv;
% Wind
windg_oandm = oandm_wt * Areawt;
windg_oandm_npv = windg_oandm * factor2;
const_wt2 = windg_oandm_npv / Areawt;
% Battery
batt_oandm = oandm_batt * Capbatt;
batt_oandm_npv = batt_oandm * factor2;
const_batt2 = batt_oandm_npv / Capbatt;
%-------------------------------------------------------------------------%
% Salvage revenues
% Solar
pv_sal = sal_pv * Areapv;
pv_sal_npv = pv_sal * factor3a;
const_pv3 = pv_sal_npv / Areapv;
% Wind
windg_sal = sal_wt * Areawt;
windg_sal_npv = windg_sal * factor3b;
const_wt3 = windg_sal_npv / Areawt;
% Battery
batt_sal = sal_batt * Capbatt;
batt_sal_npv = batt_sal * factor3c;
const_batt3 = batt_sal_npv / Capbatt;
%-------------------------------------------------------------------------%
% Total system cost
% In general: system cost = capital cost + o&m cost - salvage revenue
system_cost = ((const_pv1 + const_pv2 - const_pv3) * Areapv) + ...
((const_wt1 + const_wt2 - const_wt3) * Areawt) +((const_batt1 + const_batt2 - const_batt3) * Capbatt);
end
the ga script
nvars = 3; % number of input variables
fun = @cost; % objective function to be optimized (system cost)
lb = [16.3 6.865 16.8]; % lower bounds of input variables
ub = [32.6 22.882 33.6]; % lower bounds of input variables
nonlcon = @LPSP; % nonlinear constraint function (loss of power supply
% probability
% Optimization command
[X,fval] = ga(fun,nvars,[],[],[],[],lb,ub,nonlcon);
i can provide the exel data files you need to run the function
these are the errors i could not fix :
Error using .*
Matrix dimensions must agree.
Error in solar (line 12)
solargen = w .* x;
Error in LPSP (line 9)
s = solar(Areapv); % call solar function
Error in createAnonymousFcn>@(x)fcn(x,FcnArgs{:}) (line 11)
fcn_handle = @(x) fcn(x,FcnArgs{:});
Error in constrValidate (line 23)
[cineq,ceq] = nonlcon(Iterate.x');
Error in gacommon (line 132)
[LinearConstr, Iterate,nineqcstr,neqcstr,ncstr] = constrValidate(NonconFcn, ...
Error in ga (line 336)
NonconFcn,options,Iterate,type] = gacommon(nvars,fun,Aineq,bineq,Aeq,beq,lb,ub, ...
Error in gen (line 28)
[X,fval] = ga(fun,nvars,[],[],[],[],lb,ub,nonlcon);
Caused by:
Failure in initial user-supplied nonlinear constraint function evaluation.

Accepted Answer

Alan Weiss
Alan Weiss on 3 May 2023
For an error of that type I think that you should learn to use the debugger. See Debug MATLAB Code Files and Set Breakpoints.
Also, you would probably get better performance by not reading Excel spreadsheets within a loop or function call. Instead, read the data once somewhere in your script, and pass the resulting data using the techniques in Parameterizing Functions.
Good luck,
Alan Weiss
MATLAB mathematical toolbox documentation
  3 Comments
Alan Weiss
Alan Weiss on 5 May 2023
Did you ask for a plot function or iterative display? If not, try running again with them enabled.
Alan Weiss
MATLAB mathematical toolbox documentation
Camilo Sosapanta
Camilo Sosapanta on 12 Feb 2024
Hello, souleyman I hope everything is fine for you and your family
could you please send me the codes and neccesary files
Thanks in advance
jcsosasalas@gmail.com

Sign in to comment.

More Answers (1)

Prashant Gaidhane
Prashant Gaidhane on 23 Apr 2024
hi give me excel sheets to run the code.. prashant.gaidhane@gcoej.ac.in

Categories

Find more on Wind Power in Help Center and File Exchange

Community Treasure Hunt

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

Start Hunting!