粒子群优化算法（PSO）及其Matlab实现

github:https://github.com/Myoontyee/PSO

背景

粒子群优化（Particle Swarm Optimization, PSO），又称微粒群算法，是由J. Kennedy和R. C. Eberhart等于1995年开发的一种演化计算技术，来源于对一个简化社会模型的模拟。其中“群（swarm）”来源于微粒群匹配M. M. Millonas在开发应用于人工生命（artificial life）的模型时所提出的群体智能的5个基本原则。“粒子（particle）”是一个折衷的选择，因为既需要将群体中的成员描述为没有质量、没有体积的，同时也需要描述它的速度和加速状态。

PSO算法最初是为了图形化的模拟鸟群优美而不可预测的运动。而通过对动物社会行为的观察，发现在群体中对信息的社会共享提供一个演化的优势，并以此作为开发算法的基础。通过加入近邻的速度匹配、并考虑了多维搜索和根据距离的加速，形成了PSO的最初版本。之后引入了惯性权重来更好的控制开发（exploitation）和探索（exploration），形成了标准版本。为了提高粒群算法的性能和实用性，中山大学、（英国）格拉斯哥大学等又开发了自适应（Adaptive PSO）版本和离散（discrete）版本。

流程

标准PSO的算法流程如下：
  i)初始化一群微粒（群体规模为m），包括随机的位置和速度；
  ii)评价每个微粒的适应度；
  iii)对每个微粒，将它的适应值和它经历过的最好位置pbest的作比较，如果较好，则将其作为当前的最好位置pbest；
  iv)对每个微粒，将它的适应值和全局所经历最好位置gbest的作比较，如果较好，则重新设置gbest的索引号；
  v)变化微粒的速度和位置；
  vi)如未达到结束条件（通常为足够好的适应值或达到一个预设最大代数Gmax），回到ii）。

特点

  i)相较于传统算法计算速度非常快，全局搜索能力也很强；
  ii)PSO对于种群大小不十分敏感;
  iii)适用于连续函数极值问题，对于非线性、多峰问题均有较强的全局搜索能力。
  即

• 输入
  连续函数极值、非线性、多峰值问题
• f(x)
  PSO、MOPSO、etc…
• 输出
  全局较优解、全局最优解

变量

针对建议修改与建议默认参数均给予解释与建议值。

  1)待解目标函数

  程序：

f= @(x)x.*sin(x)+x.*sin(2.*x);          % 待解目标函数

  解释：

• 输入
  需要求解的目标函数（该程序仅针对单目标优化）
• 输出
  定义优化的目标函数
• Tip
  该处使用函数 $x\times sin(x)+x\times sin(2x)$作为目标函数检验算法，在具体的问题中修改该函数为目标函数即可，注意该处为矩阵计算，使用.* .^等进行计算

  2)待解函数上下限

  程序：

xLower = 0;                             % 待解函数下限
xTop = 30;                              % 待解函数上限

  解释：

• 输入
  目标函数上下限，2个函数定义域内的值
• 输出
  定义优化时自变量的上下限

  3)插值

  程序：

Interpolation = 0.01;                   % 插值

  解释：

• 输入
  小于函数上下限的数值
• 输出
  定义目标函数求解时的插值密度
• Tip
  理论上任意小于函数上下限的数值均可，一般取值≤0.01，越小的数值将使得目标函数连续性更好，求解时间更长

  4)最大迭代次数

  程序：

maxIterations = 100;                    % 最大迭代次数 

  解释：

• 输入
  最大迭代次数数值
• 输出
  定义最大迭代次数
• Tip
  取值范围一般为100~5000，依照实际情况，权衡计算时间与求解结果而定

  5)学习因子

  程序：

selfFactor = 3;                         % 自我学习因子
crowdFactor = 3;                        % 群体学习因子 

  解释：

• 输入
  个体（自我）与群体的学习因子
• 输出
  个体、群体学习因子定义
• Tip
  取值范围一般为0~4，依照实际情况，权衡自变量取值范围而定

  6)速度上下限

  程序：

vLower = -1;                            % 速度下限
vTop = 1;                               % 速度上限

  解释：

• 输入
  速度上下限数值
• 输出
  定义学习速度的上下限数值
• Tip
  i)过大速度将导致最优解被越过
  ii)过小速度将导致求解速度过慢

  7)初始种群个数
  程序：

InitialNum = 50;                        % 初始种群个数

  解释：

• 输入
  输入初始种群个数数值
• 输出
  定义初始化种群个数
