粒子群算法(PSO)的C++实现

粒子群算法(PSO----Particle Swarm Optimization)是常用的智能算法之一，它模拟了 鸟群觅食 行为，是一种具有随机性的 仿生算法 。PSO算法在无约束条件函数最优化问题上具有全局搜索能力强，局部收敛能力较强的优点。

本篇博文目的在于：

记录基本的PSO算法原理。
利用C++将PSO抽象算法本身加以实现，构建适合一类函数优化问题的PSO算法类。并附上相应的源码。

PSO算法基本原理

对于一个简单的 无约束最优化 问题：
max ⁡ f ( x ) , ∀ x ∈ R d i m e n s i o n \max\quad f(x)\quad ,\forall x\in R^{dimension} maxf(x),∀x∈Rdimension
首先引入粒子群算法中粒子的结构概念。

struct Partical
{double _position[_dimension];  //粒子所在位置，即为问题候选解double _velocity[_dimension];    //粒子运动速度double _fitness;        //粒子适应度double _bestPosition[_dimension]; //粒子所经历过的最优位置
}；

一个粒子携带了一个候选解_position以及该候选解对应的函数值_fitness（这里称为粒子适应度）。另外粒子本身是运动的，具有速度_dimension。同时，一个粒子在运动过程中会保留自己所经历过的候选解中最优的一个_bestPosition。

接下来引入粒子群的结构概念：

struct pso
{Partical _particalSet[_particalCount]; //一群粒子构成的集合double _localBestGuideCoe; //个体最优解引导加速度double _globalBestGuideCoe; //全局最优解引导加速度double _globalBestPosition; //当前粒子群进化中的全局最优解double _globalBestFitness; //当前粒子群进化中的全局最优解适应度
}

粒子群简单来说就是由一群粒子_particalSet构成。但是粒子群还要规定粒子运动的一些属性以及记录这群粒子的一些状态。

模拟鸟类觅食的行为特性，将粒子想象成觅食鸟群中的一只鸟，这只鸟既想向着找到食物（适应度）最多的鸟那(_globalBestPosition)飞，又不想放弃自己以前找到过最多食物的那个位置(_bestPosition)。

数学一点进行描述，就是 粒子在全局最优方向和个体局部最优方向都具有一定的加速度 ，而系数_globalBestGuideCoe和_localBestGuideCoe就在一定程度上描述了这两个加速度的大小。
用数学公式进行表达，粒子群在第 k 时刻到第 k+1 时刻的运动描述如下：
V i ( k + 1 ) = V i ( k ) + G r g [ G B P ( k ) − P i ( k ) ] + L r l [ L B P i ( k ) − P i ( k ) ] P i ( k + 1 ) = P i ( k ) + V i ( k ) ∀ i = 1 , 2 , 3 , . . . , p a r t i c a l C o u n t V_i(k+1)=V_i(k)+Gr_g[GBP(k)-P_i(k)]+Lr_l[LBP_i(k)-P_i(k)]\\P_i(k+1)=P_i(k)+V_i(k)\\ \forall i = 1,2,3,...,particalCount Vi(k+1)=Vi(k)+Grg[GBP(k)−Pi(k)]+Lrl[LBPi(k)−Pi(k)]Pi(k+1)=Pi(k)+Vi(k)∀i=1,2,3,...,particalCount
更新策略中，V代表粒子速度，P为粒子所在位置，G为全局最优引导加速度_globalBestGuideCoe，L为局部最优引导加速度_localBestGuideCoe，GBP(k)代表 k 时刻粒子群的全局最优_globalBestPosition，LBP(k)代表k时刻该粒子的个体局部最优，也即_localBestPosition。r 代表0-1均匀分布的随机数，i为粒子在粒子群中的下标序号。

对于PSO算法，不难想象，在一定的条件下，全局的粒子从四面八方向着最优解聚拢过去，从而实现函数的最优化。

