MATLAB Answers

Use of Genetic Algorithm feature for optimization of technical design

28 views (last 30 days)
Nathan-Douglas Ngumi
Nathan-Douglas Ngumi on 8 Sep 2021
Greetings MATLAB Community.
I recently came across some academic publications about using metaheuristic optimization techniques like genetic algorithms and particle swarm optimization to optimize design of hybrid renewable energy systems.
I tried to implement one using a genetic algorithm for design of a solar-wind-battery renewable energy system.
I wrote separate scripts for the solar component, wind component, battery component, loss of power supply probability (constraint function), cost (objective function) and the genetic algorithm. The objective of the optimization is to minimize the system cost, constrained by the desired loss of power supply probability limit.
The issue arises when the algorithm script calls the constraint function script.
I will paste the code to illustrate.
*******************************************************************************************************************************************************************
% 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 = readcell('data.xlsx','Sheet','WeatherData','Range','B3:B26');
y = cell2mat(x);
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 .* y;
end
*******************************************************************************************************************************************************************
% 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)
x = readcell('data.xlsx','Sheet','WeatherData','Range','C3:C26');
vact = cell2mat(x);
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 * (vact(a).^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 * (vact(a).^3) * 0.001;
end
windgen(a) = z;
end
end
*******************************************************************************************************************************************************************
% 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)
x = readcell('data.xlsx','Sheet','LoadData','Range','B2:B25');
Pl = cell2mat(x); % load power in W
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
*******************************************************************************************************************************************************************
% Loss of Power Supply Probability (LPSP) script
% This is the nonlinear constraint of the optimization.
% System designers typically set a value of LPSP that should not be exceeded.
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
LPSP_value = ((Pload - ((Pgen + battp - battp) * eff_inv)) / Pload) - 0.05; % LPSP <= 0.05
LPSP_eq = [];
end
*******************************************************************************************************************************************************************
% 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
*******************************************************************************************************************************************************************
% Genetic algorithm script for optimization of cost function
% MATLAB inbuilt function ga is used for genetic algorithm optimization
% Syntaxes of ga function:
% x = ga(fun,nvars)
% x = ga(fun,nvars,A,b)
% x = ga(fun,nvars,A,b,Aeq,beq)
% x = ga(fun,nvars,A,b,Aeq,beq,lb,ub)
% x = ga(fun,nvars,A,b,Aeq,beq,lb,ub,nonlcon)
% x = ga(fun,nvars,A,b,Aeq,beq,lb,ub,nonlcon,options)
% x = ga(fun,nvars,A,b,[],[],lb,ub,nonlcon,IntCol)
% x = ga(fun,nvars,A,b,[],[],lb,ub,nonlcon,IntCol,options)
% x = ga(problem)
% [x,fval] = ga(___)
% [x,fval,exitflag,output] = ga(___)
% [x,fval,exitflag,output,population,scores] = ga(___)
% Any element in ga(___) not used is replaced with [] as shown above
% Parameters of optimization:
nvars = 3; % number of input variables
fun = @cost; % objective function to be optimized (system cost)
nonlcon = @LPSP; % nonlinear constraint function (loss of power supply
% probability)
lb = [16.3 6.865 16.8]; % lower bounds of input variables
ub = [32.6 22.882 33.6]; % lower bounds of input variables
% Optimization command
[x,fval] = ga(fun,nvars,[],[],[],[],lb,ub,nonlcon);
*******************************************************************************************************************************************************************
When the optimization command is executed, this is what is displayed in the error window:
-----------------------------------------------------------------------------------------------------------------------------
>> [x,fval] = ga(fun,nvars,[],[],[],[],lb,ub,nonlcon)
Not enough input arguments.
Error in LPSP (line 10)
w = windpower(Areawt); % call wind 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 125)
[LinearConstr, Iterate,nineqcstr,neqcstr,ncstr] = constrValidate(NonconFcn, ...
Error in ga (line 363)
NonconFcn,options,Iterate,type] = gacommon(nvars,fun,Aineq,bineq,Aeq,beq,lb,ub, ...
Caused by:
Failure in initial user-supplied nonlinear constraint function evaluation.
-------------------------------------------------------------------------------------------------------------------------
The error points to the windpower function, yet this function when executed separately at the command prompt with an input is working okay.
What could be the issue? I am at my wits end on this. Please shed some light.
Thanks.

Answers (1)

Steven Lord
Steven Lord on 8 Sep 2021
What signature does the documentation page for the ga function state your nonlinear constraint function must have? [How many input arguments does it accept and how many output arguments does it return?]
Does the signature of your nonlinear constraint function match that signature?
The techniques described on the "Passing Extra Parameters" documentation page linked in the first Note block in the Description section on that page may be of use to you.
  1 Comment
Nathan-Douglas Ngumi
Nathan-Douglas Ngumi on 9 Sep 2021
Sir, thanks, the documentation page helped.
The issue was that I defined three independent variables (Areapv, Areawt, Capbatt) in the objective and constraint functions instead of one (x), specifying x(1) in place of Areapv, x(2) in place of Areawt and x(1) in place of Capbatt. The MATLAB solvers accept only one independent variable.
Thanks.

Sign in to comment.

Community Treasure Hunt

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

Start Hunting!