clear all

close all
clc
%% Problem definition
%hade bala
C=0.9;%price of electricity
W=0.2;%lose of load probability%
K=0.99;%renewable energy factor%
nvars=1;%only one system
LB=[0 0 1 0]; % Lower bound of problem
UB=[45 8 20 10]; % Upper bound of problem
%%
max_it=20;%100;% Maximum number of iterations
NPOP=5;%20;% Number of population
% Determine the maximum step for velocity
for d=1:4
if LB(d)>-1e20 && UB(d)<1e20
velmax=(UB(d)-LB(d))/NPOP;
else
velmax=inf;
end
end
%% PSO initial parameters
% w=0.5; % Inertia weight
% wdamp=0.99; % Inertia weight damping ratio
% c1=2; % Personal learning coefficient
% c2=2; % Global learning coefficient
phi1=2.05;
phi2=2.05;
phi=phi1+phi2;
chi=2/(phi-2+sqrt(phi^2-4*phi));
w1=chi; % Inertia weight
%wdamp=1; % Inertia weight damping ratio
c1=chi*phi1; % Personal learning coefficient
c2=chi*phi2; % Global learning coefficient
%% Initialization
tic
empty_particle.position=[];
empty_particle.velocity=[];
empty_particle.cost=[];
empty_particle.best.position=[];
empty_particle.best.cost=[];
%---------
particle=repmat(empty_particle,NPOP,1);
globalbest.cost=inf;
globalbest.position=[];
%rnwfct_best=inf;
for i=1:NPOP
cc=1;%a value for cost
ww=0.3;% a value for lose of load probability%
kkk=2;%renewable energy factor%
ff=0;
%C=0.9;%price of electricity
%W=0.2;%lose of load probability%
%K=0.99;%renewable energy factor%
while ww>=W | kkk>=K% LOLP/REfactor
particle(i).position(1,:)=unifrnd(0,45,1,nvars);%pv kW
particle(i).position(2,:)=unifrnd(0,8,1,nvars);%autonomy days
particle(i).position(3,:)=unifrnd(1,20,1,nvars);%number of houses
particle(i).position(4,:)=unifrnd(0,10,1,nvars);%number of wind turbine
for g=1:4
particle(i).velocity(g,:)=rand(1,nvars);
end
%----convert------------
 p_npv=particle(i).position(1);
ad=particle(i).position(2);
houses=round(particle(i).position(3));
nwt=round(particle(i).position(4));
%-----------------------

[LPSP,price_electricity,renewable_factor,b,ali]=techno_economic_analysis_pso(houses,p_npv,ad,nwt);
ff=ff+1;
ww(i)=LPSP;
kkk(i)=renewable_factor;
end

particle(i).cost=price_electricity;
particle(i).best.position=particle(i).position;
particle(i).best.cost=particle(i).cost;
if particle(i).best.cost<globalbest.cost
globalbest=particle(i).best;

end
end
Fminn=zeros(max_it,1);
%% PSO main loop
% disp('Iteration Reliability');
% disp('-----------------------------');
for u=1:max_it
vv=0;
for i=1:NPOP
cc=1;%a value for cost
LPSP=0.3;% a value for lose of load probability%
renewable_factor=2;%
bb=0;
%C=0.9;%price of electricity
%W=0.2;%lose of load probability%
%K=0.99;%renewable energy factor%
% ww=0.3*100;% a value for lose of load probability%
% kkk=2*100;%renewable energy factor%
while (LPSP)>=W | (renewable_factor)>=K

for y=1:4
particle(i).velocity(y,:)=w1*particle(i).velocity(y,:)+c1*rand*...
(particle(i).best.position(y,:)-particle(i).position(y,:))...
+c2*rand*(globalbest.position(y,:)-particle(i).position(y,:));

%particle(i).velocity(y,:)=min(max(particle(i).velocity(y,:),-velmax),velmax);

particle(i).position(y,:)=particle(i).position(y,:)+particle(i).velocity(y,:);

% flag=(particle(i).position(kk,:)<LB(kk) | particle(i).position(kk,:)>UB(kk));
% particle(i).velocity(flag)=-particle(i).velocity(flag);
particle(i).position(y,:)=min(max(particle(i).position(y,:),LB(y)),UB(y));

