Write MATLAB code to solve the following the minimum problem.

Given the function f(x, y)=3cos(xy)+x+ y^2 , where -4≤x≤4， -4≤y≤4 , please find the maximum value by the PSO algorithm. And plot the curve of the fitness evolution.

Answer 1

[Sam Chak]

Hi @Chihiro Ogino

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);

Local Maximum x at
    3.1614    4.0000

display('with the value of the fitness function = '); disp(Best_value(Gmax));

with the value of the fitness function = 
   22.1520

[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

[Chihiro Ogino]

Thank you so much! God bless!

[Sam Chak]

Good to hear that, Chihiro-san. If you find the code helpful, please consider accepting the Answer.