Algoritmo PSO

Se coloca  el código en Matlab de una simulación para el algoritmo PSO simple. Este algoritmo hace la evaluación de una función objetivo. 

Desde acá se evalua el el algoritmo

clear all
clc
N=250;% Numero de particulas
maxit=100; % Numero de pasadas
dim =1;
upbnd = 90; % Limite superior de inicio del swarm
lwbnd = 120; % Limite inferior de inicio del swarm
% Inicializa velocidad y posición del swarm.
x = rand(N,dim)*(upbnd-lwbnd) + lwbnd ;
v = rand(N,dim); %velocidad inicial
plot(x,v,'.r')
hold on
title('Posición Vs Velocidad')
grid on
xlabel('Posición')
ylabel('Velocidad')

[brs,kol]=size(x);
f=zeros(N,1);
rhomax=0.9;rhomin=0.4;
for it=1:maxit
    rho(it)=rhomax-((rhomax-rhomin)/maxit)*it;
end
%Evaluación de la función objetivo en cada particula
for i=1:brs
    f(i)=fungsi2(x(i,:));
end
%Inicia el proceso de comparación
it=1;
Pbest=x;
fbest=f;
[minf,idk]=min(f);
Gbest=x(idk,:);
lastbest=[0 0];
minftot=[];
while it<maxit
    r1=rand;r2=rand;
    for j=1:brs
        v(j,:)=rho(it).*v(j,:)+r1.*(Pbest(j,:)-x(j,:))+r2.*(Gbest-x(j,:));
        x(j,:)=x(j,:)+v(j,:);
        f(j)=fungsi2(x(j,:)); % Se evalua la función objetivo
     end
    %Se actualiza Pbest
    changerow = f < fbest;
    fbest=fbest.*(1-changerow)+f.*changerow;
    Pbest(find(changerow),:)=x(find(changerow),:);
    [minf,idk]=min(fbest);
    minftot=[minftot;minf];
    Gbest=Pbest(idk,:);
    if sum(var(Pbest))<1e-8
        break
    end
    it=it+1;
    lastbest=Gbest;
    figure(1)
    plot(x,v,'.');
    pause(0.01)
    legend('Poblacion Inicial','Poblacion final')
   
end
 print -deps -r100 resul
x(1,3)=0;
v(1,3)=0;
pos=x;
vel=v;
grafica(N,it,v(length(v)),pos,vel)
xopt=Gbest;
fmin=minf;
figure(3)
plot(minftot,'r*')
grid on
title('Valores optimos')
print -deps -r100 Valores1

Las funciones que acompañan este código son las siguientes:

Función objetivo

function f=fungsi2(x)
f=2.1048*(4*(100000)/(pi*x^(2)));

Función que grafica en 3D las particulas

function grafica(N,tmax,vMax,r,v)
a=zeros(N,3);  %accelerations
L=10;     % Tamaño del cubo
t2=0;     %Variable que toma el tiempo en cada paso
[r,v]=inicialize(L,N,r,v,vMax);
dt=0.01;
for i=1:tmax
    t2=t2+1;  %Así vario el tiempo
              %Acá viene el equivalente a la función velocityVerlet(double dt)
   [a]=computeAccelerations(L,N,r,v,a,vMax);
   for i=1:N
       for k=1:3
           r(i,k)= r(i,k)+v(i,k)*dt +0.5*a(i,k)*dt*dt;
           v(i,k)= v(i,k) + 0.5*a(i,k)*dt;
       end
   end   
   [a]=computeAccelerations(L,N,r,v,a,vMax);
   for i=1:N
       for k=1:3
           v(i,k)=v(i,k) + 0.5*a(i,k)*dt;
       end
   end
 

   figure(2)
   plot3(a,r,v,'.');
   grid on;
   title('Movimiento Particulas última Población')
    xlabel('Aceleración');
    ylabel('Posicion');
    zlabel('Velocidad')
    drawnow;
end
print -deps -r100 tempe3D

Función que calcula la aceleración de las particulas

function [a]=computeAccelerations(L,N,r,v,a,vMax)
a=zeros(N,3); %Primero que todo siempre pone en ceros la matriz acceleracion
for i=1:N-1 %loop over all distinct pairs i,j
    for j=i+1:N %j=i+1 es recorrer la Diagonal superior.
        rij=zeros(1,3); %un arreglo que guarda la position of i relative to j 
        rSqd=0; %posicion al cuadrado
        for k=1:3
            rij(1,k)=r(i,k)-r(j,k);
            rSqd=rSqd+rij(1,k)*rij(1,k);%Lennard Jones
        end
        %Aqui le meto un potencial de interaccion luego es un gas real
        f=24*(2*rSqd^(-7)-rSqd^(-4));
        for k=1:3
            a(i,k)=a(i,k)+rij(1,k)*f;
            a(j,k)=a(j,k)-rij(1,k)*f;
        end
    end
end

Se adjunta archivo con la explicación del algortimo y los resulatdos.