以下示意图是采用pso算法函数最大值点的运行过程。其中红点代表当前粒子群搜索得到的全局最优解。黑色箭头代表了粒子位置及其运动速度。待优化函数为：
z = s i n ( x 2 + y 2 ) / x 2 + y 2 z = sin(\sqrt{x^2+y^2})/\sqrt{x^2+y^2} z=sin(x2+y2 )/x2+y2
运行效果：

PSO算法的C++实现以及测试

本节提供并介绍PSO算法的C++代码，并附上算法的一个测试用例以说明算法类的使用方法。

与算法原理介绍一节类似，先介绍PSO粒子类：

class ZPSO_Partical
{public:int _dimension; //粒子维度double* _position; //粒子所在位置，即为问题候选解double* _velocity; //粒子运动速度double _fitness;  //粒子适应度double* _bestPosition; //粒子所经历过的最优位置void initial(int dimension); //粒子初始化，用以分配粒子所需内存void copy(ZPSO_Partical); //复制其他粒子的数据，类似赋值操作ZPSO_Partical(void); //构造~ZPSO_Partical(void); //析构
}；

代码注释对粒子类的各属性与函数都已经描述得比较全面，不再赘述。一下介绍粒子群算法类—— ZPSO_Algorithm 的几个重要函数。

/***************************************************************
* 函数名：ZPSO_Algorithm
* 函数描述：构造一个PSO算法
* 输入参数：
*  objFunction：优化目标函数指针
*  positionMinValue：解下界，其任一维度注意应低于上界
*  positionMaxValue：解上界，其任一维度注意应高于下界
*  dimension：解空间维度
*  particalCount：种群粒子个数
*  globalGuideCoe：粒子种群全局最优引导速度因子，默认为2
*  localGuideCoe：粒子种群个体最优引导速度因子，默认为2
*  maxSpeed：粒子运动最大速度
* 输出参数：
*  ZPSO_Algorithm&：构建得到的PSO算法本身
***************************************************************/
ZPSO_Algorithm(double (*objFunction)(ZPSO_Partical&),double *positionMinValue,double *positionMaxValue,int dimension,int particalCount,double globalGuideCoe = 2,double localGuideCoe = 2,double maxSpeed = 1)；

以上为PSO算法类 ZPSO_Algorithm的构造函数，每个参数的意义在注释中都有说明，其中要重点提到的是参数objFunction。

double objFunction(ZPSO_Partical& partical);

objFunction 是以粒子为输入，返回粒子适应度的一个函数。一般而言，在 objFunction 中读取 partical 的候选解_position，并调用待优化的目标函数计算_position位置的函数值作为粒子适应度返回。

以下介绍pso的优化函数：

/***************************************************************
* 函数名：findMax
* 函数描述：采用粒子群算法搜索最优解
* 输入参数：
*  times：粒子群进化次数
*  bestPartical：进化得到的最优个体
*  disturbanceRate：粒子速度扰动概率，默认为0.2
*  disturbanceVelocityCoe：速度扰动因子，表征扰动速度相对_maxSpeed大小
*                          用于扰动粒子速度以提高局部搜索能力，默认为0.05
* 输出参数：void
***************************************************************/
void findMax(int times,ZPSO_Partical& bestPartical,double disturbanceRate = 0.2,double disturbanceVelocityCoe = 0.05);

ZPSO_Algorithm 提供的最大化函数 findMax 需要用户给定粒子群进化次数（进化时间）times。需要说明的是，为了提高PSO算法的搜索能力，避免加速度选择要求过于严格的问题，本文实现的PSO算法中引入了速度扰动的概念，即粒子除了全局最优加速度与局部最优加速度，还有一个随机加速度。参数 disturbanceRate 和 disturbanceVelocityCoe则分别描述了产生随机加速度的概率与随机加速度的大小。令 disturbanceRate = 0 则算法退化为最基本的PSO。

以下提供测试样例，举例说明ZPSO_Algorithm的使用方法。

主函数文件 <main.cpp>：用于构建 objFunction 与调用ZPSO_Algorithm。

