粒子群算法(Particle swarm optimization, PSO)原理和C++实现

RobotsRuning

已于 2023-06-25 21:43:21 修改

阅读量1.7k

点赞数 3

文章标签：数学建模 c++ 算法

于 2023-06-25 10:08:47 首次发布

本文链接：https://blog.csdn.net/u014374826/article/details/131367195

版权

文章介绍了粒子群优化(PSO)算法的基本思想和更新规则，这是一种模拟生物群体行为的优化算法，适用于解决函数优化等问题。通过跟踪个体和全局最优解，粒子群中的每个成员更新其速度和位置，最终找到全局最优解。文中还给出了C++代码示例，展示如何实现PSO算法来求解特定优化问题。

摘要由CSDN通过智能技术生成

1. 算法简介

粒子群优化（Particle Swarm Optimization，PSO）是一种模拟生物群体行为的优化算法，其灵感来源于鸟群捕食的行为。每个“粒子”（或“鸟”）在搜索空间中飞翔，寻找优化目标的最优解（个体认知（Cognition-only model)）。粒子之间可以分享彼此的信息（社会行为（Social-only model)），从而找到全局最优解。它与其他优化算法的不同之处在于，它只需要目标函数，不依赖于目标的梯度或任何微分形式。它也有很少的超参数。

PSO的优势：在于简单容易实现并且没有许多参数的调节。目前已被广泛应用于函数优化、神经网络训练、模糊系统控制以及其他遗传算法的应用领域。

2. 算法介绍：

2.1 基本思想

粒子群算法通过设计一种无质量的粒子来模拟鸟群中的鸟，粒子仅具有两个属性：速度和位置，速度代表移动的快慢，位置代表移动的方向。每个粒子在搜索空间中单独的搜寻最优解，并将其记为当前个体极值，并将个体极值与整个粒子群里的其他粒子共享，找到最优的那个个体极值作为整个粒子群的当前全局最优解，粒子群中的所有粒子根据自己找到的当前个体极值和整个粒子群共享的当前全局最优解来调整自己的速度和位置。

2.2 更新规则

PSO初始化为一群随机粒子(随机解)。然后通过迭代找到最优解。在每一次的迭代中，粒子通过跟踪两个“极值”(pbest，gbest)来更新自己。在找到这两个最优值后，粒子通过下面的公式来更新自己。

公式 (1):

公式 (2):

在公式(1))中， $i$ =1,2,..., $N$ ， $N$ 是此群中粒子的总数； $v_{i}$ ：粒子速度，最大值为 $V_{max}$ ，如果 $v_{i}$ 大于 $V_{max}$ ，则 $v_{i}$ = $V_{max}$ ； $rand()$ ：介于(0,1)之间的随机数； $x_{i}$ ：粒子的当前位置； $c_{1}$ 和 $c_{2}$ ：学习因子，通常均为2。

公式(1)的第一部分称为【记忆项】，表示上次速度大小和方向的影响；公式(1)的第二部分称为【自身认知项】，是从当前点指向粒子自身最好点的一个矢量，表示粒子的动作来源于自己经验的部分；公式(1)的第三部分称为【群体认知项】，是一个从当前点指向种群最好点的矢量，反映了粒子间的协同合作和知识共享。粒子就是通过自己的经验和同伴中最好的经验来决定下一步的运动。以上面两个公式为基础，形成了PSO的标准形式。

公式(2)和公式(3)被视为标准PSO算法。

3. PSO算法的流程和伪代码

4. PSO算法举例

4.1 待优化问题

4.2 C++代码实现

#include <iostream>
#include <vector>
#include <cmath>
#include <map>
#include <algorithm>
#include <random>
#include <ctime>
#include <Eigen/Dense>
using namespace Eigen;
using namespace std;

const int dim = 1;//维数
const int p_num = 10;//粒子数量
const int iters = 100;//迭代次数
const int inf = 999999;//极大值
const double pi = 3.1415;
//定义粒子的位置和速度的范围
const double v_max = 4;
const double v_min = -2;
const double pos_max = 2;
const double pos_min = -1;
//定义位置向量和速度向量
vector<double> pos;
vector<double> spd;
//定义粒子的历史最优位置和全局最优位置
vector<double> p_best;
double g_best;
//使用eigen库定义函数值矩阵和位置矩阵
Matrix<double, iters, p_num> f_test;
Matrix<double, iters, p_num> pos_mat;

