Particle Swarm Optimization(PSO)

Particle Swarm Optimization란?

Particle Swarm method는 동물들의 군집화 행동에서 영감을 받아 1995년에 Social Psychologist 'James Kennedy'와 Electrical engineer인 Russel C.Eberhart에 의해 만들어진 알고리즘입니다.Swarm은 많은 입자들의 움직임에 의해 만들어지고, 각각 입자의 움직임은 다음과 같이 1)2)3)을 고려해 업데이트 됩니다. 

1)Current speed : Adventurous tendency

2)Its personal experience : conservative tendency to come back as its best position

3)The social experience : Panurgian tendency to follow a group of neighbors

조금 더 자세하게 살펴보도록 하겠습니다. Paricle Swarm 최적화 방법에서는 각각의  입자(particle)들이 3개의 vector를 가지고 있습니다. 속도(Velocity)는 각각의 Personal Best와 Global Best를 고려해 업데이트 됩니다.

4가지 중요한 control parameter :

• Swarm size  

• Inertial Weight W 

• Acceleration coefficeint c1,c2

• Number of interations

입자가 크고 무거울수록 더 좋은 Global search ability를 가집니다. 반면, 작고 가벼운 입자들은 Local search ability를 가집니다. Acceleration coeefficent는 c1>c2일 경우, glabal search ability가 더 뛰어납니다. 반대로 c1<c2일경우, Local search ability가 더 뛰어납니다. Number of iteration도 정확한 해를 구하기 위한 중요한 조건이 됩니다. 더 많은 loop을 돌릴 수록 Solution에 가까워 지겠죠.

 

Particle Swarm의 Algorithm 순서도

Particle Swarm Algorithm

 

 

Ragistrigin function f에 대해 minimum 을 구하도록 하겠습니다.

$ f\left(x,y\right)=20+x^2 +y^2 -10\left(\mathrm{cos}\left(2\pi x\right)+\mathrm{cos}\left(2\pi y\right)\right)$

아래와 같이 parameter를 설정했습니다.

1.  Number of particles = 40

2. Maximum interation = 90

3. C1,C2=2

$W=W_{\mathrm{max}} -\left(\frac{W_{\mathrm{max}} -W_{\mathrm{min}} }{K_{\mathrm{max}} }\right)\cdot K$

$w_{\mathrm{min}} =0\ldotp 4$

$w_{\mathrm{max}} =0\ldotp 9$

$K_{\mathrm{max}} \;\mathrm{is}\;\mathrm{the}\;\mathrm{maximum}\;\mathrm{iteration}\;\mathrm{number},K\;\mathrm{is}\;\mathrm{the}\;\mathrm{current}\;\mathrm{iterations}\;\mathrm{number}$

Simulation 결과

Simulation result

------------------------------------------------------------------------------------------------
our Local Best Position is
LBP = 2×40    
10^(-3) x

  -0.0120    0.0089    0.0151    0.0124   -0.0555   -0.0034    
   0.0076   -0.0024    0.0160   -0.0917   -0.2094    0.0061    
   0.0103   -0.0023   -0.0079   -0.0573    0.0134    0.1938    
   0.0835    0.0093   -0.0560   -0.4339    0.0064    0.0208    
   0.0035   -0.0118    0.0042   -0.0981    0.1485    0.0060   
  -0.0253    0.0064    0.0109    0.0030    0.0468    0.1617    
   0.0084    0.0090   -0.0811   -0.0033  0.0082    0.0010    
   0.0191   -0.0152    0.0089   -0.0023   -0.0035    0.0293   
  -0.0239    0.0293   -0.0927    0.0047    0.0429    0.2110    
   0.0140   -0.0612   -0.0000    0.0472    0.1607    0.0061    
   0.0149   -0.7499   -0.0095    0.0047    0.0137   -0.0076   
  -0.0594    0.0735    0.1812    0.0106   -0.0275    0.0041   
  -0.0043    0.0323   -0.0868   -0.1122    0.0340   -0.0078
  -0.2272    0.0209

with fitness values =
LBF = 1×40    
10^(-3) x

    0.0000    0.0000    0.0001    0.0001    0.0006    0.0000    
    0.0000    0.0002    0.0002    0.0018    0.0104    0.0000    
    0.0004    0.0088    0.0001    0.0014    0.0000    0.0079    
    0.0065    0.0000    0.0007    0.1489    0.0000    0.0001    
    0.0000    0.0000    0.0007    0.0030    0.0109    0.0000    
    0.0003    0.0000    0.0000    0.0002    0.0019    0.0077    
    0.0002    0.0000    0.0115    0.0001

------------------------------------------------------------------------------------------------

Matlab으로 구현한 코드입니다.

clear all;
close all;
clc;
 
syms x y % we made symbolic variable x1 and x2
fun = 20 + x^2 + y^2 - 10*(cos(2*pi*x) + cos(2*pi*y))
 
figure(1)
[x1, y1] = meshgrid(-2:0.05:2); 
surf(x1, y1, 20 + x1^2 + y1^2 - 10*(cos(2*pi*x1) + cos(2*pi*y1))); %we plot function to see what it looks like
 
