function [gbest,gbestval,fitcount]= PSO_func(fhd,Max_Gen,Max_FES,Particle_Number,Dimension,VRmin,VRmax,inertia_gen,varargin)
%[gbest,gbestval,fitcount]= PSO_func('hemant5',500,2000,30,2,0,1)
%rand('twister',sum(100*clock));
rand('twister',8399);
me=Max_Gen;
ps=Particle_Number;
D=Dimension;
cc=[1.4945 1.4945];%acceleration constants
%cc=[2.5 1.5];
mvcrit=1;
cnt3=0;
cnt2=0;
ergrd=1e-9;
ergrdep=40;
%cc=[2.5 1.5];
iw1=0.9;iw2=0.4;iwe=inertia_gen;
%iwt=0.9-(1:me).*(0.5./me);
if length(VRmin)==1
VRmin=repmat(VRmin,1,D);
VRmax=repmat(VRmax,1,D);
end
%mv=(VRmax-VRmin);
mvmax=VRmax;
mvmin=-mvmax;
VRmin=repmat(VRmin,ps,1);
VRmax=repmat(VRmax,ps,1);
%Vmin=repmat(-mv,ps,1);
Vmin=repmat(mvmin,ps,1);
Vmax=repmat(mvmax,ps,1);
%Vmax=-Vmin;
%Vmin=VRmin;
%Vmax=VRmax;
pos=VRmin+(VRmax-VRmin).*rand(ps,D);
ccmax=2.0;ccmin=1.3;
for i=1:ps;
e(i,1)=feval(fhd,pos(i,:),varargin{:});
end
fitcount=ps;
vel=Vmin+(Vmax-Vmin).*rand(ps,D);%initialize the velocity of the particles
pbest=pos;
pbestval=e; %initialize the pbest and the pbest's fitness value
[gbestval,gbestid]=min(pbestval);
gbest=pbest(gbestid,:);%initialize the gbest and the gbest's fitness value
gbestrep=repmat(gbest,ps,1);
tr(1)=gbestval;
average(1)=mean(e);
for i=1:me %start iteration loop
%iwt(i)=0.5+rand(1,1)/2;
if i<=iwe
iwt(i)=((iw2-iw1)/(iwe-1))*(i-1)+iw1; % get inertia weight, just a linear funct w.r.t. epoch parameter iwe
else
iwt(i)=iw2;
end
Vmax1=((0.1*Vmax-Vmax)/(me-1))*(i-1)+Vmax;
Vmin1=-Vmax1;
% cc1(i)=(ccmin-ccmax)/(me-1)*(i-1)+ccmax;
%cc2(i)=(ccmax-ccmin)/(me-1)*(i-1)+ccmin;
r1=rand(1,D);
if r1<=0.05,T=-1;else T=1;end
%if r1>0.05,sign(r1)=1;end
for k=1:ps %start particle loop
aa(k,:)=cc(1).*rand(1,D).*(pbest(k,:)-pos(k,:))+cc(2).*rand(1,D).*(gbestrep(k,:)-pos(k,:));
vel(k,:)=T.*iwt(i).*vel(k,:)+aa(k,:);
for dimcnt=1:D
if vel(k,dimcnt)>=Vmax1(k,dimcnt)
vel(k,dimcnt)=Vmax1(k,dimcnt);
end
if vel(k,dimcnt)<=Vmin1(k,dimcnt)
vel(k,dimcnt)=Vmin1(k,dimcnt);
end
end
dimcnt=1:D;
vel(k,:)=min(Vmax1(k,dimcnt),max(Vmin1(k,dimcnt),vel(k,:)));
%vel(k,:)=(vel(k,:)>mv).*mv+(vel(k,:)<=mv).*vel(k,:);
%vel(k,:)=(vel(k,:)>mvmax).*mvmax+(vel(k,:)<=mvmax).*vel(k,:);
%vel(k,:)=(vel(k,:)<(-mv)).*(-mv)+(vel(k,:)>=(-mv)).*vel(k,:);
%vel(k,:)=(vel(k,:)<(mvmin)).*(mvmin)+(vel(k,:)>=(mvmin)).*vel(k,:);
pos(k,:)=pos(k,:)+vel(k,:);
pos(k,:)=((pos(k,:)<=VRmax(k,:))&(pos(k,:)>=VRmin(k,:))).*pos(k,:)...
+(pos(k,:)>VRmax(k,:)).*VRmax(k,:)+(pos(k,:)<VRmin(k,:)).*VRmin(k,:);
%pos(k,:)=((pos(k,:)<=VRmax(1,:))&(pos(k,:)>=VRmin(1,:))).*pos(k,:)+(1-(pos(k,:)<=VRmax(1,:))&(pos(k,:)>=VRmin(1,:))).*(VRmin(1,:)+(VRmax(1,:)-VRmin(1,:)).*rand(1,D));
% pos(k,:)=((pos(k,:)>=VRmin(1,:))&(pos(k,:)<=VRmax(1,:))).*pos(k,:)...
% +(pos(k,:)<VRmin(1,:)).*(VRmin(1,:)+0.25.*(VRmax(1,:)-VRmin(1,:)).*rand(1,D))+(pos(k,:)>VRmax(1,:)).*(VRmax(1,:)-0.25.*(VRmax(1,:)-VRmin(1,:)).*rand(1,D));
%for dimcnt=1:D
%if pos(k,dimcnt)>VRmax(k,dimcnt)
% pos(k,dimcnt)=VRmax(k,dimcnt);
% vel(k,dimcnt)=-rand(1)*vel(k,dimcnt);
%vel(k,dimcnt)=Vmax(k,dimcnt);
%end
% if pos(k,dimcnt)<VRmin(k,dimcnt)
% pos(k,dimcnt)=VRmin(k,dimcnt);
% vel(k,dimcnt)=-rand(1)*vel(k,dimcnt);
%vel(k,dimcnt)=Vmin(k,dimcnt);
%end
%end
% if (sum(pos(k,:)>VRmax(k,:))+sum(pos(k,:)<VRmin(k,:)))==0;
e(k,1)=feval(fhd,pos(k,:),varargin{:});
fitcount=fitcount+1;
tmp=(pbestval(k)<e(k));
temp=repmat(tmp,1,D);
pbest(k,:)=temp.*pbest(k,:)+(1-temp).*pos(k,:);
pbestval(k)=tmp.*pbestval(k)+(1-tmp).*e(k);%update the pbest
if pbestval(k)<gbestval
gbest=pbest(k,:);
gbestval=pbestval(k);
gbestrep=repmat(gbest,ps,1);%update the gbest
end
% end
tr(i+1)=gbestval;
end %end particle loop
iteration=i
gbestval
average(i+1)=mean(e);
%% this section is for modifying the maximum allowable velocity
%
%% dynamically decrease maximum velocity if gbestval
%% doesn't change over mvcrit iterations (set at top of function)
% Vmax=0.97*Vmax;
%Vmin=-Vmax;
%if tr(i)~=gbestval
%cnt3=0;
%elseif tr(i)==gbestval
% cnt3=cnt3+1;
% if cnt3>=mvcrit
% Vmax=Vmax*.999;
% Vmin=Vmin*.999;
% cnt3=0;
%end
%end
% check for stopping criterion based on speed of convergence to desired error
%if tr(i)~=gbestval
%cnt2=0;
%elseif tr(i)==gbestval
% cnt2=cnt2+1;
% if cnt2>=ergrdep
% disp('***************************** global error gradient too small for too long');
% break
%end
%end
%tmp1=abs(tr(i)-gbestval);
%if tmp1>ergrd
% cnt2=0;
%elseif tmp1<=ergrd
% cnt2=cnt2+1;
%if cnt2>=ergrdep
% disp('***************************** global error gradient too small for too long');
% break
%end
%end
% if round(i/20)==i/20
% plot(pos(:,1),pos(:,2),'b*');hold on
% plot(pbest(:,1),pbest(:,2),'r*');hold off
% title(['PSO: ',num2str(i),' generations, Gbestval=',num2str(gbestval)]);
% axis([-100,100,-100,100]);
% drawnow
% end
%if fitcount>=Max_FES
% break;
%end
if gbestval==0
break
end
end %end iteration loop
figure(1);
plot([0,1:iteration],tr);
figure(2);
plot([1:iteration],iwt);
