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.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. 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. genetic algorithm MATLAB Answers — New Questions

​