• Tip
  取值范围一般为500~1000，PSO算法对种群大小不敏感，该处设定为50

  8)惯性权重
  程序：

weightFactor = 0.8;                     % 惯性权重

  解释：

• 输入
  惯性权重数值
• 输出
  定义惯性权重
• Tip
  取值范围一般为0.5~1，该参数反映个体历史成绩对现有成绩的影响

  9)空间维数
  程序：

d = 1;                                  % 空间维数

  解释：

• 输入
  1
• 输出
  设定空间维数为1
• Tip
  该参数指代自变量个数，1意味着这是个单目标优化问题

代码

全部实现代码如下：

clc;clear;close all;
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Author: Myoontyee.Chen
%% Data：20181221
%% License：BSD 3.0
%% PSO

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% Initialization初始化种群

%% 可修改参数

f= @(x)x.*sin(x)+x.*sin(2.*x);          % 待解目标函数
xLower = 0;                             % 待解函数下限
xTop = 30;                              % 待解函数上限
Interpolation = 0.01;                   % 插值
maxIterations = 100;                    % 最大迭代次数     
selfFactor = 3;                         % 自我学习因子
crowdFactor = 3;                        % 群体学习因子 
vLower = -1;                            % 速度下限
vTop = 1;                               % 速度上限
InitialNum = 50;                        % 初始种群个数
weightFactor = 0.8;                     % 惯性权重

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% 建议默认参数
d = 1;                                  % 空间维数
xLimit = [xLower, xTop];                % 位置参数限制
vLimit = [vLower, vTop];                % 设置速度限制
figure(1);
% 线型颜色
Figure1 = ezplot(f,[xLower, Interpolation, xTop]);
set(Figure1,'Color','k');


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 位置初始化及其图像

for i = 1:d
    x = xLimit(i, 1) + (xLimit(i, 2) - xLimit(i, 1)) * rand(InitialNum, d);
end                                     %初始种群的位置

vInitial = rand(InitialNum, d);         % 初始种群的速度
xm = x;                                 % 每个个体的历史最佳位置
ym = zeros(1, d);                       % 种群的历史最佳位置
fxm = zeros(InitialNum, 1);             % 每个个体的历史最佳适应度
fym = -inf;                             % 种群历史最佳适应度
hold on
plot(xm, f(xm), 'ro');title('参数初始化结果');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% 群体更新
figure(2);
iter = 1;
record = zeros(maxIterations, 1);       % 记录器
while iter <= maxIterations
     fx = f(x) ;                        % 个体当前适应度   
     for i = 1:InitialNum      
        if fxm(i) < fx(i)
            fxm(i) = fx(i);             % 更新个体历史最佳适应度
            xm(i,:) = x(i,:);           % 更新个体历史最佳位置
        end 
     end
if fym < max(fxm)
        [fym, nmax] = max(fxm);         % 更新群体历史最佳适应度
        ym = xm(nmax, :);               % 更新群体历史最佳位置
end                                     % 速度更新
    vInitial = vInitial * weightFactor + selfFactor * rand * (xm - x) + crowdFactor * rand * (repmat(ym, InitialNum, 1) - x);
                                        % 边界速度处理
    vInitial(vInitial > vLimit(2)) = vLimit(2);
    vInitial(vInitial < vLimit(1)) = vLimit(1);
    x = x + vInitial;                   % 位置更新
                                        % 边界位置处理
    x(x > xLimit(2)) = xLimit(2);
    x(x < xLimit(1)) = xLimit(1);
    record(iter) = fym;                 %最大值记录
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % 迭代过程及其图像
    % 若想直接观看结果，可注释掉该部分
    x0 = xLower : Interpolation : xTop;
    plot(x0, f(x0), 'k-', x, f(x), 'ro');title('状态位置变化');
    pause(0.1)
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    iter = iter+1;
end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 收敛过程记录图像
figure(3);
plot(record,'k-');
title('收敛过程记录');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% 最终状态位置图像
x0 = xLower : Interpolation : xTop;
figure(4);
plot(x0, f(x0), 'k-', x, f(x), 'ro');title('最终状态位置');
disp(['最大值：',num2str(fym)]);
disp(['变量取值：',num2str(ym)]);

结果

运行程序
输入

PSO.m

输出

最大值：45.8982
变量取值：26.0839

参考

[1]杨维, 李歧强. 粒子群优化算法综述[J]. 中国工程科学, 2004, 6(5):87-94.
[2]粒子群优化
[3]粒子群算法的matlab实现（一）