Particle Swarm Optimization （粒子群优化算法）
一、算法简介
粒子群算法（ Particle Swarm Optimization, PSO）最早是由Eberhart和Kennedy于1995年提出，它的基本概念源于对鸟群觅食行为的研究，其基本思想在于通过群体中个体之间的协作和信息共享来寻找最优解。设想这样一个场景:一群鸟在随机搜寻食物，在这个区域里只有一块食物，所有的鸟都不知道食物在哪里，但是它们知道当前的位置离食物还有多远。那么找到食物的最优策略就是搜寻目前离食物最近的鸟的周围区域。
PSO的优势在于简单，容易实现，无需梯度信息，参数少，特别是其天然的实数编码特点特别适合于处理实优化问题。
粒子群算法当中的基本概念：
粒子：优化问题的候选解
位置：候选解所在的位置、速度：候选解移动的速度
适应度：评价粒子优劣的值，一般设置为目标函数值
个体最佳位置：单个粒子迄今为止找到的最佳位置、群体最佳位置：所有粒子迄今为止找到的最佳位置
二、算法原理
粒子群算法初始化为一群随机的粒子（随机解），然后根据迭代找到最优解。每一次迭代中，粒子通过跟踪两个极值来更新自己：
- 第1个是粒子本身所找到的最优解，这个称为个体极值
- 第2个是整个种群目前找到的最优解，这个称为全局极值（也可以只用种群其中的一部分作为粒子的邻居，称为局部极值）
整体流程：
-
**初始化粒子群：**为每个粒子随机分配初始位置和速度。这些位置和速度应在定义的问题空间范围内。初始化每个粒子的个体最佳位置（pbest）为其初始位置，全局最佳位置（gbest）为整个粒子群中的最佳位置。
-
**评估粒子性能：**对每个粒子，使用目标函数评估其当前位置的性能。如果当前位置比该粒子的pbest性能更好，则更新该粒子的pbest。如果当前位置比所有粒子的gbest性能更好，则更新gbest。
-
**更新速度和位置：**根据pbest、gbest、当前速度以及两个学习因子（个体和社会学习因子）来更新每个粒子的速度。 更新每个粒子的位置。
-
**重复或终止：**如果满足终止条件（达到最大迭代次数或gbest的改进小于阈值），则停止算法。
速度调整：
-
寻找个体极值：
如果它在pbestx的右边，则vx[]=vx[] - rand() * p_increment。
如果它在pbestx的左边，则vx[]=vx[] + rand() * p_increment。
-
寻找种群极值：
速度调整基于一个粗糙的不等式检验：如果presentx > bestx，则使其更小；如果presentx < bestx，那就把它变大。
一些实验表明，进一步修改算法使其更容易理解，并提高了其性能。不是简单地测试不平等的符号，而是根据最佳位置的每个维度的差异来调整速度：
没有很好的方法来猜测p_increment及g_increment是否应该更大，所以进行简化
简化版本：
标准版本速度调整公式：
位置调整公式：
三、测试代码
测试问题：
测试代码：
%% 粒子群算法PSO: 求解函数y = x1^2+x2^2-x1*x2-10*x1-4*x2+60在[-15,15]内的最小值
clear; clc
%% 绘制函数的图形
x1 = -15:1:15;
x2 = -15:1:15;
[x1,x2] = meshgrid(x1,x2);
y = x1.^2 + x2.^2 - x1.*x2 - 10*x1 - 4*x2 + 60;
mesh(x1,x2,y)
xlabel('x1'); ylabel('x2'); zlabel('y'); % 加上坐标轴的标签
axis vis3d % 冻结屏幕高宽比，使得一个三维对象的旋转不会改变坐标轴的刻度显示
hold on % 不关闭图形，继续在上面画图
%% 粒子群算法中的预设参数（参数的设置不是固定的，可以适当修改）
n = 30; % 粒子数量
narvs = 2; % 变量个数
c1 = 2; % 每个粒子的个体学习因子，也称为个体加速常数
c2 = 2; % 每个粒子的社会学习因子，也称为社会加速常数
w = 0.9; % 惯性权重
K = 100; % 迭代的次数
vmax = [6 6]; % 粒子的最大速度
x_lb = [-15 -15]; % x的下界
x_ub = [15 15]; % x的上界
%% 初始化粒子的位置和速度
x = zeros(n,narvs);
for i = 1: narvs
x(:,i) = x_lb(i) + (x_ub(i)-x_lb(i))*rand(n,1); % 随机初始化粒子所在的位置在定义域内
end
v = -vmax + 2*vmax .* rand(n,narvs); % 随机初始化粒子的速度（这里我们设置为[-vmax,vmax]）
%% 计算适应度(注意，因为是最小化问题，所以适应度越小越好)
fit = zeros(n,1); % 初始化这n个粒子的适应度全为0
for i = 1:n % 循环整个粒子群，计算每一个粒子的适应度
fit(i) = Obj_fun2(x(i,:)); % 调用Obj_fun2函数来计算适应度
end
pbest = x; % 初始化这n个粒子迄今为止找到的最佳位置（是一个n*narvs的向量）
ind = find(fit == min(fit), 1); % 找到适应度最小的那个粒子的下标
gbest = x(ind,:); % 定义所有粒子迄今为止找到的最佳位置（是一个1*narvs的向量）
%% 在图上标上这n个粒子的位置用于演示
h = scatter3(x(:,1),x(:,2),fit,'*r'); % scatter3是绘制三维散点图的函数（这里返回h是为了得到图形的句柄，未来我们对其位置进行更新）
%% 迭代K次来更新速度与位置
fitnessbest = ones(K,1); % 初始化每次迭代得到的最佳的适应度
for d = 1:K % 开始迭代，一共迭代K次
for i = 1:n % 依次更新第i个粒子的速度与位置
v(i,:) = w*v(i,:) + c1*rand(1)*(pbest(i,:) - x(i,:)) + c2*rand(1)*(gbest - x(i,:)); % 更新第i个粒子的速度
% 如果粒子的速度超过了最大速度限制，就对其进行调整
for j = 1: narvs
if v(i,j) < -vmax(j)
v(i,j) = -vmax(j);
elseif v(i,j) > vmax(j)
v(i,j) = vmax(j);
end
end
x(i,:) = x(i,:) + v(i,:); % 更新第i个粒子的位置
% 如果粒子的位置超出了定义域，就对其进行调整
for j = 1: narvs
if x(i,j) < x_lb(j)
x(i,j) = x_lb(j);
elseif x(i,j) > x_ub(j)
x(i,j) = x_ub(j);
end
end
fit(i) = Obj_fun2(x(i,:)); % 重新计算第i个粒子的适应度
if fit(i) < Obj_fun2(pbest(i,:)) % 如果第i个粒子的适应度小于这个粒子迄今为止找到的最佳位置对应的适应度
pbest(i,:) = x(i,:); % 那就更新第i个粒子迄今为止找到的最佳位置
end
if fit(i) < Obj_fun2(gbest) % 如果第i个粒子的适应度小于所有的粒子迄今为止找到的最佳位置对应的适应度
gbest = pbest(i,:); % 那就更新所有粒子迄今为止找到的最佳位置
end
end
fitnessbest(d) = Obj_fun2(gbest); % 更新第d次迭代得到的最佳的适应度
pause(0.1) % 暂停0.1s
h.XData = x(:,1); % 更新散点图句柄的x轴的数据（此时粒子的位置在图上发生了变化）
h.YData = x(:,2); % 更新散点图句柄的y轴的数据（此时粒子的位置在图上发生了变化）
h.ZData = fit; % 更新散点图句柄的z轴的数据（此时粒子的位置在图上发生了变化）
end
figure(2)
plot(fitnessbest) % 绘制出每次迭代最佳适应度的变化图
xlabel('迭代次数');
disp('最佳的位置是：'); disp(gbest)
disp('此时最优值是：'); disp(Obj_fun2(gbest))
function [y] = Obj_fun2(x)
y = x(1)^2 + x(2)^2 - x(1)*x(2) - 10*x(1) - 4*x(2) + 60;
end
测试结果：
最佳的位置是：7.9991 6.0162， 此时最优值是：8.0003
