clc;clear;close all;
%% initialization
%
swarm_size = 64; % number of the swarm particles
maxIter = 70; % maximum number of iterations
inertia = 1;
correction_factor = 1.5;
% set the position of the initial swarm
a = 1:8;
[X,Y] = meshgrid(a,a);
C = cat(2,X',Y');
D = reshape(C,[],2);
Position= D; % set the position of the particles in 2D
Velocity= zeros(size(D)); % set initial velocity for particles
Pbest_Pos = zeros(size(D));
Pbest_Value=ones(size(D,1),1)*inf; % set the best value so far
plotObjFcn = 1; % set to zero if you do not need a final plot
%% define the objective funcion here (vectorized form)
objfcn = @(x)(x(:,1) - 20).^2 + (x(:,2) - 25).^2;
change_Pos=@(old,velo)(old+velo);
change_Velocity=@(a,b,c,velo,oldpos,pbestpos,gbestpos)...
(a*rand(size(velo)).*velo+...
b*rand(size(pbestpos)).*(pbestpos-oldpos)+...
c*rand(size(gbestpos)).*(gbestpos-oldpos));
%% The main loop of PSO
for iter = 1:maxIter
Position= change_Pos(Position,Velocity); %update position with the velocity
fval = objfcn(Position); % evaluate the function using the position of the particle
% compare the function values to find the best ones
Pbest_Pos(fval<Pbest_Value,:) =Position(fval<Pbest_Value,:);
Pbest_Value(fval<Pbest_Value,:) =fval(fval<Pbest_Value,:);
[~, gbestMark] = min(Pbest_Value); % find the best function value in total
Gbest_Pos=Pbest_Pos(gbestMark*ones(swarm_size,1),:);
% update the velocity of the particles
Velocity=change_Velocity(inertia,correction_factor,correction_factor,...
Velocity,Position,Pbest_Pos,Gbest_Pos);
% plot the particles
clf;plot(Position(:,1), Position(:,2), 'bx'); % drawing swarm movements
axis([-2 40 -2 40]);
pause(.1); % un-comment this line to decrease the animation speed
disp(['iteration: ' num2str(iter)]);
end
%% plot the function
if plotObjFcn
ub = 40;
lb = 0;
npoints = 1000;
x = (ub-lb) .* rand(npoints,2) + lb;
for ii = 1:npoints
f = objfcn([x(ii,1) x(ii,2)]);
plot3(x(ii,1),x(ii,2),f,'.r');hold on
end
plot3(Pbest_Pos(:,1),Pbest_Pos(:,2),Pbest_Value,'xb','linewidth',5,'Markersize',5);grid
end