end
%p_npv(i,:)=round(particle(i).position(1,:));
oo=0;
%ad(i,:)=round(particle(i).position(2,:));
p_npv=round(particle(i).position(1));
ad=round(particle(i).position(2));
houses=round(particle(i).position(3));
nwt=round(particle(i).position(4));

%[cc ww]=constraints(c,w,nnn(i,:),zz(i,:),zmax,nmax,N);
[LPSP,price_electricity,renewable_factor,b,ali]=techno_economic_analysis_pso(houses,p_npv,ad,nwt);
bb=bb+1;

end
%----convert------------
% az(i,:)=round(particle(i).position(1,:));
% z(i,:)=round(particle(i).position(2,:));
% n(i,:)=round(particle(i).position(3,:));
% sol(i).pos(1,:)=round(particle(i).position(1,:));
% sol(i).pos(2,:)=round(particle(i).position(2,:));
% sol(i).pos(3,:)=round(particle(i).position(3,:));

%-----------------------
%[LPSP,price_electricity,renewable_factor]=techno_economic_analysis_pso(houses,p_npv,ad,nwt);
particle(i).cost=price_electricity;
%rnwfct=renewable_factor;
%particle(i).cost=objective_function(sol(i).pos);
vv=vv+1;
if particle(i).cost<particle(i).best.cost
particle(i).best.cost=particle(i).cost;
particle(i).best.position=particle(i).position;
if particle(i).best.cost<globalbest.cost%& rnwfct<rnwfct_best
globalbest=particle(i).best;
% rnwfct_best=rnwfct;
end
end

end

Fminn(u)=globalbest.cost;
Xmin=globalbest.position;
p_npv=round(globalbest.position(1));
ad=round(globalbest.position(2));
houses=round(globalbest.position(3));
nwt=round(globalbest.position(4));

