粒子群算法
粒子群算法基础
算法简介
粒子群算法，其全称为粒子群优化算法(ParticleSwarmOptimization,PSO)。它是通过模拟鸟群觅食行为而发展起来的一种基于群体协作的搜索算法。
算法策略
粒子群算法的目标是使所有粒子在多维超体中找到最优解。
首先给空间中的所有粒子分配初始随机位置和初始随机速度。然后根据每个粒子的速度、已知空间中的最优全局位置和粒子最优位置依次推进每个粒子的位置。随着计算的推移迭代，通过探索和利用搜索空间中已知的有利位置，粒子将聚集在一个或多个最优点附近。
它的核心思想是利用群体中的个体对信息的共享使整个群体的运动在问题求解空间中产生从无序到有序的演化过程，从而获得问题的可行解。
算法流程图
名词阐释
粒子：优化问题的候选解
位置：候选解所在的位置
速度：候选解移动的速度
适应度：评价粒子优劣的值，一般设置为目标函数值
个体最佳位置：单个粒子迄今为止找到的最佳位置
群体最佳位置：所有粒子迄今为止找到的最佳位置
符号说明
| 符号 | 含义 |
|---|---|
| n | 粒子个数 |
| c1 | 粒子的个体学习因子，也称为个体加速因子 |
| c2 | 粒子的社会学习因子，也称为社会加速因子 |
| w | 速度的惯性权重 |
| vid | 第d次迭代时，第i个粒子的速度 |
| xid | 第d次迭代时，第i个粒子所在的位置 |
| f(x) | 在位置x时的适应度值（一般取目标函数值） |
| pbestid | 到第d次迭代为止，第i个粒子经过的最好的位置 |
| gbestd | 到第d次迭代为止，所有粒子经过的最好的位置 |
图像解释
核心公式
这只鸟第d步的速度=上一步自身的速度惯性+ 自我认知部分+社会认知部分
νid=wνid-1+c1 r1 (pbestid- xid )+c2 r2 (gbestd- xid )
这只鸟第d+1步所在的位置= 第d步所在的位置+ 第d步的速度* 运动的时间
xid+1=xid+νid
参数分析
- 学习因子[1]
νid=wνid-1+c1 r1 (pbestid- xid )+c2 r2 (gbestd- xid )
学习因子c1和c2代表每个粒子推向 pi和pg的统计权重。低的值允许粒子在被拉回之前可以在目标区域外徘徊，而高的值则导致粒子突然冲向或者越过目标区域。
在初提出粒子群算法的论文中指出，个体学习因子和社会(或群体) 学习因子取2比较合适。
（注意：最初提出粒子群算法的这篇论文没有惯性权重）
- 惯性权重[2]
νid=wνid-1+c1 r1 (pbestid- xid )+c2 r2 (gbestd- xid )
惯性权重使粒子保持运动惯性，使其有扩展搜索空间的趋势，并有能力搜索新的区域。
论文中得到的结论：惯性权重取0.9‐1.2是比较合适的
- 位置更新方程各部分
- νid=wνid-1+c1 r1 (pbestid- xid )+c2 r2 (gbestd- xid )
如果没有第一部分，即 w = 0，则速度只取决于粒子当前位置和历史最好位置 pi和pg，速度本身没有记忆性。假设一个粒子位于全局最好位置，它将保持静止。而其他粒子则飞向它本身最好位置pi和全局最好位置pg的加权中心。在这种情况下，粒子群将收敛到当前的全局最好位置，更像一个局部最优算法。在加上第一部分后，粒子有扩展搜索空间的趋势，即第一部分有全局搜索能力。 - νid=wνid-1+c1 r1 (pbestid- xid )+c2 r2 (gbestd- xid )
如果没有第二部分，即 c1 = 0，则粒子没有认知能力。在粒子的相互作用下，有能力达到新的搜索空间。他的收敛速度比标准版本要快，但对复杂问题，则比标准版本更容易陷入局部最优点。 - νid=wνid-1+c1 r1 (pbestid- xid )+c2 r2 (gbestd- xid )
如果没有第三部分，即 c2 = 0,则粒子间没有社会共享信息能力。个体间没有交互，一个规模为M的群体等价于M个单个粒子的运行，因而得到最优解的机率很小。
Matlab实现
求解函数y = 7cos(5x) + 4*sin(x)在[-5,5]内的最大值
%% 绘制函数的图形
x = -5:0.01:5;
y = 7*cos(5*x) + 4*sin(x);
figure(1)
plot(x,y,'b-')
title('y = 7*cos(5*x) + 4*sin(x)')
hold on % 不关闭图形，继续在上面画图
%% 粒子群算法中的预设参数（参数的设置不是固定的，可以适当修改）
n = 10; % 粒子数量
narvs = 1; % 变量个数
c1 = 2; % 每个粒子的个体学习因子
c2 = 2; % 每个粒子的社会学习因子
w = 0.9; % 惯性权重
K = 50; % 迭代的次数
vmax = 1.2; % 粒子的最大速度
x_lb = -3; % x的下界
x_ub = 3; % x的上界
%% 初始化粒子的位置和速度
x = zeros(n,narvs);
x = x_lb + (x_ub-x_lb).*rand(n,narvs) % 随机初始化粒子所在的位置
v = -vmax + 2*vmax .* rand(n,narvs); % 随机初始化粒子的速度
%% 计算适应度
fit = zeros(n,1); % 初始化这n个粒子的适应度全为0
for i = 1:n % 循环整个粒子群，计算每一个粒子的适应度
fit(i) = Obj_fun1(x(i,:)); % 调用Obj_fun1函数来计算适应度
end
pbest = x; % 初始化这n个粒子迄今为止找到的最佳位置（是一个n*narvs的向量）
ind = find(fit == max(fit), 1); % 找到适应度最大的那个粒子的下标
gbest = x(ind,:); % 定义所有粒子迄今为止找到的最佳位置（是一个1*narvs的向量）
%% 在图上标上这n个粒子的位置用于演示
h = scatter(x,fit,80,'*r'); % scatter是绘制二维散点图的函数
%% 迭代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_fun1(x(i,:)); % 重新计算第i个粒子的适应度
if fit(i) > Obj_fun1(pbest(i,:)) % 如果第i个粒子的适应度大于这个粒子迄今为止找到的最佳位置对应的适应度
pbest(i,:) = x(i,:); % 更新第i个粒子迄今为止找到的最佳位置
end
if fit(i) > Obj_fun1(gbest) % 如果第i个粒子的适应度大于所有的粒子迄今为止找到的最佳位置对应的适应度
gbest = pbest(i,:); % 更新所有粒子迄今为止找到的最佳位置
end
end
fitnessbest(d) = Obj_fun1(gbest); % 更新第d次迭代得到的最佳的适应度
pause(0.1) % 暂停0.1s
h.XData = x; % 更新散点图句柄的x轴的数据
h.YData = fit; % 更新散点图句柄的y轴的数据
end
figure(2)
plot(fitnessbest) % 绘制出每次迭代最佳适应度的变化图
xlabel('迭代次数');
disp('最佳的位置是：'); disp(gbest)
disp('此时最优值是：'); disp(Obj_fun1(gbest))
function y = Obj_fun1(x)
y = 7*cos(5*x) + 4*sin(x);
end
求解结果
粒子群算法改进--线性递减惯性权重
惯性权重回顾
惯性权重w体现的是粒子继承先前的速度的能力
Shi,Y最先将惯性权重w引入到粒子群算法中,并且分析指出
-
一个较大的惯性权值有利于全局搜索
-
一个较小的权值则更利于局部搜索
线性递减惯性权重
- 在搜索初期，增强全局搜索能力可以更大可能遍历解空间，避免陷入局部最优解 >广撒网
- 在搜索后期，增强局部搜索能力可以更大可能的锁定最优解 >精准打击
为了更好地平衡算法的全局搜索以及局部搜索能力，Shi,Y提出了线性递减惯性权重LDIW(LinearDecreasingInertiaWeight)
公式如下
wd=wstart-(wstart - wend)x (d/K)
d是当前迭代的次数,K是迭代总次数
wstart一般取0.9，wend一般取0.4
与原来的相比，现在惯性权重和迭代次数有关
范例代码
%% 线性递减惯性权重的粒子群算法PSO: 求解函数y = 7*cos(5*x) + 4*sin(x)在[-5,5]内的最大值
clear; clc
%% 绘制函数的图形
x = -5:0.01:5;
y = 7*cos(5*x) + 4*sin(x);
figure(1)
plot(x,y,'b-')
title('y = 7*cos(5*x) + 4*sin(x)')
hold on % 不关闭图形，继续在上面画图
%% 粒子群算法中的预设参数（参数的设置不是固定的，可以适当修改）
n = 10; % 粒子数量
narvs = 1; % 变量个数
c1 = 2; % 每个粒子的个体学习因子
c2 = 2; % 每个粒子的社会学习因子
w_start = 0.9; % 初始惯性权重
w_end = 0.4; % 粒子群最大迭代次数时的惯性权重
K = 50; % 迭代的次数
vmax = 1.2; % 粒子的最大速度
x_lb = -3; % x的下界
x_ub = 3; % x的上界
%% 初始化粒子的位置和速度
x = zeros(n,narvs);
x = x_lb + (x_ub-x_lb).*rand(n,narvs) % 随机初始化粒子所在的位置
v = -vmax + 2*vmax .* rand(n,narvs); % 随机初始化粒子的速度
%% 计算适应度
fit = zeros(n,1); % 初始化这n个粒子的适应度全为0
for i = 1:n % 循环整个粒子群，计算每一个粒子的适应度
fit(i) = Obj_fun1(x(i,:)); % 调用Obj_fun1函数来计算适应度
end
pbest = x; % 初始化这n个粒子迄今为止找到的最佳位置（是一个n*narvs的向量）
ind = find(fit == max(fit), 1); % 找到适应度最大的那个粒子的下标
gbest = x(ind,:); % 定义所有粒子迄今为止找到的最佳位置（是一个1*narvs的向量）
%% 在图上标上这n个粒子的位置用于演示
h = scatter(x,fit,80,'*r'); % scatter是绘制二维散点图的函数
%% 迭代K次来更新速度与位置
fitnessbest = ones(K,1); % 初始化每次迭代得到的最佳的适应度
for d = 1:K % 开始迭代，一共迭代K次
for i = 1:n % 依次更新第i个粒子的速度与位置 **注意这边的变化**
w = w_start - (w_start - w_end) * d / K; % 更新惯性权重
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_fun1(x(i,:)); % 重新计算第i个粒子的适应度
if fit(i) > Obj_fun1(pbest(i,:)) % 如果第i个粒子的适应度大于这个粒子迄今为止找到的最佳位置对应的适应度
pbest(i,:) = x(i,:); % 更新第i个粒子迄今为止找到的最佳位置
end
if fit(i) > Obj_fun1(gbest) % 如果第i个粒子的适应度大于所有的粒子迄今为止找到的最佳位置对应的适应度
gbest = pbest(i,:); % 更新所有粒子迄今为止找到的最佳位置
end
end
fitnessbest(d) = Obj_fun1(gbest); % 更新第d次迭代得到的最佳的适应度
pause(0.1) % 暂停0.1s
h.XData = x; % 更新散点图句柄的x轴的数据
h.YData = fit; % 更新散点图句柄的y轴的数据
end
figure(2)
plot(fitnessbest) % 绘制出每次迭代最佳适应度的变化图
xlabel('迭代次数');
disp('最佳的位置是：'); disp(gbest)
disp('此时最优值是：'); disp(Obj_fun1(gbest))
function y = Obj_fun1(x)
y = 7*cos(5*x) + 4*sin(x);
end
运行结果如下
对比未经优化的粒子群算法
可知：求得的最优解更接近实际最优解
参考论文：Shi, Y. and Eberhart, R.C. (1999) Empirical Study of Particle Swarm Optimization. Proceedings of the 1999 Congress on Evolutionary Computation, Washington DC, 6‐9 July 1999, 1945‐1950.
粒子群算法改进--自适应惯性权重
适应度回顾
适应度用于评价粒子优劣，一般设置为目标函数值
众所周知
- 一个较大的惯性权重有利于全局搜索
- 一个较小的惯性权重有利于局部搜索
基本粒子群算法的粒子速度迭代公式：
vid=wvid-1+c1r1(pbestid-xid)+c2r2(gbestd-xid)
w是定值，代表粒子自身的惯性权重，但随着迭代次数的增加，问题的求解细节也会有所改变，固定值在整体求解的过程中存在不少缺陷。因而，引入变动的惯性权重，以动态适应问题的求解流程。以下是两种用自适应惯性权重求解问题的范例。
求解最小值问题
1、wmin和wmax是预设的最小与最大惯性系数，一般wmin取0.4，wmax取0.9
2、faveraged为第d次迭代时所有粒子的平均适应度
3、fmind=min{f(x1d),f(x2d),…，f(xnd)}，即第d次迭代时所有粒子的最小适应度
- 适应度越小，说明距离最优解越近，此时更需要局部搜索
- 适应度越大，说明距离最优解越远，此时跟需要全局搜索
求解最大值问题
1、wmin和wmax是预设的最小与最大惯性系数，一般wmin取0.4，wmax取0.9
2、faveraged为第d次迭代时所有粒子的平均适应度
3、fmaxd=max{f(x1d),f(x2d),…，f(xnd)}，即第d次迭代时所有粒子的最大适应度
- 适应度越小，说明距离最优解越近，此时更需要局部搜索
- 适应度越大，说明距离最优解越远，此时跟需要全局搜索
与基本粒子群算法相比，现在的惯性权重与迭代次数和每个粒子的适应度都有关
举个栗子
- 求解最大值问题
假设目前有五个粒子ABCDE
它们此刻各自的适应度分别为9，8，7，6，5
取最小惯性权重0.4，最大惯性权重0.9
那么，带入公式可以求得
faveraged=（9+8+7+6+5）/5 = 7
fmaxd = max{9,8,7,6,5} = 9
此刻五个粒子的惯性权重分别为：
A: 0.4+ (0.9-0.4) (9-9)/(9-7) = 0.4
B: 0.4+(0.9-0.4) (9-8)/(9-7) = 0.65
C: 0.4 + (0.9-0.4) (9-7)/(9-7) = 0.9
D: 0.9
E: 0.9
怎么去理解这个公式呢？
因为要求解的是最大值问题（求解最大适应度），那么当粒子的适应度小于平均适应度时，表示粒子距最大值相对较远，需要扩大搜索范围来寻找最值；当粒子的适应度大于平均适应度时，表示粒子距离最大值较近，需要缩小搜索范围进行局部精确搜索。
代码实现
clear
%% 绘制函数的图形
x = -5:0.01:5;
y = 7*cos(5*x) + 4*sin(x);
figure(1)
plot(x,y,'b-')
title('y = 7*cos(5*x) + 4*sin(x)')
hold on % 不关闭图形，继续在上面画图
%% 粒子群算法中的预设参数（参数的设置不是固定的，可以适当修改）
n = 10; % 粒子数量
narvs = 1; % 变量个数
c1 = 2; % 每个粒子的个体学习因子
c2 = 2; % 每个粒子的社会学习因子
w_max = 0.9; % 最大惯性权重，通常取0.9
w_min = 0.4; % 最小惯性权重，通常取0.4
K = 50; % 迭代的次数
vmax = 1.2; % 粒子的最大速度
x_lb = -3; % x的下界
x_ub = 3; % x的上界
%% 初始化粒子的位置和速度
x = zeros(n,narvs);
x = x_lb + (x_ub-x_lb).*rand(n,narvs) % 随机初始化粒子所在的位置
v = -vmax + 2*vmax .* rand(n,narvs); % 随机初始化粒子的速度
%% 计算适应度
fit = zeros(n,1); % 初始化这n个粒子的适应度全为0
for i = 1:n % 循环整个粒子群，计算每一个粒子的适应度
fit(i) = Obj_fun1(x(i,:)); % 调用Obj_fun1函数来计算适应度
end
pbest = x; % 初始化这n个粒子迄今为止找到的最佳位置（是一个n*narvs的向量）
ind = find(fit == max(fit), 1); % 找到适应度最大的那个粒子的下标
gbest = x(ind,:); % 定义所有粒子迄今为止找到的最佳位置（是一个1*narvs的向量）
%% 在图上标上这n个粒子的位置用于演示
h = scatter(x,fit,80,'*r'); % scatter是绘制二维散点图的函数
%% 迭代K次来更新速度与位置
fitnessbest = ones(K,1); % 初始化每次迭代得到的最佳的适应度
for d = 1:K % 开始迭代，一共迭代K次
for i = 1:n % 依次更新第i个粒子的速度与位置
f_i = fit(i); % 取出第i个粒子的适应度
f_avg = sum(fit)/n; % 计算此时适应度的平均值
f_max = max(fit); % 计算此时适应度的最大值
if f_i >= f_avg
if f_avg ~= f_max % 如果分母为0，我们就干脆令w=w_max
w = w_min + (w_max - w_min)*(f_max - f_i)/(f_max - f_avg);
else
w = w_max;
end
else
w = w_max;
end
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_fun1(x(i,:)); % 重新计算第i个粒子的适应度
if fit(i) > Obj_fun1(pbest(i,:)) % 如果第i个粒子的适应度大于这个粒子迄今为止找到的最佳位置对应的适应度
pbest(i,:) = x(i,:); % 更新第i个粒子迄今为止找到的最佳位置
end
if fit(i) > Obj_fun1(gbest) % 如果第i个粒子的适应度大于所有的粒子迄今为止找到的最佳位置对应的适应度
gbest = pbest(i,:); % 更新所有粒子迄今为止找到的最佳位置
end
end
fitnessbest(d) = Obj_fun1(gbest); % 更新第d次迭代得到的最佳的适应度
pause(0.1) % 暂停0.1s
h.XData = x; % 更新散点图句柄的x轴的数据
h.YData = fit; % 更新散点图句柄的y轴的数据
end
figure(2)
plot(fitnessbest) % 绘制出每次迭代最佳适应度的变化图
xlabel('迭代次数');
disp('最佳的位置是：'); disp(gbest)
disp('此时最优值是：'); disp(Obj_fun1(gbest))
function y = Obj_fun1(x)
y = 7*cos(5*x) + 4*sin(x);
end
运行结果如下
对比未经优化的粒子群算法
可知：求得的最优解更接近实际最优解
粒子群算法改进--随机惯性权重
随机惯性权重
随机惯性权重通过随机分布的方式获取惯性权重，这样既可以避免在迭代前期局部搜索能力的不足； 也可以避免在迭代后期全局搜索能力的不足，其在收敛速度和全局收敛性方面相比基本粒子群算法有明显提高。
随机惯性权重公式
vid = wvid-1 + c1r1(pbestid-xid) + c2r2(gbestd-xid)
w = μmin + (μmax - μmin) · rand() + σ · randn()
μmin是随机惯性权重的最小值；
μmax是随机惯性权重的最大值；
rand()为[0,1]均匀分布的随机数；
randn()为正态分布的随机数；
σ(标准差)用来度量随机惯性权重w与其数学期望之间的偏离程度（一般取0.2~0.5之间的一个数），该项是为了控制取值中的权重误差，使权重w有利于向期望权重方向进化（一般情况下实验误差服从正态分布）。
范例代码
求解函数y = 7cos(5x) + 4*sin(x)在[-5,5]内的最大值
%% 绘制函数的图形
x = -5:0.01:5;
y = 7*cos(5*x) + 4*sin(x);
figure(1)
plot(x,y,'b-')
title('y = 7*cos(5*x) + 4*sin(x)')
hold on % 不关闭图形，继续在上面画图
%% 粒子群算法中的预设参数（参数的设置不是固定的，可以适当修改）
n = 10; % 粒子数量
narvs = 1; % 变量个数
c1 = 2; % 每个粒子的个体学习因子
c2 = 2; % 每个粒子的社会学习因子
w_max = 0.9; % 最大惯性权重，通常取0.9
w_min = 0.4; % 最小惯性权重，通常取0.4
sigma = 0.3; % 正态分布的随机扰动项的标准差（一般取0.2-0.5之间）
K = 50; % 迭代的次数
vmax = 1.2; % 粒子的最大速度
x_lb = -3; % x的下界
x_ub = 3; % x的上界
%% 初始化粒子的位置和速度
x = zeros(n,narvs);
x = x_lb + (x_ub-x_lb).*rand(n,narvs) % 随机初始化粒子所在的位置
v = -vmax + 2*vmax .* rand(n,narvs); % 随机初始化粒子的速度
%% 计算适应度
fit = zeros(n,1); % 初始化这n个粒子的适应度全为0
for i = 1:n % 循环整个粒子群，计算每一个粒子的适应度
fit(i) = Obj_fun1(x(i,:)); % 调用Obj_fun1函数来计算适应度
end
pbest = x; % 初始化这n个粒子迄今为止找到的最佳位置（是一个n*narvs的向量）
ind = find(fit == max(fit), 1); % 找到适应度最大的那个粒子的下标
gbest = x(ind,:); % 定义所有粒子迄今为止找到的最佳位置（是一个1*narvs的向量）
%% 在图上标上这n个粒子的位置用于演示
h = scatter(x,fit,80,'*r'); % scatter是绘制二维散点图的函数
%% 迭代K次来更新速度与位置
fitnessbest = ones(K,1); % 初始化每次迭代得到的最佳的适应度
for d = 1:K % 开始迭代，一共迭代K次
for i = 1:n % 依次更新第i个粒子的速度与位置
w = w_min + (w_max - w_min)*rand(1) + sigma*randn(1);
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_fun1(x(i,:)); % 重新计算第i个粒子的适应度
if fit(i) > Obj_fun1(pbest(i,:)) % 如果第i个粒子的适应度大于这个粒子迄今为止找到的最佳位置对应的适应度
pbest(i,:) = x(i,:); % 更新第i个粒子迄今为止找到的最佳位置
end
if fit(i) > Obj_fun1(gbest) % 如果第i个粒子的适应度大于所有的粒子迄今为止找到的最佳位置对应的适应度
gbest = pbest(i,:); % 更新所有粒子迄今为止找到的最佳位置
end
end
fitnessbest(d) = Obj_fun1(gbest); % 更新第d次迭代得到的最佳的适应度
pause(0.1) % 暂停0.1s
h.XData = x; % 更新散点图句柄的x轴的数据
h.YData = fit; % 更新散点图句柄的y轴的数据
end
figure(2)
plot(fitnessbest) % 绘制出每次迭代最佳适应度的变化图
xlabel('迭代次数');
disp('最佳的位置是：'); disp(gbest)
disp('此时最优值是：'); disp(Obj_fun1(gbest))
function y = Obj_fun1(x)
y = 7*cos(5*x) + 4*sin(x);
end
求解结果
对比未经优化的粒子群算法
可知：求得的最优解更接近实际最优解
参考文献：基于随机惯性权重的简化粒子群优化算法[J]. 计算机应用研究, 2014, 031(002):361-363,391.
粒子群算法改进--压缩因子法
前言概述
粒子速度更新公式如下：
vid = wvid-1 + c1r1(pbestid-xid)+ c2r2(gbestd-xid)
在研究完粒子群算法中有关惯性权重的优化之后，我们把目光转向速度更新公式的后两项，根据之前所学可知：个体学习因子c1和社会学习因子c2决定了粒子本身经验信息和其他粒子的经验信息对粒子运行轨迹的影响，其反映了粒子群之间的信息交流。
较大的c1值，会使粒子过多地在自身的局部范围内搜索
较大的c2值，则又会促使粒子过早收敛到局部最优值。
为了有效地控制粒子的飞行速度，使算法达到全局搜索与局部搜索两者间的有效平衡，Clerc构造了引入收缩因子的PSO模型，采用了压缩因子，这种调整方法通过合适选取参数，可确保PSO算法的收敛性，并可取消对速度的边界限制。
参数设置
- 个体学习因子c1 = 2.05 （应用较多）
- 社会学习因子 c2 = 2.05（应用较多）
- 惯性权重 w = 0.9
- C = c1+c2 = 4.1
- 收缩因子
综上，速度更新公式如下
示例代码
求解函数y = 7cos(5x) + 4*sin(x)在[-5,5]内的最大值
…… ……
%% 初始化参数
c1 = 2.05; % 每个粒子的个体学习因子
c2 = 2.05; % 每个粒子的社会学习因子
C = c1 + c2;
w = 0.9; % 惯性权重
fai = 2/abs((2-C-sqrt(C^2-4*C)));% 收缩因子
%% 初始化粒子的位置和速度
…… ……
%% 计算适应度
…… ……
%% 迭代K次来更新速度与位置
fitnessbest = ones(K,1); % 初始化每次迭代得到的最佳的适应度
for d = 1:K % 开始迭代，一共迭代K次
for i = 1:n % 依次更新第i个粒子的速度与位置
v(i,:) = fai * (w*v(i,:) + c1*rand(1)*(pbest(i,:) - x(i,:)) + c2*rand(1)*(gbest - x(i,:))); % 更新第i个粒子的速度
x(i,:) = x(i,:) + v(i,:); % 更新第i个粒子的位置
% 如果粒子的位置超出了定义域，就对其进行调整
…… ……
…… ……
end
…… ……
function y = Obj_fun1(x)
y = 7*cos(5*x) + 4*sin(x);
end
运行结果如下
对比未经优化的粒子群算法
可知，求得的最优解更接近实际最优解
参考文献：Eberhart R C . Comparing inertia weights and
constriction factors in optimization[C]//
Proceedings of the 2000 IEEE Congress on
Evolutionary Computation, La Jolla, CA. IEEE, 2000.
粒子群算法改进--自动退出迭代循环
前言
- 当粒子已经找到最佳位置后，再增加迭代次数只会浪费计算时间，那么我们能否设计一个策略，能够自动退出迭代呢?
循环跳出策略
（1）初始化最大迭代次数、计数器以及最大计数值（例如分别取100, 0, 20）
（2）定义“函数变化量容忍度”，一般取非常小的正数；
（3）在迭代的过程中，每次计算出来最佳适应度后，都计算该适应度和上一次迭代时最佳适应度的变化量(取绝对值)；
（4）判断这个变化量和“函数变化量容忍度”的相对大小，如果前者小，则计数器加1；否则计数器清0；
（5）不断重复这个过程，有以下两种可能:
① 此时还没有超过最大迭代次数，计数器的值超过了最大计数值，那么我们就跳出迭代循环，搜索结束。
② 此时已经达到了最大迭代次数，那么直接跳出循环，搜索结束。
代码实现
求解函数y = 7cos(5x) + 4*sin(x)在[-5,5]内的最大值
%% 改进：自动退出迭代循环
Count = 0; % 计数器初始化为0
max_Count = 30; % 最大计数值初始化为30
tolerance = 1e-2; % 函数变化量容忍度，取10^(-2）
%% 绘制函数的图形
x = -5:0.01:5;
y = 7*cos(5*x) + 4*sin(x);
figure(1)
plot(x,y,'b-')
title('y = 7*cos(5*x) + 4*sin(x)')
hold on % 不关闭图形，继续在上面画图
%% 粒子群算法中的预设参数（参数的设置不是固定的，可以适当修改）
n = 10; % 粒子数量
narvs = 1; % 变量个数
c1 = 2; % 每个粒子的个体学习因子
c2 = 2; % 每个粒子的社会学习因子
w = 0.9; % 惯性权重
K = 100; % 迭代的次数
vmax = 1.2; % 粒子的最大速度
x_lb = -3; % x的下界
x_ub = 3; % x的上界
%% 初始化粒子的位置和速度
x = zeros(n,narvs);
x = x_lb + (x_ub-x_lb).*rand(n,narvs) % 随机初始化粒子所在的位置
v = -vmax + 2*vmax .* rand(n,narvs); % 随机初始化粒子的速度
%% 计算适应度
fit = zeros(n,1); % 初始化这n个粒子的适应度全为0
for i = 1:n % 循环整个粒子群，计算每一个粒子的适应度
fit(i) = Obj_fun1(x(i,:)); % 调用Obj_fun1函数来计算适应度
end
pbest = x; % 初始化这n个粒子迄今为止找到的最佳位置（是一个n*narvs的向量）
ind = find(fit == max(fit), 1); % 找到适应度最大的那个粒子的下标
gbest = x(ind,:); % 定义所有粒子迄今为止找到的最佳位置（是一个1*narvs的向量）
%% 在图上标上这n个粒子的位置用于演示
h = scatter(x,fit,80,'*r'); % scatter是绘制二维散点图的函数
%% 迭代K次来更新速度与位置
% fitnessbest = ones(K,1); % 初始化每次迭代得到的最佳的适应度
for d = 1:K % 开始迭代，一共迭代K次
tem = gbest; % 将上一步找到的最佳位置保存为临时变量
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_fun1(x(i,:)); % 重新计算第i个粒子的适应度
if fit(i) > Obj_fun1(pbest(i,:)) % 如果第i个粒子的适应度大于这个粒子迄今为止找到的最佳位置对应的适应度
pbest(i,:) = x(i,:); % 更新第i个粒子迄今为止找到的最佳位置
end
if fit(i) > Obj_fun1(gbest) % 如果第i个粒子的适应度大于所有的粒子迄今为止找到的最佳位置对应的适应度
gbest = pbest(i,:); % 更新所有粒子迄今为止找到的最佳位置
end
end
fitnessbest(d) = Obj_fun1(gbest); % 更新第d次迭代得到的最佳的适应度
delta_fit = abs(Obj_fun1(gbest) - Obj_fun1(tem)); % 计算相邻两次迭代适应度的变化量
if delta_fit < tolerance % 判断这个变化量和“函数变化量容忍度”的相对大小，如果前者小，则计数器加1
Count = Count + 1;
else
Count = 0; % 否则计数器清0
end
if Count > max_Count % 如果计数器的值达到了最大计数值
break; % 跳出循环
end
pause(0.1) % 暂停0.1s
h.XData = x; % 更新散点图句柄的x轴的数据
h.YData = fit; % 更新散点图句柄的y轴的数据
end
figure(2)
plot(fitnessbest) % 绘制出每次迭代最佳适应度的变化图
xlabel('迭代次数');
disp('最佳的位置是：'); disp(gbest)
disp('此时最优值是：'); disp(Obj_fun1(gbest))
function y = Obj_fun1(x)
y = 7*cos(5*x) + 4*sin(x);
end
迭代一百次的情况下优化与未优化的对比如下
自动退出迭代循环
基本粒子群算法
对比可知，未经优化的粒子群算法会浪费大量时间在无明显变化的搜索中。
浙公网安备 33010602011771号