You are advised to do some reading on the PSO algorithm to understand the code.
% Step 1: Initialize PSO
Xmin = -4; % lower bound of candidate solutions
Xmax = 4; % upper bound of candidate solutions
Vmax = 1; % maximum moving distance of a particle in a loop
Vmin = -1; % minimum moving distance of a particle in a loop
c1 = 1.3; % Local learning factor
c2 = 1.7; % Global learning factor
wmin = 0.10; % min weight value
wmax = 0.90; % max weight value
Gmax = 100; % number of generations (fancy term for computational iterations)
Size = 50; % number of particles (fancy term for candidate solutions)
for i = 1:Gmax
w(i) = wmax - ((wmax - wmin)/Gmax)*i;
end
for i = 1:Size
for j = 1:2
x(i, j) = Xmin + (Xmax - Xmin)*rand(1);
v(i, j) = Vmin + (Vmax - Vmin)*rand(1);
end
end
% Step 2: Calculte fitness
for i = 1:Size
p(i) = fitfun(x(i, :));
y(i, :) = x(i, :);
if i == 1
plocal(i, :) = findlbest(x(Size, :), x(i, :), x(i+1, :));
elseif i == Size
plocal(i, :) = findlbest(x(i-1, :), x(i, :), x(1, :));
else
plocal(i, :) = findlbest(x(i-1, :), x(i, :), x(i+1, :));
end
end
BestS = x(1,:);
for i = 2:Size
if fitfun(x(i, :)) > fitfun(BestS)
BestS = x(i, :);
end
end
% Step 3: Main loop
for kg = 1:Gmax
for i = 1:Size
% local adaptive mutation operator
M = 1;
if M == 1
v(i, :) = w(kg)*v(i, :) + c1*rand*(y(i, :) - x(i, :)) + c2*rand*(plocal(i, :) - x(i, :)); % Local optimization
elseif M == 2
v(i, :) = w(kg)*v(i, :) + c1*rand*(y(i, :) - x(i, :)) + c2*rand*(BestS - x(i, :)); % Global optimization
end
for j = 1:2 % Evaluate the velocity
if v(i, j) < Vmin
v(i, j) = Vmin;
elseif x(i, j) > Vmax
v(i, j) = Vmax;
end
end
x(i, :) = x(i, :) + v(i, :)*1; % Update position
for j = 1:2 % Check the limit
if x(i, j) < Xmin
x(i, j) = Xmin;
elseif x(i,j) > Xmax
x(i,j) = Xmax;
end
end
% Adaptive mutation
if rand > 0.60
k = ceil(2*rand);
x(i, k) = Xmin + (Xmax - Xmin)*rand(1);
end
% Step 4: Evaluate and update
if i == 1
plocal(i, :) = findlbest(x(Size, :), x(i, :), x(i+1, :));
elseif i == Size
plocal(i, :) = findlbest(x(i-1, :), x(i, :), x(1, :));
else
plocal(i, :) = findlbest(x(i-1, :), x(i, :), x(i+1, :));
end
if fitfun(x(i, :)) > p(i) % Evaluate and update
p(i) = fitfun(x(i, :));
y(i, :) = x(i, :);
end
if p(i) > fitfun(BestS)
BestS = y(i, :);
end
end
Best_value(kg) = fitfun(BestS);
end
figure(1);
kg = 1:Gmax;
plot(kg, Best_value, 'r', 'linewidth', 2, 'Color', "#F7A7B7"), grid on
xlabel('Generations');
ylabel('Fitness function');
display('Local Maximum x at'); disp(BestS);
display('with the value of the fitness function = '); disp(Best_value(Gmax));
[X, Y] = meshgrid(-4:8/40:4);
Z = (3*cos(X.*Y) + X + Y.^2);
surf(X, Y, Z), hold on
plot3(BestS(1), BestS(2), Best_value(Gmax), 'pentagram', 'MarkerFaceColor','red', 'MarkerSize', 15),
xlabel('x'), ylabel('y'), zlabel('f(x,y)')
% Fitness Function
function f = fitfun(x)
f = 3*cos(x(1)*x(2)) + x(1) + x(2)^2;
end
% Find Local Best
function f = findlbest(x1, x2, x3)
xSet = [x1; x2; x3];
fSet = [fitfun(x1), fitfun(x2), fitfun(x3)];
[maxvalue index] = max(fSet);
plocalbest = xSet(index,:);
f = plocalbest;
end