% w=wdamp*w;
% disp(['Iteration ',num2str(u),': Best cost= ',num2str(Fminn(u))]);
end
Time=toc;
%% Results
close all;
% figure;
% plot(Fminn,'Linewidth',3);
% ylabel('price of electricity');
% xlabel('Number of iterations');
% ylim([0 0.4])
% xlim([1 max_it ])
Fmin=min(Fminn);
Xmin;
[LPSP,price_electricity,renewable_factor,b,ali]=techno_economic_analysis_pso(houses,p_npv,ad,nwt);
format long
result=[p_npv,ad,houses,nwt,LPSP,price_electricity,renewable_factor,b(1),b(2),b(3),b(4),b(5)];
% bar(result)
% legend('pv(kW)','days of autonomy','number of houses','number of wind
turbines','PV','WIND','BATTERY','DIESEL','loss of load probability','price of
electricity($/kW)','renewable factor')
figure,subplot(2,1,1),x=1:168;bar(x,ali(:,1:4),'group');colormap
jet,subplot(2,1,2),x=1:168;bar(x,ali(:,5:6),'group');colormap jet
p_npv
ad
houses
nwt
LPSP
price_electricity
renewable_factor

function
[LPSP,price_electricity,renewable_factor,b,ali,Edump]=techno_economic_analysis_pso(houses,p_npv,ad,nwt,
nPng)
%@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
% %% 3)load inputs RAFSANJAN/KERMAN
% load('rafsanjan.mat');
% load('WindTurbines.mat');
% wind_speed=rafsanjan(:,1);[wind_speed]=min10tohourly(wind_speed);%hourly
% temperature=rafsanjan(:,2);[temperature]=min10tohourly(temperature);%hourly
% solar_radiation=rafsanjan(:,3)/1000;[solar_radiation]=min10tohourly(solar_radiation);%hourly%kw
% clear rafsanjan
%@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
% %% 2)load inputs KHASH/SISTAN BALOCHESTAN
% load('sistan_khash.mat');
% load('WindTurbines.mat');
% wind_speed=sistan_khash(:,1);[wind_speed]=min10tohourly(wind_speed);%hourly
% temperature=sistan_khash(:,2);[temperature]=min10tohourly(temperature);%hourly
% solar_radiation=sistan_khash(:,3)/1000;[solar_radiation]=min10tohourly(solar_radiation);%hourly%kw
% clear sistan_khash
%@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
%% 1)load inputs NAHAVAND
load('nahavand.mat');nahavand = circshift(nahavand,60);
load('WindTurbines.mat');
wind_speed=nahavand(:,1);[wind_speed]=min10tohourly(wind_speed);%hourly
temperature=nahavand(:,2);[temperature]=min10tohourly(temperature);%hourly
solar_radiation=nahavand(:,3)/1000;[solar_radiation]=min10tohourly(solar_radiation);%hourly%kw
clear nahavand
%@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@
%% solar power
%#######inputs##############inputs####################inputs##
%ambient temperature-hararate mohit
tamb2=temperature;%temperature of nahavand on each hour;
%p_npv=7.3;%kW >>0.130w each pv,14 pv in string;%rated power at reference condition
g=solar_radiation;%hourly_solar_radiation_nahavand';%kW
%############################################################
gref=1 ;%1000kW/m^2
tref=25;%temperature at reference condition
kt=-3.7e-3;% temperature coefficient of the maximum power(1/c0)
tc=tamb2+(0.0256).*g;
upv=0.986;%efficiency of pv with tilted angle>>98.6%
p_pvout_hourly=upv*(p_npv.*(g/gref)).*(1+kt.*(tc-tref));%output power(kw)_hourly
clear tc
clear tamb2
clear g
%% battery
%#######inputs##############inputs####################inputs##
%load demand >> hourly typical rural household load profile Kw
load1=[1.5 1 0.5 0.5 1 2 2 2.5 2.5 3 3 5 4.4 4 3.43 3 1.91 2.48 3 3 3.42 3.44 2.51 2];
load1=load1/5; load1=load1*2;% maximum would be 2kW mean is 1kW
%the load curve of a typical complete day's consumption, palestin case study
%load1=[1 1 1 1 1.5 2 2.5 2.5 1.5 1.5 1.5 2 2 2.5 2.5 3 3.5 4 4 3.5 2.5 1.5 1 1];
%houses=1;%number of houses in a village
load2=houses.*load1;%total load in a day for the whole village
%hourly load data for one year
a=0;
for i=1:1:360
a=[a,load2];
end
a(1)=[];
Pl1=a;
%ad=3;%daily autonomy
uinv=0.92;
ub=0.85;
dod=0.8;%depth of discharge 0.5 0.7 in the article 80%
el=mean(load2);
%bcap=40;%battery capacity 40 kWh
% ############################################################
cwh=(el*ad)/(uinv*ub*dod);%storage capacity for battery,bmax,kW
%% wind turbine
%#######inputs##############inputs####################inputs##
MCH=4;%choose the model>>model two
rw=cell2mat(WindTurbines(MCH,2));% rw=4;%blades diameter(m)
aw=cell2mat(WindTurbines(MCH,3));% aw=pi*(rw)^2;%Swept Area>>pi x Radius� = Area Swept by the Blades
uw=cell2mat(WindTurbines(MCH,4));% uw=0.95;%
vco=cell2mat(WindTurbines(MCH,5));% vco=25;%cut out
vci=cell2mat(WindTurbines(MCH,6));% vci=3;%cut in
vr=cell2mat(WindTurbines(MCH,7));% vr=8;%rated speed(m/s)
pr=cell2mat(WindTurbines(MCH,8));% pr=2;%rated power(kW)
pmax=cell2mat(WindTurbines(MCH,10));% % pmax=2.5;%maximum output power(kW)
pfurl=cell2mat(WindTurbines(MCH,9));% pfurl=2.5;%output power at cut-out speed9kW)
% ############################################################
v2=wind_speed;
%% diesel generator
Png=4;Png=nPng*Png;%kW output power of diesel generator
Bg=0.08145;%1/kW
Ag=0.246;%1/kW
Pg=4;Pg=nPng*Pg;%nominal power kW
% %fuel consumption of the diesel generator
Fg=Bg*Pg+Ag*Png;
%% MAIN PROGRAM
contribution=zeros(5,8640);%pv,wind, battery, diesel contribution in each hour
Ebmax=cwh;%40kWh%battery capacity 40 kWh
Ebmin=cwh*(1-dod);%40kWh
SOCb=0.2;%state of charge of the battery>>20%
Eb=zeros(1,8640);
time1=zeros(1,8640);
diesel=zeros(1,8640);
Edump=zeros(1,8640);
Edch=zeros(1,8640);
Ech=zeros(1,8640);
Eb(1,1)=SOCb*Ebmax;%state of charge for starting time
%^^^^^^^^^^^^^^START^^^^^^^^^^^^^^^^^^^^^^^^
Pl=Pl1;
clear Pl1;
%^^^^^^^^^^Out put power calculation^^^^^^^^
%solar power calculation
Pp=p_pvout_hourly;%output power(kw)_hourly
for i=1:1:8640
if Pp(i)>p_npv
Pp(i)=p_npv;%if the power output of pv exceed the maximum
end
end
% wind power calculation
for t=1:1:8640
%pr *((v2(t)-vci)/(vr-vci))^3pr+(((pfurl-pr)/(vco-vr))*(v2(t)-vr));
if v2(t)<vci %v2>>hourly_wind_speed;
pwtg(t)=0;
elseif vci<=v2(t)&& v2(t)<=vr
pwtg(t)=(pr/(vr^3-vci^3))*(v2(t))^3-(vci^3/(vr^3-vci^3))*(pr);
elseif vr<=v2(t) &&v2(t)<=vco
pwtg(t)=pr;
else
pwtg(t)=0;
end
Pw(t)=pwtg(t)*uw*nwt;%electric power from wind turbine
end