#include "ZPSOAlgorithm.h"
#include <math.h>
#include <iostream>
using namespace std;//待优化的目标函数
double objFunction(ZPSO_Partical& partical);int main(void)
{//定义求解范围,x在[-10,10]范围，y在[-10,10]范围double minPos[] = {-10,-10};double maxPos[] = {10,10};//定义问题描述参数int dimension = 2; //候选解维度int particalCount = 200; //粒子群粒子个数double globalGuideCoe = 1; //全局引导加速度double localGuideCoe = 1; //局部引导加速度double maxSpeed = 4; //粒子最大速度//构建pso算法ZPSO_Algorithm pso(objFunction,minPos,maxPos,dimension,particalCount,globalGuideCoe,localGuideCoe,maxSpeed);//运行pso算法ZPSO_Partical bestPartical; //粒子群最终进化结果int generation = 200; //粒子群进化代数pso.findMax(generation,bestPartical); //获取最终进化结果//输出最终结果cout<<"the best position for the objFunction is:"<<endl;cout<<"x="<<bestPartical._position[0]<<endl;cout<<"y="<<bestPartical._position[1]<<endl;cout<<"the best fitness for the objFunction is:"<<bestPartical._fitness<<endl;return(0);
}//优化目标函数定义
double objFunction(ZPSO_Partical &partical)
{//从partical中读取候选解double x = partical._position[0]-0.3;double y = partical._position[1]+0.1;//计算候选解对应的函数值double r = sqrt(x*x+y*y);double rtn;if(r<1e-8)rtn = 1;elsertn = sin(r)/r;return(rtn);
}

PSO算法类定义头文件 <ZPSOAlgorithm.h>

