200字范文 > 粒子群优化算法（PSO）及其Matlab实现

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

时间：2022-01-28 09:53:58

github:/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×sin(x)+x×sin(2x)x\times sin(x)+x\times sin(2x)x×sin(x)+x×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输出
设定空间维数为1Tip
该参数指代自变量个数，1意味着这是个单目标优化问题

代码

全部实现代码如下：

clc;clear;close all;%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% Author: Myoontyee.Chen%% Data：1221%% 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:dx = 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 onplot(xm, f(xm), 'ro');title('参数初始化结果');%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 群体更新figure(2);iter = 1;record = zeros(maxIterations, 1); % 记录器while iter <= maxIterationsfx = f(x) ; % 个体当前适应度 for i = 1:InitialNumif fxm(i) < fx(i)fxm(i) = fx(i); % 更新个体历史最佳适应度xm(i,:) = x(i,:); % 更新个体历史最佳位置end endif 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)]);