optimising hybrid energy design using genetic algorithm

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.

Answer 1

[Alan Weiss]

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

[souleyman]

i hve made some chnges and when i run the program it just says bussy it's take too long and nothing has show either on command window or workspace what is the solution of that

[Alan Weiss]

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]

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

[e-mail]