/************************************************************* 头文件名：ZPSOAlgorithm* 头文件描述：用户自定义粒子群算法* 作者：chikuo* 日期：20190706* 版本号：ZPSOAlgorithm_2019070601
************************************************************//************************************************************* 修改日期：20190708* 修改人：chikuo* 修改信息：引入扰动速度，提高pso算法局部搜索能力
************************************************************/
#ifndef _ZPSOALGORITHM_H
#define _ZPSOALGORITHM_H#include <stdlib.h>
#include <time.h>
#include <math.h>//粒子群算法例子个体
class ZPSO_Partical
{public:int _dimension; //粒子维度double *_position; //粒子所在位置数组指针double *_velocity; //例子速度double *_bestPosition; //当前粒子搜索得到过的最优位置double _fitness; //例子适应度double _bestFitness; //当前粒子搜索得到过的最优适应度//构造函数，粒子维度初始化为0ZPSO_Partical(void){_dimension = 0;}//析构函数，释放粒子内存~ZPSO_Partical(void){if(_dimension){delete []_position;delete []_velocity;delete []_bestPosition;}}//初始化函数，用于为粒子开辟内存空间void initial(int dimension){if(_dimension!=dimension && dimension){//需要重新分配内存if(_dimension){//消除已有内存delete []_position;delete []_velocity;delete []_bestPosition;}//开辟新内存_dimension = dimension;_position = new double[_dimension];_velocity = new double[_dimension];_bestPosition = new double[_dimension];}}//复制函数，用于粒子间的复制操作void copy(ZPSO_Partical& partical){this->initial(partical._dimension);for(int i=0;i<_dimension;i++){_position[i] = partical._position[i];_velocity[i] = partical._velocity[i];_bestPosition[i] = partical._bestPosition[i];_fitness = partical._fitness;_bestFitness = partical._bestFitness;}}
};//PSO算法
class ZPSO_Algorithm
{public:int _dimension; //粒子群维度int _particalCount; //种群粒子数量double _globalGuideCoe; //全局最优引导系数double _localGuideCoe; //局部最优引导系数ZPSO_Partical _globalBestPartical; //搜索过程得到的全局最优粒子double *_positionMinValue; //粒子位置的最小界double *_positionMaxValue; //粒子位置的最大界double _maxSpeed; //粒子允许最大速度double (*_fitnessFunction)(ZPSO_Partical&); //粒子适应度函数ZPSO_Partical *_particalSet; //粒子集/**************************************************************** 函数名：ZPSO_Algorithm* 函数描述：构造一个PSO算法* 输入参数：*  objFunction：优化目标函数指针*  positionMinValue：解下界，其任一维度注意应低于上界*  positionMaxValue：解上界，其任一维度注意应高于下界*  dimension：解空间维度*  particalCount：种群粒子个数*  globalGuideCoe：粒子种群全局最优引导速度因子，默认为2*  localGuideCoe：粒子种群个体最优引导速度因子，默认为2*  maxSpeed：粒子运动最大速度* 输出参数：*  ZPSO_Algorithm&：构建得到的PSO算法本身***************************************************************/ZPSO_Algorithm(double (*objFunction)(ZPSO_Partical&),double *positionMinValue,double *positionMaxValue,int dimension,int particalCount,double globalGuideCoe = 2,double localGuideCoe = 2,double maxSpeed = 1){//初始化类内参数并分配内存_fitnessFunction = objFunction;_dimension = dimension;_positionMinValue = new double[_dimension];_positionMaxValue = new double[_dimension];for(int i=0;i<_dimension;i++){_positionMinValue[i] = positionMinValue[i];_positionMaxValue[i] = positionMaxValue[i];}_particalCount = particalCount;_globalGuideCoe = globalGuideCoe;_localGuideCoe = localGuideCoe;_maxSpeed = maxSpeed;_particalSet = new ZPSO_Partical[_particalCount];for(int i=0;i<_particalCount;i++)_particalSet[i].initial(_dimension);_globalBestPartical.initial(_dimension);//配置随机数种子srand((unsigned int)time(NULL));}/**************************************************************** 函数名：~ZPSO_Algorithm* 函数描述：析构一个PSO算法，释放算法内存* 输入参数：void* 输出参数：void***************************************************************/~ZPSO_Algorithm(void){//释放内存delete []_positionMinValue;delete []_positionMaxValue;delete []_particalSet;}/**************************************************************** 函数名：rand0_1* 函数描述：生成一个0-1的均匀分布随机数* 输入参数：void* 输出参数：*  double：在0-1上均匀分布的随机数***************************************************************/double rand0_1(void){return((1.0*rand())/RAND_MAX);}/**************************************************************** 函数名：refresh* 函数描述：计算粒子适应度并更新粒子个体最优位置与全局最优位置* 输入参数：void* 输出参数：void***************************************************************/void refresh(void){int globalBestParticalIndex = -1;for(int i=0;i<_particalCount;i++){//循环遍历所有粒子，更新粒子适应度_particalSet[i]._fitness = this->_fitnessFunction(_particalSet[i]);if(_particalSet[i]._fitness > _particalSet[i]._bestFitness){//更新粒子的个体最优位置for(int j=0;j<_dimension;j++)_particalSet[i]._bestPosition[j] = _particalSet[i]._position[j];_particalSet[i]._bestFitness = _particalSet[i]._fitness;//是否更新全局最优解if(_particalSet[i]._bestFitness > _globalBestPartical._bestFitness)globalBestParticalIndex = i;}}//更新全局最优粒子位置if(globalBestParticalIndex != -1)_globalBestPartical.copy(_particalSet[globalBestParticalIndex]);}/**************************************************************** 函数名：randomlyInitial* 函数描述：随机初始化种群中粒子位置* 输入参数：void* 输出参数：void***************************************************************/void randomlyInitial(void){int globalBestParticalIndex = -1;//遍历所有粒子//初始化第0个粒子与全局最优粒子//初始化粒子位置与速度double velocityMod = 0;//遍历粒子的任一维度for(int j=0;j<_particalSet[0]._dimension;j++){//随机初始化粒子位置与最佳位置double tempVal = _positionMinValue[j];tempVal += rand0_1()*(_positionMaxValue[j]-_positionMinValue[j]);_particalSet[0]._position[j] = tempVal;_particalSet[0]._bestPosition[j] = tempVal;//随机初始化粒子速度_particalSet[0]._velocity[j] = rand0_1();velocityMod += _particalSet[0]._velocity[j]*_particalSet[0]._velocity[j];}//粒子速度归化为随机大小v_moddouble v_mod = rand0_1()*_maxSpeed;velocityMod = sqrt(velocityMod);for(int j=0;j<_particalSet[0]._dimension;j++)_particalSet[0]._velocity[j] *= (v_mod/velocityMod);//更新粒子初代适应度值与最佳适应度值_particalSet[0]._fitness = _fitnessFunction(_particalSet[0]);_particalSet[0]._bestFitness = _particalSet[0]._fitness;_globalBestPartical.copy(_particalSet[0]);//初始化1~_particalCount-1个粒子for(int i=1;i<_particalCount;i++){velocityMod = 0;//初始化粒子位置与速度//遍历粒子的任一维度for(int j=0;j<_particalSet[i]._dimension;j++){//随机初始化粒子位置与最佳位置double tempVal = _positionMinValue[j];tempVal += rand0_1()*(_positionMaxValue[j]-_positionMinValue[j]);_particalSet[i]._position[j] = tempVal;_particalSet[i]._bestPosition[j] = tempVal;//随机初始化粒子速度_particalSet[i]._velocity[j] = rand0_1();velocityMod += _particalSet[i]._velocity[j]*_particalSet[i]._velocity[j];}//粒子速度归化为随机大小v_modv_mod = rand0_1()*_maxSpeed;velocityMod = sqrt(velocityMod);for(int j=0;j<_particalSet[i]._dimension;j++)_particalSet[i]._velocity[j] *= (v_mod/velocityMod);//更新粒子初代适应度值与最佳适应度值_particalSet[i]._fitness = _fitnessFunction(_particalSet[i]);_particalSet[i]._bestFitness = _particalSet[i]._fitness;if(_particalSet[i]._bestFitness > _globalBestPartical._bestFitness)globalBestParticalIndex = i;}//更新粒子群全局最佳数据if(globalBestParticalIndex != -1)_globalBestPartical.copy(_particalSet[globalBestParticalIndex]);}/**************************************************************** 函数名：disturbance* 函数描述：对粒子速度进行给定大小的扰动* 输入参数：*  partical：被扰动的粒子对象*  relativeVelocityRate：扰动速度大小上限相对于_maxSpeed的比例，默认为0.05* 输出参数：void***************************************************************/void disturbance(ZPSO_Partical &partical,double relativeVelocityRate = 0.05){//生成扰动速度double *disturbanceVelocity = new double[_dimension];//随机生成扰动速度大小double disturbanceVelocityMod = relativeVelocityRate*_maxSpeed*rand0_1();double v_mod = 0;for(int i=0;i<_dimension;i++){disturbanceVelocity[i] = rand0_1();v_mod += disturbanceVelocity[i]*disturbanceVelocity[i];}v_mod = sqrt(v_mod);//扰动速度大小归化到disturbanceVelocityModfor(int i=0;i<_dimension;i++)disturbanceVelocity[i] *= (disturbanceVelocityMod/v_mod);//扰动粒子速度v_mod = 0;for(int i=0;i<_dimension;i++){partical._velocity[i] += disturbanceVelocity[i];v_mod += partical._velocity[i]*partical._velocity[i];}v_mod = sqrt(v_mod);//粒子速度受限if(v_mod > _maxSpeed)for(int i=0;i<_dimension;i++)partical._velocity[i] *= (_maxSpeed/v_mod);delete[]disturbanceVelocity;}/**************************************************************** 函数名：update* 函数描述：更新粒子群粒子位置与适应度* 输入参数：*  disturbanceRate：粒子速度扰动概率，默认为0.2*  disturbanceVelocityCoe：速度扰动因子，表征扰动速度相对_maxSpeed大小*                          用于扰动粒子速度以提高局部搜索能力，默认为0.05* 输出参数：void***************************************************************/void update(double disturbanceRate = 0.2,double disturbanceVelocityCoe = 0.05){double v_mod;//遍历所有粒子for(int i=0;i<_particalCount;i++){//遍历所有维度v_mod = 0;double r1 = rand0_1();double r2 = rand0_1();for(int j=0;j<_particalSet[i]._dimension;j++){//速度更新//全局最优位置加速度_particalSet[i]._velocity[j] += _globalGuideCoe*r1*(_globalBestPartical._bestPosition[j]-_particalSet[i]._position[j]);//个体局部最优位置加速度_particalSet[i]._velocity[j] += _localGuideCoe*r2*(_particalSet[i]._bestPosition[j]-_particalSet[i]._position[j]);//粒子速度模二v_mod += _particalSet[i]._velocity[j]*_particalSet[i]._velocity[j];}//粒子速度受限v_mod = sqrt(v_mod);if(v_mod > _maxSpeed)for(int j=0;j<_particalSet[i]._dimension;j++)_particalSet[i]._velocity[j] *= (_maxSpeed/v_mod);//对粒子速度进行扰动，提高算法局部搜索能力if(rand0_1()<disturbanceRate)this->disturbance(_particalSet[i],disturbanceVelocityCoe);//位置更新for(int j=0;j<_particalSet[i]._dimension;j++){_particalSet[i]._position[j] += _particalSet[i]._velocity[j];//粒子位置受限if(_particalSet[i]._position[j] < _positionMinValue[j])_particalSet[i]._position[j] = _positionMinValue[j];else if(_particalSet[i]._position[j] > _positionMaxValue[j])_particalSet[i]._position[j] = _positionMaxValue[j];}}//更新粒子群适应度this->refresh();}/**************************************************************** 函数名：findMax* 函数描述：采用粒子群算法搜索最优解* 输入参数：*  times：粒子群进化次数*  bestPartical：进化得到的最优个体*  disturbanceRate：粒子速度扰动概率，默认为0.2*  disturbanceVelocityCoe：速度扰动因子，表征扰动速度相对_maxSpeed大小*                          用于扰动粒子速度以提高局部搜索能力，默认为0.05* 输出参数：void***************************************************************/void findMax(int times,ZPSO_Partical& bestPartical,double disturbanceRate = 0.2,double disturbanceVelocityCoe = 0.05){this->randomlyInitial();for(int i=0;i<times;i++)this->update(disturbanceRate,disturbanceVelocityCoe);bestPartical.copy(_globalBestPartical);}/**************************************************************** 函数名：findMax* 函数描述：采用粒子群算法搜索最优解* 输入参数：*  times：粒子群进化次数*  bestParticalInEachLoop：每一次进化中的最优个体数组，*                          长度为times+1，由外部调用者提供内存空间*  disturbanceRate：粒子速度扰动概率，默认为0.2*  disturbanceVelocityCoe：速度扰动因子，表征扰动速度相对_maxSpeed大小*                          用于扰动粒子速度以提高局部搜索能力，默认为0.05* 输出参数：void***************************************************************/void findMax(int times,ZPSO_Partical *bestParticalInEachLoop,double disturbanceRate = 0.2,double disturbanceVelocityCoe = 0.05){this->randomlyInitial();for(int i=1;i<=times;i++){this->update(disturbanceRate,disturbanceVelocityCoe);bestParticalInEachLoop[i].copy(_globalBestPartical);}}
};#endif // _ZPSOALGORITHM_H

将以上两个文件直接放到一个文件夹中即可运行，为纯CPP文件。ZPSOAlgorithm.h则更类似于一个库文件，其定义的PSO算法类ZPSO_Algorithm可以解决一类最优化问题，而不限于某一个具体的待优化函数。

算法改进意见

PSO算法实现代码中，需要用户给定优化过程的粒子群进化次数，可以进一步改进PSO算法的收敛条件，改为当全局最优适应度在指定代数内都不发生改变时，判断算法收敛，退出循环。

[笔者注] PSO算法由J. Kennedy和R. C. Eberhart等在1995年提出，本文代码的设计与实现由作者个人设计与实现，仅用于PSO算法的记录、学习与分享。