粒子群算法介绍、matlab实现及相关改进

粒子群算法介绍、matlab实现及相关改进

1. 参数

   • N：粒子数量

     $⟹$一般取[20,40]，对于较难或者特定类别的问题，可以取[100,200];

   • D：决策变量维度

   • iter_Max：最大迭代次数

   • X：决策变量

   • V：每个粒子个体对应速度集合

   • $x_{\text{j,i}}^L$,$x_{\text{j,i}}^U$：决策变量 边界

   • $v_{\text{j,i}}^L$,$v_{\text{j,i}}^U$：速度边界

     $⟹$v_max较大时，探索能力强，但是粒子容易飞过最优解；

     $⟹$v_max较小时，开发能力强，但是容易陷入局部最优解；

     $⟹$v_max一般设为每一维变脸变化范围的10%~20%；

   • w：惯性权重

     $⟹$w越大，全局搜索能力越强；w越小，局部搜索能力越强。

     $⟹$根据经验，一般最大取0.9，最小取0.4

   • c1：学习因子1$⟹$个人认知学习因子，反映了粒子对自身历史经验的记忆和会议，代表粒子向自身历史最佳位置逼近的趋势。

   • c2：学习因子2$⟹$个人认知学习因子，反映了粒子建协同合作与知识共享的群体历史经验，代表粒子向群体或邻近历史最佳位置逼近的趋势。

     上述三个参数均可以采取不同的方式进行设置，在后面的算法改进中，会详细介绍。

   • r1、r2：两个位于0~1之间的随机数，用来增加系统的粒子运动趋势的随机性。

2. 算法流程

   1. 参数及粒子初始化

   2. 进入迭代更新粒子的位置和速度(可以按照标号4中的不同方法进行更新)

   3. 更新个体最优解和全局最优解

      在一迭代次数下，每更新一个粒子的位置和速度，便更新一次最优解(亲测，要是一次性更新的话，很难找到全局最优解)

   4. 完成迭代，对结果进行展示

3. 流程总结

4. 不同方法对于粒子群算法的改进

   1. 基础粒子群算法
      $$\begin{gathered} v_{\text{i,j}}(t+1)=v_{\text{i,j}}(t)+c_1{r}_1[p_{\text{i,best}}(t)-x_{\text{i,j}}(t)]+c_2{r}_2[g_{\text{i,best}}(t)-x_{\text{i,j}}(t)] \\ {x}_{\text{i,j}}(t+1)=x_{\text{i,j}}(t)+v_{\text{i,j}}(t+1) \end{gathered}$$

      • 优点：参数少，速度快，逻辑简单，编写代码方便
      • 缺点：全局搜索能力差，容易陷入局部最优解，搜索精度低

   2. 标准粒子群算法
      $$\begin{array}{l}\omega ={\omega }_{\text{max}}-({\omega }_{\text{max}}-{\omega }_{\text{min}})\frac{\text{iter}}{\text{iter\_Max}}\\ v_{\text{i,j}}(t+1)=\omega \cdot v_{\text{i,j}}(t)+c_1{r}_1[p_{\text{i,best}}(t)-x_{\text{i,j}}(t)]+c_2{r}_2[g_{\text{i,best}}(t)-x_{\text{i,j}}(t)]\\ x_{\text{i,j}}(t+1)=x_{\text{i,j}}(t)+v_{\text{i,j}}(t+1)\\ {\omega }_{\text{max}}=0.9;\\ {\omega }_{\text{max}}=0.4;\end{array}$$
      通过增加随迭代次数改变而改变的惯性因子，使得开始时具有较强的全局搜索能力，后期具备较强的局部搜索能力。

   3. 压缩因子粒子群算法
      $$\begin{array}{l}\omega =\frac{2}{|2-\phi -\sqrt{{\phi }^2-4\phi }|}\\ v_{\text{i,j}}(t+1)=\omega \cdot v_{\text{i,j}}(t)+c_1{r}_1[p_{\text{i,best}}(t)-x_{\text{i,j}}(t)]+c_2{r}_2[g_{\text{i,best}}(t)-x_{\text{i,j}}(t)]\\ x_{\text{i,j}}(t+1)=x_{\text{i,j}}(t)+v_{\text{i,j}}(t+1)\\ \phi =c_1+c_2\end{array}$$
      通过这种方式，收敛速度更快

   4. 结合自适应惯性权重的混合粒子群算法
      $$\begin{array}{l}\omega =({\omega }_{\text{max}}-{\omega }_{\text{min}})\cdot \tan (0.875(1-\frac{\text{iter}}{\text{iter\_Max}}{)}^k)+{\omega }_{\text{min}};\\ v_{\text{i,j}}(t+1)=\omega \cdot v_{\text{i,j}}(t)+c_1{r}_1[p_{\text{i,best}}(t)-x_{\text{i,j}}(t)]+c_2{r}_2[g_{\text{i,best}}(t)-x_{\text{i,j}}(t)];\\ x_{\text{i,j}}(t+1)=x_{\text{i,j}}(t)+v_{\text{i,j}}(t+1);\\ k=0.6;控制因子，根据经验，一般取该值\end{array}$$
      该算法将模拟退火中的Metropolis准则引入到了粒子群算法中，极大程度的避免寻优过程中陷入局部最优问题，而且收敛速度快，稳定性好。

   5. 动态调整惯性权重改进粒子群算法(IDWPSO)
      $$\begin{array}{l}\omega ={\omega }_{\text{min}}+({\omega }_{max}-{\omega }_{min})\cdot e^{-\frac{\text{iter}}{\text{iter\_Max}}}+\sigma \cdot betarand(p,q);\\ v_{\text{i,j}}(t+1)=\omega \cdot v_{\text{i,j}}(t)+c_1{r}_1[p_{\text{i,best}}(t)-x_{\text{i,j}}(t)]+c_2{r}_2[g_{\text{i,best}}(t)-x_{\text{i,j}}(t)];\\ x_{\text{i,j}}(t+1)=x_{\text{i,j}}(t)+v_{\text{i,j}}(t+1);\\ \sigma =0.1;惯性调整因子，根据经验，一般取该值\\ p=1;\\ q=3;\end{array}$$
      使得$\omega$随着迭代次数的增加而非线性减小，同时在最后面添加一个生成贝塔分布规律的随机数，使得算法在迭代后期可以增加算法的全局搜索能力，减少算法陷入局部最优解的可能性。

   6. 不同改进算法的横向性能对比

      为了方便比较不同算法在面对同一优化问题时性能表现，这里针对某一优化问题采用不同的方法进行500次求解，具体数据见下两表：

