PSO optimization for rastrigin function

Could anyone give me a hint on how to optimize the reuslts of the code below, of PSO applied to rastrigin function? For now I'm getting an average of 14.5 as result of 10 runs:

pso.m

clear vars;
close all;
clc;

graph = 1;

maxFX = 1e3;

% ***** CHANGE IN SEARCH AREA, LARGER
xmin = -10.24;
xmax = 10.24;

numVAR = 5; % 5 dimensions.
numPOP = 20;

% start of TIC TOC function to calculate processing time:
tic

POP = xmin + rand (numPOP, numVAR) * (xmax - xmin);
FX = trace (POP);
numFX = numPOP;

% random speed, where to go and how fast:
V = 0.5 * rand (numPOP, numVAR) * (xmax - xmin); % ***** CAN CHANGE THIS VARIABLE
% at most 50% of the size of the search space in each dimension.

% ***** CAN CHANGE THIS VARIABLE
Vmax = 0.2 * (xmax - xmin);

% Definition of PBEST
POPp = POP;
FXp = FX;

% GBEST Definition - JUST ONE, THE LOWEST EVER.
[~, ind] = min (FXp); % Position of lowest F (x)
POPg = POPp (ind, :);
FXg = FXp (ind);

% attraction values:
alpha = [1 0.2 0.2];

% the population will walk along the search space using a
% vector of speed that will be updated each season.


% V = Alpha + V_Current + (Beta1 * DistPbest) + (Beta2 * DistGbest)
while (numFX <maxFX)
    alpha (1) = 1 - 0.8 * numFX / maxFX;
    
    % updating population displacement:
    POP = POP + V;
    POP = max (POP, xmin);
    POP = min (POP, xmax);
    
    FX = trace (POP);
    numFX = numFX + numPOP;
    
    % for each individual, you need to verify that the generated FX value is
    % better than Pbest
    for i = 1: numPOP
        if (FX (i) <= FXp (i))
            POPp (i, :) = POP (i, :);
            FXp (i) = FX (i);
            
            % checks whether to update Gbest as well:
            if (FX (i) <= FXg)
                POPg = POP (i, :);
                FXg = FX (i);
            end
        end
        
        % updating vector V:
        V (i, :) = alpha (1) * V (i, :) ...% Inertia coefficient
            + alpha (2) * (POPp (i, :) - POP (i, :)) ...% Distance to PBEST
            + alpha (3) * (POPg - POP (i, :)); % Distance to GBEST
        % alpha (2) = Beta1, alpha (3) = Beta2
        % [ith] Beta1 = distance from Pbest to me: me (POPp) - where I am.
        % [ith] Beta2 = distance from Gbest to me: Gbest (POPg) - where am I
    
    end
    % max between a vector (V) and a scalar (Vmax) = each cell of V is
    % compared with Vmax
    V = min (Vmax, V);
    V = max (-Vmax, V);
    Final% of TIC TOC function to calculate processing time in seconds:
    tempoPROCESS = toc
    if (graphic == 1)
        clf;
        hold on;
        plot (POP (:, 1), POP (:, 2), 'bo'); % blue balls (particles, POP)
        plot (POPp (:, 1), POPp (:, 2), 'rx'); % red crosses (PBest)
        plot (POPg (1), POPg (2), 'gv'); % green triangle (GBest)
        axis ([xmin xmax xmin xmax]);
        xlabel (num2str (FXg));
        drawnow;
    end
end

rastrigin.m

function FX = trace (POP)
     [nPOP, dim] = size (POP);

     FX = ones (nPOP, 1) * 10 * dim;

     for d = 1: dim
         FX = FX + POP (:, d). ^ 2-10 * cos (2 * pi * POP (:, d));
     end
     FX = FX;
end

Image result of one of the runs: