粒子群算法实战：基于MATLAB

求解$z=x^2+y^2-10cos(2\pi x)-10cos(2\pi y)-20$在x、y属于[-5,-5]的最大值

主函数 PSO.m

1 绘制函数图像

clear; clc
% 绘制函数图像
f1 = figure('color', [1, 1, 1]);
[x, y] = meshgrid(-5 : 0.1 : 5, -5 : 0.1 : 5);
z = x .^ 2 + y .^ 2 - 10 * cos(2 * pi * x) - 10 * cos(2 * pi * y) - 20;
mesh(x, y, z)

2 初始化参数

% 设置种群参数
maxgen = 200; % 最大运行代数
sizepop = 50; % 种群大小
nVar = 2; % 粒子搜索空间维数
% 省略的r1、r2的赋值不影响代码的运行
w = 1; % 惯性权重
c1 = 2; % 认知学习因子
c2 = 2; % 社会学习因子
popmax = 5; % 粒子最大位置(自变量的最大值)
popmin = -5; % 粒子最小位置(自变量的最小值)
Vmax = 1; % 粒子最大运行速度
Vmin = -1; % 粒子最小运行速度

pop = zeros(sizepop, nVar); % 粒子位置
V = zeros(sizepop, nVar); % 粒子速度
fitness = zeros(sizepop, 1); % 粒子适应度

3 初始化种群

for i = 1 : sizepop % 初始化种群(放鸟)
    % rand(0 - 1) rands(-1 - 1)
    pop(i, :) = (popmax - popmin) * rand(1, nVar) + popmin; % 初始化粒子位置
    V(i, :) = (Vmax - Vmin) * rand(1, nVar) + Vmin; % 初始化粒子速度
    fitness(i) = fun(pop(i, :));
end

[bestfitness, bestindex] = max(fitness);
gbest = pop(bestindex, :); % 全局最优个体
fitnessgbest = bestfitness; % 全局最优值 
pbest = pop; % 初始个体

fitnesspbest = fitness; % 初始个体最优值

4 主题函数：实现搜索过程

yy = zeros(maxgen, 1); % 每一代最优值
for i = 1 : maxgen
    % 更新粒子的速度和位置
    for j = i : sizepop
        V(j, :) = w .* V(j, :) + c1 * rand * (pbest(j, :) - pop(j, :)) + ...
                    c2 * rand * (gbest - pop(j, :)); % 速度更新公式
        if V(j, :) > Vmax
            V(j, :) = Vmax;
        elseif V(j, :) < Vmin
            V(j, :) = Vmin;
        end
        
        pop(j, :) = pop(j, :) + V(j, :); % 位置更新公式
        for n = 1 : nVar
            if pop(j, n) > popmax
                pop(j, n) = popmax;
            elseif pop(j, n) < popmin
                pop(j, n) = popmin;
            end
        end
     
        fitness(j) = fun(pop(j, :)); % 计算粒子目标函数值
    end

    for j = 1 : sizepop
        % 更新个体最优
       if fitness(j) > fitnesspbest(j)
          pbest(j, :) = pop(j, :);
          fitnesspbest(j) = fitness(j);
       end
       % 更新全局最优
       if fitness(j) > fitnessgbest
          gbest = pop(j, :);
          fitnessgbest = fitness(j);
       end
    end
    yy(i) = fitnessgbest;
end

5 结果展示

%% 输出结果
disp(['目标函数的最大值为：', num2str(fitnessgbest)])
disp(['自变量的取值为：', num2str(gbest)])
f2 = figure('color', [1, 1, 1]);
mesh(x, y, z)
hold on % 叠加绘图
plot3(gbest(1), gbest(2), fitnessgbest(1), 'ro', 'MarkerFaceColor', 'r', ...
        'MarkerEdgeColor', 'k', 'MarkerSize', 5)
title('粒子寻优结果')
xlabel('X轴'); ylabel('Y轴'); zlabel('Z轴')
f3 = figure('color', [1, 1, 1]);
plot(yy)
title('适应值变化曲线')
xlabel('迭代次数'); ylabel('适应值')

适应度函数 fun.m

function y = fun(x)
y = x(1) .^ 2 + x(2) .^ 2 - 10 * cos(2 * pi * x(1)) - 10 * cos(2 * pi * x(2)) - 20;