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
该处使用函数作为目标函数检验算法，在具体的问题中修改该函数为目标函数即可，注意该处为矩阵计算，使用.* .^等进行计算
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实现（一）
网友评论