表1 500次运算结果展示

方法	迭代次数平均值	迭代次数方差	运行时间平均值	运行时间方差	最优适应度平均值	最优适应度方差	单次迭代时间
基础粒子群算法	93	3168	0.1476	0.00008326	-6.40186	0.00001752	0.0015871
标准粒子群	134	2178	0.1545	0.00003303	-6.40720	0.00000029	0.0011530
压缩因子粒子群算法	92	3294	0.1492	0.00004816	-6.40126	0.00002318	0.0016217
自适应惯性权重混合粒子群算法	146	2811	1.7072	0.0025975	-6.39763	0.0000888	0.0116931
IDWPSO	187	0	1.3836	0.00371457	-6.40611	0.00000191	0.0073989

P.S. 上述数据转载自智能算法之粒子群算法及改进

运行效率排序：
$$压缩因子粒子群算法>基础粒子群算法>标准粒子群>自适应惯性权重混合粒子群算法>IDWPSO$$
寻找最优解速度排序：
$$基础粒子群算法>压缩因子粒子群算法>标准粒子群>IDWPSO>自适应惯性权重混合粒子群算法$$
单次迭代时间排序：
$$标准粒子群>基础粒子群算法>压缩因子粒子群算法>IDWPSO>自适应惯性权重混合粒子群算法$$
寻找最优解能力：
$$自适应惯性权重混合粒子群算法>压缩因子粒子群算法>基础粒子群算法>IDWPSO>标准粒子群$$
综合以上考虑，在使用粒子群算法时：

• 时间要求不高：采用自适应惯性权重混合粒子群算法（能有效避免局部最优解）和IDWPSO；
• 时间要求比较高时：采用压缩因子粒子群算法或者标准粒子群算法

5. 示例代码

close all
clear
clc

N = 50; % 粒子个数
D = 10; % 决策变量维度
X_Max = 20; X_Min = -20; % 决策变量上下限
V_Max = (X_Max - X_Min) * 0.15; V_Min = -V_Max; % 速度上下限
w_Max = 0.9; w_Min = 0.4; % 惯性因子上下限
c1 = 1.5; c2 = 1.5; % 两个学习因子
iter_Max = 1000; % 最大迭代次数
WaitbarInter = iter_Max / 100; % 一个和进度条有关的参数

% 粒子初始化
X = X_Min + (X_Max - X_Min) .* rand(N, D);
V = V_Min + (V_Max - V_Min) .* rand(N, D);

% 写入当前全局最优解和个体最优解
GBestFitness = zeros(1, iter_Max + 1);
PBestFitness = zeros(N, 1);
PBest = X;

for i = 1:N
    PBestFitness(i) = Fitness(X(i, :));
end

[GBestFitness(1), GBestIndex] = min(PBestFitness);
GBest = PBest(GBestIndex, :);

tic
h = waitbar(0, ['已完成:0%   运算中...用时:', num2str(toc)]);