//定义适应度函数
double fun_test(double x)
{
    double res = x * x + 1;
    return res;
}

//初始化粒子群的位置和速度
void init()
{
    //矩阵中所有元素初始化为极大值
    f_test.fill(inf);
    pos_mat.fill(inf);
    //生成范围随机数
    static std::mt19937 rng;
    static std::uniform_real_distribution<double> distribution1(-1, 2);
    static std::uniform_real_distribution<double> distribution2(-2, 4);
    for (int i = 0; i < p_num; ++i)
    {
        pos.push_back(distribution1(rng));
        spd.push_back(distribution2(rng));
    }
    vector<double> vec;
    for (int i = 0; i < p_num; ++i)
    {
        auto temp = fun_test(pos[i]);//计算函数值
        //初始化函数值矩阵和位置矩阵
        f_test(0, i) = temp;
        pos_mat(0, i) = pos[i];
        p_best.push_back(pos[i]);//初始化粒子的历史最优位置
    }
    std::ptrdiff_t minRow, minCol;
    f_test.row(0).minCoeff(&minRow, &minCol);//返回函数值矩阵第一行中极小值对应的位置
    g_best = pos_mat(minRow, minCol);//初始化全局最优位置
}

void PSO()
{
    static std::mt19937 rng;
    static std::uniform_real_distribution<double> distribution(0, 1);
    for (int step = 1; step < iters; ++step)
    {
        for (int i = 0; i < p_num; ++i)
        {
            //更新速度向量和位置向量
            spd[i] = 0.5 * spd[i] + 2 * distribution(rng) * (p_best[i] - pos[i]) +
                2 * distribution(rng) * (g_best - pos[i]);
            pos[i] = pos[i] + spd[i];
            //如果越界则取边界值
            if (spd[i] < -2 || spd[i] > 4)
                spd[i] = 4;
            if (pos[i] < -1 || pos[i] > 2)
                pos[i] = -1;
            //更新位置矩阵
            pos_mat(step, i) = pos[i];
        }
        //更新函数值矩阵
        for (int i = 0; i < p_num; ++i)
        {
            auto temp = fun_test(pos[i]);
            f_test(step, i) = temp;
        }
        for (int i = 0; i < p_num; ++i)
        {
            MatrixXd temp_test;
            temp_test = f_test.col(i);//取函数值矩阵的每一列
            std::ptrdiff_t minRow, minCol;
            temp_test.minCoeff(&minRow, &minCol);//获取每一列的极小值对应的位置
            p_best[i] = pos_mat(minRow, i);//获取每一列的极小值，即每个粒子的历史最优位置
        }
        g_best = *min_element(p_best.begin(), p_best.end());//获取全局最优位置
    }
    cout << fun_test(g_best);
}

int main()
{
    init();
    PSO();
    system("pause");
    return 0;
}

#include <iostream>
#include<time.h>
#include <random>
 
using namespace std;
 
//使用c++实现粒子群算法
//需要寻优的非线性函数为：
//f(x,y) = sin(sqrt(x^2+y^2))/(sqrt(x^2+y^2)) + exp((cos(2*PI*x)+cos(2*PI*y))/2) - 2.71289
//该函数有很多局部极大值点，而极限位置为(0,0),在(0,0)附近取得极大值
 
// Author: yuzewei
 
const int N = 20;             //粒子群数量
const int dim = 2;            //粒子群维度
const int gen = 300;//100;    //迭代次数
const double PI = 3.1415926;
const double w = 0.9;//0.6;   //惯性权重
const double c1 = 2;// 1.49445;     //自学习，学习因子
const double c2 = 2; // 1.49445;    //群体最优学习，学习因子
const double pmin = -2.0;     //粒子的移动区间范围
const double pmax = 2.0;
const double vmin = -0.8;// -0.5;   //粒子的移动速度范围
const double vmax = 0.8;// 0.5;
 
double P[N][dim];       //粒子位置
double V[N][dim];       //粒子速度
double fitness[N];      //粒子适应度
double pbest[N][dim];   //粒子最优位置
 
double gbest[dim];      //粒子群最优位置
int gbestGen = -1;      //记录达到最优位置时的迭代次数
 
