粒子群优化算法(Particle Swarm Optimization, PSO)
粒子群优化算法是进化计算的一个分支,是一种模拟自然界的生物活动 的随机搜索算法。模拟了自然界鸟群捕食和鱼群捕食的过程。通过群体中的协作寻找到问题的全局最优解。
文章目录
粒子更新公式
v i d = ω × v i d + c 1 × r 1 d × ( p B e s t i d − x i d ) + c 2 × r 2 d × ( g B e s t d − x i d ) x i d = x i d + v i d v_i^d=\omega\times v_i^d+c_1\times r_1^d\times(pBest_i^d-x_i^d)+c_2\times r_2^d\times(gBest^d-x_i^d)\\ x_i^d=x_i^d+v_i^d vid=ω×vid+c1×r1d×(pBestid−xid)+c2×r2d×(gBestd−xid)xid=xid+vid
设个体位置 X i = ( x i 1 , x i 2 , ⋯ , x i n ) , V i = ( v i 1 , v i 2 , ⋯ , v i n ) X_i=(x_i^1,x_i^2,\cdots,x_i^n),V_i=(v_i^1,v_i^2,\cdots,v_i^n) Xi=(xi1,xi2,⋯,xin),Vi=(vi1,vi2,⋯,vin)
- 下标 i i i——种群的第 i i i个个体
- 上标 d d d——个体中第 d d d个分量
- ω \omega ω——惯性权重,控制着前一速度对当前速度的影响, 用于平衡算法的探索和开发能力
- c 1 , c 2 c_1,c_2 c1,c2——加速系数,代表粒子向自身极值 p B e s t pBest pBest和全局极值 g B e s t gBest gBest推进的加速权值
- r 1 , r 2 r_1,r_2 r1,r2——随机数[0, 1],确保随机搜索能力
- p B e s t pBest pBest——自身极值坐标
- g B e s t gBest gBest——全局极值坐标
终止条件
- 设置最大循环次数或函数评估次数
- 求得可接受范围的解时结束
- 算法在很长一段迭代中没有任何改善时结束
全局和局部PSO
- 全局PSO速度快,但容易陷入局部最优
- 局部PSO收敛慢,但不易陷入局部最优
同步和异步更新
区别在于对全局 g B e s t gBest gBest和自身 p B e s t pBest pBest的更新方法
- 同步更新中,在每一代,当所有粒子都采取当前的 g B e s t gBest gBest进行速度和位置更新之后才对粒子进行评估,更新各自的 p B e s t pBest pBest,在选择最好的 g B e s t gBest gBest。
- 异步更新中,在每一代,粒子采用当前的 g B e s t gBest gBest进行速度和位置更新,然后马上评估,更新自己的 p B e s t pBest pBest,如果由于当前 g B e s t gBest gBest,则更新 g B e s t gBest gBest,迅速将更优的 g B e s t gBest gBest用于 后面的粒子更新。
一般而言,异步更新的PSO具高效的信息传播能 力,具有有更快的收敛速度。
算法流程
简单实现粒子群算法
使用Eigen库。
粒子类
// Particle.h
#pragma once
#include <Eigen/Eigen>
class Particle
{
public:
Eigen::VectorXf position;
Eigen::VectorXf velocity;
float score; // 当前评估值
Eigen::VectorXf pBest; // 最优解
float pBestScore; // 最优评估值
public:
Particle(int n, float xMin, float xMax);
Particle():score(std::numeric_limits<float>::max()), pBestScore(std::numeric_limits<float>::max()) {};
friend bool operator<(const Particle& p1, const Particle& p2) {
return p1.score < p2.score;
}
};
// Particle.cpp
#include "Particle.h"
#include <iostream>
#include <ctime>
#include <Eigen/Eigen>
#include <limits>
// n——粒子长度(x1...xn), xMin < x < xMax
Particle::Particle(int n, float xMin, float xMax)
:score(std::numeric_limits<float>::max())
{
position = Eigen::VectorXf(n);
velocity = Eigen::VectorXf(n);
for (int i = 0; i < n; i++) {
position[i] = rand() % static_cast<int>(xMax - xMin) - (xMax - xMin) / 2.0f;
velocity[i] = (rand() % static_cast<int>(xMax - xMin) - (xMax - xMin) / 2.0) * ((rand() % 20 + 1) / 100.0 + 0.2);
// 当前速度使解超出范围,重新随机
while (position[i] + velocity[i] > xMax || position[i] + velocity[i] < xMin)
{
velocity[i] = (rand() % static_cast<int>(xMax - xMin) - (xMax - xMin) / 2.0) * ((rand() % 20 + 1) / 100.0 + 0.2);
}
}
// 初始化局部最优
pBest = position;
pBestScore = score;
}
粒子群类
// ParticleGroup.h
#pragma once
#include <Eigen/Eigen>
#include <vector>
#include "Particle.h"
class ParticleGroup
{
public:
std::vector<Particle> Group; // 粒子群
Eigen::VectorXf gBest; // 全局最优解
float gBestScore; // 全局最优评估
private:
int groupSize; // 种群大小
int particleSize; // 粒子大小
float xMin;
float xMax;
float weight;
float c1, c2;
int Gmax; // 迭代次数
float (*ObjectFunc)(const Particle&); // 目标函数
public:
ParticleGroup(int N, int n, float x_min, float x_max, float w, float c1, float c2, int gmax, float (*func)(const Particle&));
void Optimize(); // 优化
private:
float Eval(const Particle& p); // 评估粒子
void UpdateParticle(Particle& p); // 更新粒子速度位置
void UpdateBest(); // 更新最优
inline float GetRandomC() // 随机[0,1]
{
return static_cast<float>(rand() % 101) / 100.0;
}
};
// ParticleGroup.cpp
#include "ParticleGroup.h"
#include <iostream>
ParticleGroup::ParticleGroup(int N, int n, float x_min, float x_max, float w, float _c1, float _c2, int gmax, float(*func)(const Particle&))
:groupSize(N),particleSize(n),xMin(x_min),xMax(x_max),weight(w),c1(_c1),c2(_c2),Gmax(gmax),ObjectFunc(func)
{
for (int i = 0; i < N; i++)
{
Particle p(n, x_min, x_max); // 随机粒子
p.score = Eval(p); // 评估
Group.push_back(p);
}
UpdateBest(); // 更新全局最优
}
float ParticleGroup::Eval(const Particle& p)
{
return ObjectFunc(p);
}
void ParticleGroup::UpdateParticle(Particle& p)
{
float r1 = GetRandomC();
float r2 = GetRandomC();
Eigen::VectorXf newVelocity = weight * p.velocity + c1 * r1 * (p.pBest - p.position) + c2 * r2 * (gBest - p.position);
// 更新位置速度
p.position = p.position + newVelocity;
for (int i=0;i<p.position.size();i++)
{
if (p.position[i] > xMax)
p.position[i] = xMax;
else if (p.position[i] < xMin)
p.position[i] = xMin;
}
p.velocity = newVelocity;
}
// 优化——异步更新
void ParticleGroup::Optimize()
{
for (int i = 0; i < Gmax; i++)
{
for (auto& p : Group)
{
UpdateParticle(p);
p.score = Eval(p);
// 更新个体最优
if (p.score < p.pBestScore)
{
p.pBest = p.position;
p.pBestScore = p.score;
}
// 更新全局最优
if (p.pBestScore < gBestScore)
{
gBest = p.pBest;
gBestScore = p.pBestScore;
}
}
std::cout << "[Step " << i + 1 << "] Score: " << gBestScore<< std::endl;
}
}
void ParticleGroup::UpdateBest()
{
Particle best;
for (auto p : Group)
{
if (p < best)
{
best = p;
}
}
gBest = best.position;
gBestScore = best.score;
}
main函数
优化两个连续函数,求最小值
f
1
=
∑
i
=
1
n
x
i
2
,
−
100
<
=
x
<
=
100
f
2
=
∑
i
=
1
n
[
x
i
2
−
10
cos
(
2
π
x
i
)
+
10
]
,
−
5.12
<
=
x
<
=
5.12
\begin{aligned} &f_1=\sum_{i=1}^nx_i^2\quad,-100<=x<=100\\ &f_2=\sum_{i=1}^n[x_i^2-10\cos(2\pi x_i)+10]\quad,-5.12<=x<=5.12 \end{aligned}
f1=i=1∑nxi2,−100<=x<=100f2=i=1∑n[xi2−10cos(2πxi)+10],−5.12<=x<=5.12
#include <iostream>
#include <cmath>
#include <ctime>
#include "Particle.h"
#include "ParticleGroup.h"
const float pi = 3.14159;
float ObjectFunc1(const Particle& p)
{
float sum = 0;
for (int i = 0; i < p.position.size(); i++)
{
sum += p.position[i] * p.position[i];
}
return sum;
}
float ObjectFunc2(const Particle& p)
{
float sum = 0;
for (int i = 0; i < p.position.size(); i++)
{
sum += p.position[i] * p.position[i] - 10 * cos(2 * pi * p.position[i]) + 10;
}
return sum;
}
int main()
{
srand(time(0));
int N = 20;
int n = 30;
float x_min = -100;
float x_max = 100;
float w = 0.2;
float c1 = 3.0;
float c2 = 3.0;
int Gmax = 2000;
float(*func)(const Particle&) = ObjectFunc1;
ParticleGroup G(N, n, x_min, x_max, w, c1, c2, Gmax, func);
G.Optimize();
std::cout << std::endl;
for (int i = 0; i < G.gBest.size(); i++)
{
std::cout << G.gBest[i];
if (i % 4 == 0 && i != 0)
std::cout << std::endl;
else std::cout << "\t";
}
}
实验结果
种群大小N = 20,粒子大小 n = 30,最大迭代Gmax=2000
对目标函数1:
f
1
=
∑
i
=
1
n
x
i
2
f_1=\sum_{i=1}^nx_i^2
f1=i=1∑nxi2
不同参数组合30次平均最优值,
$ c$ \ ω \omega ω | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.6 | 0.7 | 0.8 | 0.9 |
---|---|---|---|---|---|---|---|---|---|
0.5 | 4630.84 | 5773.61 | 5403.27 | 5227.46 | 4467.92 | 3208.11 | 2841.99 | 2704.21 | 2166.46 |
1.0 | 5445.2 | 4888.12 | 3751.11 | 2864.6 | 3029.65 | 2641.12 | 1919.3 | 1379.56 | 1446.35 |
1.5 | 6666.63 | 4796.14 | 3926.47 | 2571.03 | 2151.45 | 1380.62 | 1361.32 | 1963.52 | 1.54317 |
2.0 | 3066.41 | 776.496 | 646.178 | 1248.21 | 1616.63 | 1666.16 | 1026.29 | 5.40241 | 361.723 |
2.5 | 0.051045 | 0.0834625 | 0.147702 | 0.465709 | 0.554524 | 667.047 | 1.03938 | 30.6529 | 4420.42 |
3.0 | 0.0407009 | 0.0135592 | 0.0195125 | 0.0972921 | 0.76264 | 9.22703 | 595.156 | 3638.07 | 15697.8 |
3.5 | 0.15986 | 0.138149 | 0.3788 | 3.35501 | 31.7988 | 881.197 | 3541.52 | 15576.4 | 27478.4 |
可以看到参数组合 ω = 0.2 , c 1 = c 2 = 3.0 \omega=0.2,c_1=c_2=3.0 ω=0.2,c1=c2=3.0具有最好的效果,在该参数组合下迭代2000的30次运行结果如下:
运行次数 | 最优函数值 | 运行次数 | 最优函数值 |
---|---|---|---|
1 | 0.0157651 | 16 | 0.0155469 |
2 | 0.00023543 | 17 | 0.00468773 |
3 | 0.0444555 | 18 | 0.00939207 |
4 | 0.0654469 | 19 | 0.0172933 |
5 | 0.00541587 | 20 | 0.0847485 |
6 | 0.0251924 | 21 | 0.152291 |
7 | 0.00320896 | 22 | 0.158301 |
8 | 0.0115285 | 23 | 0.0184512 |
9 | 0.0276425 | 24 | 0.0120654 |
10 | 0.00361924 | 25 | 0.00193475 |
11 | 0.0126247 | 26 | 0.00755276 |
12 | 0.0030805 | 27 | 0.00507301 |
13 | 0.000822947 | 28 | 0.00929038 |
14 | 0.00999743 | 29 | 0.00145645 |
15 | 0.0173169 | 30 | 0.00806818 |
而对另一组参数 ω = 0.6 , c 1 = c 2 = 2.0 \omega=0.6,c_1=c_2=2.0 ω=0.6,c1=c2=2.0,在迭代15000次的运行结果如下:
运行次数 | 最优函数值 | 运行次数 | 最优函数值 |
---|---|---|---|
1 | 0.390762 | 16 | 0.237949 |
2 | 0.6735 | 17 | 0.563062 |
3 | 0.228124 | 18 | 0.269186 |
4 | 0.499667 | 19 | 0.628782 |
5 | 2.31534 | 20 | 0.996604 |
6 | 5.38956 | 21 | 5.84938 |
7 | 0.0977704 | 22 | 0.489698 |
8 | 1.55453 | 23 | 4.67167 |
9 | 0.609252 | 24 | 0.662991 |
10 | 1.36529 | 25 | 2.36874 |
11 | 0.386664 | 26 | 3.75168 |
12 | 0.26491 | 27 | 2.53652 |
13 | 1.55851 | 28 | 0.304391 |
14 | 1.37107 | 29 | 1.6134 |
15 | 1.34155 | 30 | 0.642401 |
可以看出最优参数组合不仅稳定,而且收敛速度很快,而 ω = 0.6 , c 1 = c 2 = 2.0 \omega=0.6,c_1=c_2=2.0 ω=0.6,c1=c2=2.0时要迭代15000次才能有少数几次得到较满意的结果,而且十分不稳定。
对于目标函数2:
f
2
=
∑
i
=
1
n
[
x
i
2
−
10
cos
(
2
π
x
i
)
+
10
]
f_2=\sum_{i=1}^n[x_i^2-10\cos(2\pi x_i)+10]
f2=i=1∑n[xi2−10cos(2πxi)+10]
不同参数组合30次平均最优值,
$ c$ \ ω \omega ω | 0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.6 | 0.7 | 0.8 | 0.9 |
---|---|---|---|---|---|---|---|---|---|
0.5 | 251.441 | 248.744 | 230.428 | 221.246 | 220.616 | 216.748 | 199.621 | 192.256 | 195.869 |
1.0 | 236.34 | 227.494 | 222.604 | 220.547 | 206.073 | 196.203 | 185.437 | 191.831 | 190.182 |
1.5 | 227.224 | 220.173 | 204.188 | 190.05 | 182.723 | 191.814 | 188.572 | 181.372 | 144.751 |
2.0 | 190.135 | 187.269 | 189.321 | 191.638 | 182.22 | 174.156 | 165.532 | 143.316 | 206.345 |
2.5 | 114.52 | 123.144 | 158.977 | 132.679 | 128.274 | 122.812 | 131.472 | 198.868 | 254.179 |
3.0 | 27.0794 | 51.3114 | 103.024 | 130.466 | 167.03 | 205.635 | 252.067 | 258.16 | 267.313 |
3.5 | 40.5218 | 41.5576 | 109.442 | 169.86 | 227.575 | 255.501 | 267.308 | 262.852 | 266.164 |
参数组合 ω = 0.1 , c 1 = c 2 = 3.0 \omega=0.1,c_1=c_2=3.0 ω=0.1,c1=c2=3.0时效果最好,在该组合下迭代2000次的30次运行如下:
运行次数 | 最优函数值 | 运行次数 | 最优函数值 |
---|---|---|---|
1 | 1.07258 | 16 | 28.9406 |
2 | 28.9494 | 17 | 55.812 |
3 | 0.0187387 | 18 | 78.8774 |
4 | 2.296 | 19 | 0.0182142 |
5 | 0.0502567 | 20 | 28.9255 |
6 | 0.0360365 | 21 | 24.9017 |
7 | 89.3454 | 22 | 28.9662 |
8 | 34.4231 | 23 | 29.0229 |
9 | 60.3701 | 24 | 57.8556 |
10 | 1.77543 | 25 | 35.0938 |
11 | 2.308 | 26 | 1.05813 |
12 | 37.9278 | 27 | 28.9334 |
13 | 57.8725 | 28 | 53.109 |
14 | 28.5563 | 29 | 24.8787 |
15 | 24.9157 | 30 | 24.8876 |
可以看出对目标函数2的求解容易陷入局部最优解,但也有几次运行能接近全局最优。
存在的问题
- 容易陷入局部最优
- 没有对更新过程中可能出现的越界进行处理