智能计算——粒子群优化算法(PSO)(C++实现)

粒子群优化算法(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×(pBestidxid)+c2×r2d×(gBestdxid)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具高效的信息传播能 力,具有有更快的收敛速度。

算法流程

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=1nxi2,100<=x<=100f2=i=1n[xi210cos(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=1nxi2

不同参数组合30次平均最优值,

$ c$ \ ω \omega ω0.10.20.30.40.50.60.70.80.9
0.54630.845773.615403.275227.464467.923208.112841.992704.212166.46
1.05445.24888.123751.112864.63029.652641.121919.31379.561446.35
1.56666.634796.143926.472571.032151.451380.621361.321963.521.54317
2.03066.41776.496646.1781248.211616.631666.161026.295.40241361.723
2.50.0510450.08346250.1477020.4657090.554524667.0471.0393830.65294420.42
3.00.04070090.01355920.01951250.09729210.762649.22703595.1563638.0715697.8
3.50.159860.1381490.37883.3550131.7988881.1973541.5215576.427478.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次运行结果如下:

运行次数最优函数值运行次数最优函数值
10.0157651160.0155469
20.00023543170.00468773
30.0444555180.00939207
40.0654469190.0172933
50.00541587200.0847485
60.0251924210.152291
70.00320896220.158301
80.0115285230.0184512
90.0276425240.0120654
100.00361924250.00193475
110.0126247260.00755276
120.0030805270.00507301
130.000822947280.00929038
140.00999743290.00145645
150.0173169300.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次的运行结果如下:

运行次数最优函数值运行次数最优函数值
10.390762160.237949
20.6735170.563062
30.228124180.269186
40.499667190.628782
52.31534200.996604
65.38956215.84938
70.0977704220.489698
81.55453234.67167
90.609252240.662991
101.36529252.36874
110.386664263.75168
120.26491272.53652
131.55851280.304391
141.37107291.6134
151.34155300.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=1n[xi210cos(2πxi)+10]
不同参数组合30次平均最优值,

$ c$ \ ω \omega ω0.10.20.30.40.50.60.70.80.9
0.5251.441248.744230.428221.246220.616216.748199.621192.256195.869
1.0236.34227.494222.604220.547206.073196.203185.437191.831190.182
1.5227.224220.173204.188190.05182.723191.814188.572181.372144.751
2.0190.135187.269189.321191.638182.22174.156165.532143.316206.345
2.5114.52123.144158.977132.679128.274122.812131.472198.868254.179
3.027.079451.3114103.024130.466167.03205.635252.067258.16267.313
3.540.521841.5576109.442169.86227.575255.501267.308262.852266.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次运行如下:

运行次数最优函数值运行次数最优函数值
11.072581628.9406
228.94941755.812
30.01873871878.8774
42.296190.0182142
50.05025672028.9255
60.03603652124.9017
789.34542228.9662
834.42312329.0229
960.37012457.8556
101.775432535.0938
112.308261.05813
1237.92782728.9334
1357.87252853.109
1428.55632924.8787
1524.91573024.8876

可以看出对目标函数2的求解容易陷入局部最优解,但也有几次运行能接近全局最优。

存在的问题

  • 容易陷入局部最优
  • 没有对更新过程中可能出现的越界进行处理
评论 1
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

吉拉尔

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值