figure(2)
[x2, y2] = meshgrid(-0.1:0.01:0.1); 
surf(x2, y2, 20 + x2^2 + y2^2 - 10*(cos(2*pi*x2) + cos(2*pi*y2))); %we plot function to see what it looks like
 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%% PSO %%%%%%%%%%%%%%%%%%%%%%%%%%%%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
 
%%%%%%%%%%%%%%%%%%%%% PSO parameters %%%%%%%%%%%%%%%%%%%%%%%%
Nbr_part = 40;
Max_it = 80;
C1 = 2;
C2 = 2;
 
Wmin = 0.4;
Wmax = 0.9;
%%%%%%%%%%%%%%%%%%%%% RESTRICTION %%%%%%%%%%%%%%%%%%%%%%%%%%
%%  fun = 20 + x^2 + y^2 - 10*(cos(2*pi*x) + cos(2*pi*y)) %%
Xres = [-0.1 0.1];
Yres = [-0.1 0.1];
 
%%%%%%%% Generating random for initial values %%%%%%%%%%%%%
 
R1 = rand(1,Nbr_part);
R1 = R1-0.5; %TO DISTRIBUTE POSITIVE AND NEGATIVE VALUES 
 
R2 = rand(1,Nbr_part);
R2 = R2-0.5;
 
%%%%%%%%%%%%%%%%%% DEFINE FITNESS %%%%%%%%%%%%%%
%%%%%%%%%%%%%%%% PARTICLESS INITIAL POS %%%%%%%%%%%%%% 
%%%%%%%%%%%%%%%%% INITIAL VELOCITY%%%%%%%%
 
FITNESS = zeros(1,Nbr_part); %fitness set to zero at beginning
 
CP(1,:) = 1/10 * R1; %we multiply by 1/10 to be in the restriction range
CP(2,:) = 1/10 * R2;
 
V(1,:) = R1;
V(2,:) = R2;
 
%%%%%%%%%%%%%%%%%%%%% ITERATION 1 %%%%%%%%%%%%%%%%%%%%%%%
 
FITNESS = double(subs(fun,{x y},{CP(1,:) CP(2,:)}));
 
%%%%%%%%%%% LBP, LBF GBP AND GBF CALCULATION %%%%%%%%%%%%%%%%
LBP = CP;
LBF = FITNESS;
 
[GBF,index] = min(FITNESS);
GBP(1,:) = repelem(CP(1,index),Nbr_part);
GBP(2,:) = repelem(CP(2,index),Nbr_part);
 
%%%%%%%%%%%%%%%%% ITERATION 2 %%%%%%%%%%%%%%%%%%%%%%%%%
iteration = 2;
%disp('-------------------------------------------------')
 
while(iteration<Max_it)
 
R1 = rand(1,Nbr_part);
R2 = rand(1,Nbr_part);
    
W = calcW(iteration,Max_it,Wmin,Wmax);
V= W*V + C1*R1.*(LBP-CP)+ C2*R2.*(GBP-CP);
 
%%%%%%%% confinement algorithm using Thales theorem %%%%%
xy_temp = CP + V;
 
    for i=1:Nbr_part    
        if (xy_temp(1,i)>max(Xres) || xy_temp(1,i)<min(Xres) || xy_temp(2,i)>max(Yres) || xy_temp(2,i)<min(Yres))
            abs_x = abs(xy_temp(1,i));
            abs_y = abs(xy_temp(2,i)) ;
            if(abs_x>abs_y)
                if(xy_temp(1,i)>0)
                    xy_temp(1,i)=max(Xres);
                else 
                    xy_temp(1,i)=min(Xres);
                end
                xy_temp(2,i)= (((max(Yres)-min(Yres))/2)*xy_temp(2,i))/abs_x;
            else
                if(xy_temp(2,i)>0)
                        xy_temp(2,i)=max(Yres);
                else 
                    xy_temp(2,i)=min(Yres);
                end
                xy_temp(1,i)= (((max(Xres)-min(Xres))/2)*xy_temp(1,i))/abs_y;
            end
        end   
    end
    
    CP = xy_temp; %new CP with all particules within the area
    
    FITNESS = double(subs(fun,{x y},{CP(1,:) CP(2,:)}));
    
    LBF_prev_it = LBF;
    
    LBF = min(FITNESS,LBF); %LBF updated to still keep the best result
    [GBF,index]= min(LBF);
    
    GBP(1,:) = repelem(CP(1,index),Nbr_part);
    GBP(2,:) = repelem(CP(2,index),Nbr_part);
    
    for i=1:Nbr_part
        if (LBF(i)<LBF_prev_it(i));
            LBP(1,i)=CP(1,i);
            LBP(2,i)=CP(2,i);
        end
    end
    iteration = iteration +1;
%     LBP
    %disp('-------------------------------------------------')
 
end
 
hold on;
plot3(GBP(1,1), GBP(2,1), GBF,'.y','markersize',40);
title('function representation and last Global Best Position')
 
disp('---------------------------------')
disp('our Local Best Position is')
LBP
disp('with fitness values =')
LBF
 
disp('---------------------------------')
disp('our Global Best Position is')
GBP(:,1)
disp('with a fitness value =')
GBF
 
 

 

출처 : Paris Saclay university Computational method class