double func(double *p)
{
    double x = *p;
    double y = *(p + 1);
    double result = sin(sqrt(x * x + y * y)) / (sqrt(x * x + y * y)) 
                    + exp((cos(2 * PI * x) + cos(2 * PI * y)) / 2) - 2.71289;
    return result;
}
 
int getGbestIndex(double *fit)
{
    int index = -1;
    double max = *fit;
    for (int i = 0; i < N; i++)
    {
        if (*(fit + i) > max)
        {
            max = *(fit + i);
            index = i;
        }
    }
    return index;
}
 
void init()
{
    int index = -1; //记录位置最佳的粒子编号
 
    default_random_engine e;
    uniform_real_distribution<double> rand1(pmin, pmax);
    uniform_real_distribution<double> rand2(vmin, vmax);
    e.seed(time(0));
 
    for (int i = 0; i < N; i++)
    {
        for (int j = 0; j < dim; j++)
        {
            P[i][j] = rand1(e);
            V[i][j] = rand2(e);
            pbest[i][j] = P[i][j];
        }
        fitness[i] = func(P[i]);
    }
 
    index = getGbestIndex(fitness);
    gbest[0] = pbest[index][0];
    gbest[1] = pbest[index][1];
    gbestGen = 0;
}
 
 
void printParticles()
{
    /*for (int i = 0; i < N; i++)
    {
        cout << "第" << i << "个粒子" << endl;
        cout << "  位置：(" << P[i][0] << "," << P[i][1] << ")" << endl;
        cout << "  速度：(" << V[i][0] << "," << V[i][1] << ")" << endl;
        cout << "  适应度：" << fitness[i] << endl;
        cout << "  个体最优：(" << pbest[i][0] << "," << pbest[i][1] << ")" << endl;
    }*/
    //cout << endl;
    cout<<"群体最优：(" << gbest[0] << "," << gbest[1] << ")" << endl;
    cout << "适应度为：" << func(gbest) << endl;
}
 
void PSOiterator()
{
    cout << "开始迭代！" << endl;
    cout << "第1代粒子群：" << endl;
    printParticles();
 
    for (int g = 0; g < gen - 1; g++)
    {
        cout << "第" << g + 2 << "代粒子群：" << endl;
 
        int index = -1; //记录最优粒子个体编号
 
        default_random_engine e;
        uniform_real_distribution<double> rand(0, 1.0);
        e.seed(time(0));
 
        for (int i = 0; i < N; i++)
        {
            for (int j = 0; j < dim; j++)
            {
                //更新速度
                V[i][j] = w * V[i][j] + rand(e) * (pbest[i][j] - P[i][j])
                    + rand(e) * (gbest[j] - P[i][j]);
                if (V[i][j] < vmin)
                    V[i][j] = vmin;
                if (V[i][j] > vmax)
                    V[i][j] = vmax;
 
                //更新位置
                P[i][j] = P[i][j] + V[i][j];
                if (P[i][j] < pmin)
                    P[i][j] = pmin;
                if (P[i][j] > pmax)
                    P[i][j] = pmax;
            }
 
            fitness[i] = func(P[i]); //更新适应度
 
            //更新个体最优
            if (fitness[i] > func(pbest[i]))
            {
                pbest[i][0] = P[i][0];
                pbest[i][1] = P[i][1];
            }
 
        }
 
        //更新整体最优
        index = getGbestIndex(fitness);
        if (func(P[index]) > func(gbest))
        {
            gbest[0] = P[index][0];
            gbest[1] = P[index][1];
            gbestGen = g + 1;
        }
 
        printParticles();
    }
}
 
void printResult()
{
    cout << "迭代结束！" << endl;
    cout << "总共进行了" << gen << "次迭代，";
    cout<<"第"<<gbestGen<<"代得到最优位置:(" << gbest[0] << "," << gbest[1] << ")" << endl;
    cout << "此位置对应的函数值为：" << func(gbest) << endl;
}
 
int main()
{
    init();
    PSOiterator();
    printResult();
    return 0;
}

算法优点：

粒子群算法是一门新兴算法，此算法与遗传算法有很多相似之处，其收敛于全局最优解的概率很大（全局最优概率大）。

相较于传统算法计算速度非常快，全局搜索能力很强（收敛速度快）；

PSO对于种群大小不十分敏感，所以初始种群设为500-1000，速度影响也不大；

粒子群算法适用于连续函数极值问题，对于非线性、多峰问题均有较强的全局搜索能力；

设置参数少。