clc;
clear;
disp('Particle Swarm Optimization Algorithm');
pause(1)
disp(' ');
npop=20;
nvar=3;
w=1;
wdamp=0.99;
c1=2;
c2=2;
pmin=[ 60 70 1];
pmax=[ 150 100 26 ];
LB=pmin;
UB=pmax;
dx=pmax-pmin;
vmax=0.1*dx;
maxit=100;
empty_particle.position=[];
empty_particle.velocity=[];
empty_particle.cost=[];
empty_particle.pbest=[];
empty_particle.pbestcost=[];

particle=repmat(empty_particle,npop,1);

gbest=zeros(maxit,nvar);
gbestcost=zeros(maxit,1);

for it=1:maxit
    if it==1
        gbestcost(1)=inf;
        for i=1:npop
            particle(i).velocity=zeros(1,nvar);
            particle(i).position=LB+rand.*(UB-LB);
            
            [F ]=tabe2(particle(i).position);
            particle(i).cost=F;
            
            particle(i).pbest=particle(i).position;
            particle(i).pbestcost=particle(i).cost;
            
            if particle(i).pbestcost<gbestcost(it)
                gbest(it,:)=particle(i).pbest;
                gbestcost(it)=particle(i).pbestcost;
            end
        end
    else
        gbest(it,:)=gbest(it-1,:);
        gbestcost(it)=gbestcost(it-1);
        for i=1:npop
            particle(i).velocity=w*particle(i).velocity...
                                +c1*rand*(particle(i).pbest-particle(i).position)...
                                +c2*rand*(gbest(it,:)-particle(i).position);
                            
            particle(i).velocity=min(max(particle(i).velocity,-vmax),vmax);
            
            particle(i).position=particle(i).position+particle(i).velocity;
            
               for j=1:nvar
            if j<169
            if particle(i).position(j)<LB(j)
                particle(i).position(j)=LB(j);
            end
            if particle(i).position(j)>UB(j)
                particle(i).position(j)=UB(j);
            end 
            else
              if particle(i).position(j)<LB(j)
                particle(i).position=LB(j);
            end
            if particle(i).position(j)>UB(j)
                particle(i).position(j)=UB(j);
            end   
              particle(i).position(j)=round(particle(i).position(j));  
            end
        end
            
            [F ]=tabe2(particle(i).position);
           
            FZ(it,:)=[F];
            particle(i).cost=F;
            
            if particle(i).cost<particle(i).pbestcost
                particle(i).pbest=particle(i).position;
                particle(i).pbestcost=particle(i).cost;

                if particle(i).pbestcost<gbestcost(it)
                    gbest(it,:)=particle(i).pbest;
                    gbestcost(it)=particle(i).pbestcost;
                end
            end
        end

    end
    
    disp(['Iteration ' num2str(it) ':   Best Cost = ' num2str(gbestcost(it))]);
    
    w=w*wdamp;
end

plot(gbestcost,'*');