for t=2:1:8640
%^^^^^^^^^^^^^^READ INPUTS^^^^^^^^^^^^^^^^^^
%^^^^^^^^^^^^^^COMPARISON^^^^^^^^^^^^^^^^^^^
if Pw(t)+Pp(t)>=(Pl(t)/uinv)
%^^^^^^RUN LOAD WITH WIND TURBINE AND PV^^^^^^

if Pw(t)+Pp(t)>Pl(t)
%^^^^^^^^^^^^^^CHARGE^^^^^^^^^^^^^^^^^^^^^^^^^^

[Edump,Eb,Ech] = charge(Pw,Pp,Eb,Ebmax,uinv,Pl,t,Edump,Ech);
time1(t)=1;

contribution(1,t)=Pp(t);contribution(2,t)=Pw(t);contribution(3,t)=Edch(t);contribution(4,t)=diesel(t);c
ontribution(5,t)=Edump(t);
else
Eb(t)=Eb(t-1);
return
end

else
%^^^^^^^^^^^^^^DISCHARGE^^^^^^^^^^^^^^^^^^^
[Eb,Edump,Edch,diesel,time1,t] =
dicharge(Pw,Pp,Eb,Ebmax,uinv,Pl,t,Pg,Ebmin,Edump,Edch,Ech,diesel,time1);

contribution(1,t)=Pp(t);contribution(2,t)=Pw(t);contribution(3,t)=Edch(t);contribution(4,t)=diesel(t);c
ontribution(5,t)=Pl(t);
end
end
%% plotting
% figure
a=contribution';
b=sum(a);
renewable_factor=b(4)/(b(1)+b(2));
% h=pie(b);
% colormap jet;
% legend('PV','WIND','BATTERY','DIESEL');
%reliability
%lose of load probability=sum(load-pv-wind+battery)/sum(load)
total_loss=0;
for t=1:1:8640
% aa(t)=Pl(t)-Pp(t)-Pw(t)+Eb(t);
if Pl(t)>(Pp(t)+Pw(t)+(Eb(t)-Ebmin)+diesel(t))
total_loss=total_loss+(Pl(t)-(Pp(t)+Pw(t)+(Eb(t)-Ebmin)+diesel(t)));
end

end
LPSP=total_loss/(sum(Pl));
% reliability=sum(aa)/sum(Pl);
% [price_electricity] = economic(diesel,Pl,Fg,cwh);
[price_electricity] = economic_fast_iran(diesel,Pl,Fg,cwh,p_npv,nwt,houses);
ali=[Pp(1:168)',Pw(1:168)',Eb(1:168)',diesel(1:168)',Pl(1:168)',Edump(1:168)'];
Edump=sum(Edump);
end