for iter = 1:iter_Max
    % 初始化当前迭代次数下的最优解
    CurrentIterGbest = GBest;
    CurrentIterGbestFitness = GBestFitness(iter);

    k = 0.6; % 控制因子，一般取0.6
    w = (w_Max - w_Min) * tan(0.875 * (1 - (iter / iter_Max)^k)) + w_Min;

    for i = 1:N
        % 更新个体的位置和速度
        V(i, :) = w * V(i, :) + c1 * rand * (PBest(i, :) - X(i, :)) + c2 * rand * (CurrentIterGbest - X(i, :));
        X(i, :) = X(i, :) + V(i, :);
        % 边界条件限制
        [X(i, :), V(i, :)] = BoundaryLimit(X(i, :), V(i, :), X_Max, X_Min, V_Max, V_Min);
        % 计算当前迭代次数下当前粒子的适应度
        CurrentFitness = Fitness(X(i, :));
        % 更新当前迭代次数下该粒子的最优解！！！
        % 事实证明，粒子个体最优解的更新应该在一个粒子位置苏都更新完后立即更新！！！
        % 若在当前迭代次数下所有粒子的位置和速度更新完毕后，再更新个体最优解，那么及其难以收敛！！！
        if CurrentFitness < PBestFitness(i)
            PBestFitness(i) = CurrentFitness;
            PBest(i, :) = X(i, :);
            % 更新当前迭代次数下的最优解
            if CurrentIterGbestFitness > CurrentFitness
                CurrentIterGbestFitness = CurrentFitness;
                CurrentIterGbest = X(i, :);
            end

        end

    end

    % 依据当前迭代次数下的最优解更新全局最优解
    GBest = CurrentIterGbest;
    GBestFitness(iter + 1) = CurrentIterGbestFitness;

    % 展示进度条
    if mod(iter, WaitbarInter) == 0
        waitbar(iter / iter_Max, h, ['已完成:' num2str(iter / iter_Max * 100) ...
        '%   运算中...用时:', num2str(toc),'/',num2str(toc/(iter / iter_Max))])
    end

end

close(h)

% 结果展示
semilogy(GBestFitness(1:iter - 1), 'Color', 'r', 'linewidth', 1.2)
title('Convergence Curve')
xlabel('Iteration');
ylabel('Best score obtained so far');
axis tight
grid on
box on
legend('PSO')
display(['The best solution obtained by PSO is : ', num2str(GBest)]);
display(['The best optimal value of the objective funciton found by PSO is : ', num2str(GBestFitness(iter - 1))]);

% 目标函数
function result = Fitness(X)
    result = sum(X.^2);

end

% 边界限制函数
function [result_X, result_V] = BoundaryLimit(X, V, X_Max, X_Min, V_Max, V_Min)

    for i_temp = 1:length(X(:, 1))

        if X(i_temp) > X_Max
            X(i_temp) = X_Max;
        end

        if X(i_temp) < X_Min
            X(i_temp) = X_Min;
        end

    end

    for i_temp = 1:length(V(:, 1))

        if V(i_temp) > V_Max
            V(i_temp) = V_Max;
        end

        if V(i_temp) < V_Min
            V(i_temp) = V_Min;
        end

    end

    result_X = X;
    result_V = V;

end

%% 关于粒子群算法的改进方法

%% 1.基础粒子群算法      基本不用，效果太差
% V(i, :) = V(i, :) + c1 * rand * (PBest(i, :) - X(i, :)) + c2 * rand * (CurrentIterGbest - X(i, :));
% X(i, :) = X(i, :) + V(i, :);

%% 2.标准粒子群算法    常用的粒子群更新方法
% w = w_Max - (w_Max - w_Min) * iter / iter_Max;
% V(i, :) = w * V(i, :) + c1 * rand * (PBest(i, :) - X(i, :)) + c2 * rand * (CurrentIterGbest - X(i, :));
% X(i, :) = X(i, :) + V(i, :);

%% 3.压缩因子粒子群算法   需要迭代的次数更多，全局搜索性更高
% phi = c1 + c2;  % 压缩因子粒子群算法相关参数
% lambda = 2/abs(2-phi-sqrt(phi^2-4*phi)); % 压缩因子粒子群算法相关参数
% lambda = 2/abs(2-phi-sqrt(phi^2-4*phi));
% V(i, :) = lambda * V(i, :) + c1 * rand * (PBest(i, :) - X(i, :)) + c2 * rand * (CurrentIterGbest - X(i, :));
% X(i, :) = X(i, :) + V(i, :);

%% 4.结合自适应惯性权重的混合粒子群算法(上面所使用的)
% 结合了模拟退火的思想(运用了Metropolis准则)，增加了全局搜索能力
% k = 0.6;  % 控制因子，一般取0.6
% w = (w_Max-w_Min)*tan(0.875*(1-(iter/iter_Max)^k))+w_Min;
% V(i, :) = w * V(i, :) + c1 * rand * (PBest(i, :) - X(i, :)) + c2 * rand * (CurrentIterGbest - X(i, :));
% X(i, :) = X(i, :) + V(i, :);

%% 5.动态调整惯性权重改进粒子群算法(IDWPSO)
% sigma = 0.1;  % 惯性调整因子
% p = 1;  q = 3;  % 两个beta分布的参数
% w = w_Min + (w_Max-w_Min)*exp(-iter/iter_Max)+sigma*betarnd(p,q);