粒子群算法简单实现

快来打我* 2022-03-07 01:22 216阅读 0赞

PSO算法相关定义

PSO中的每个粒子P包含两个向量: x, v

o位置向量x：粒子在解空间中的当前位置: x\[x1, x2, ..., xdim\]

o速度向量v：粒子在解空间中的飞行速度: v\[v1, v2, ..., vdim\]

pBest ：粒子自身的历史最优位置

gBest ：群体全局最优向量

粒子速度与位置更新：

P.v=omega\*P.v+c1\*r1\*(P.pBest-P.x)+c2\*r2\*(gBest-P.x);

P.x=P.x+P.v;

其中omega为惯性权重，c1,c2为学习因子，r1,r2为\[0, 1\]上的随机数

算法流程：

![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L0RYSDkyNA_size_16_color_FFFFFF_t_70][]![watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L0RYSDkyNA_size_16_color_FFFFFF_t_70 1][]

实例测试：

// 测试实例1
    // f=x1*x1+x2*x2+x3*x3, 求fmin
    // N=40, G=500
    // -10<=x<=-4
    #include<cstdio>
    #include<algorithm>
    #include<vector>
    #include<cstdlib>
    #include<ctime>
    #include<cmath>
    using namespace std;
    const int dim=3;
    const int N=40;
    const int G=500;
    const int low=-10, high=-4, wid=high-low;;
    const double INF=99999.9;
    double omega=0.729;
    double c1=1.49445, c2=c1;
    vector<double> gBest;
    // 定义粒子
    struct pp{
        vector<double> x;
        vector<double> v;
        vector<double> pBest;
        double value;
    }P[N];
    // 比较函数
    bool cmp(struct pp A, struct pp B) {
        return A.value<B.value;
    }
    int random() {
        int temp=rand()%(wid+1);
        return low+temp;
    }
    // 目标函数（适应度函数）
    double f(struct pp &P) {
        double vl=0;
        for (int i=0; i<dim; i++)
            vl+=P.x[i]*P.x[i];
        return vl;
    }
    // 更新操作
    void update(struct pp &P) {
        double r1=rand()/double(RAND_MAX);
        double r2=rand()/double(RAND_MAX);
        //P.v=omiga*P.v+c1*r1*(P.pBest-P.x)+c2*r2*(gBest-P.x);
        for (int i=0; i<dim; i++)
            P.v[i]=omega*P.v[i]+c1*r1*(P.pBest[i]-P.x[i])+c2*r2*(gBest[i]-P.x[i]);
        //P.x=P.x+P.v;
        for (int i=0; i<dim; i++) {
            P.x[i]=P.x[i]+P.v[i];
            if (P.x[i]>high)
                P.x[i]=high;
            else if (P.x[i]<low)
                P.x[i]=low;
        }
    }
    // 初始化
    void Init() {
        srand((unsigned)time(NULL));
        for (int i=0; i<N; i++) {
            for (int j=0; j<dim; j++) {
                P[i].v.push_back(random());
                P[i].x.push_back(random());
            }
            P[i].pBest=P[i].x;
            P[i].value=INF;
            P[i].value=min(f(P[i]), P[i].value);
        }
        sort(P, P+N, cmp);
        gBest=P[0].pBest;
    }
    // 运行
    void Run() {
        for (int i=0; i<G; i++) {
            for (int j=0; j<N; j++) {
                update(P[j]);
                double temp=f(P[j]);
                if (temp<P[j].value) {
                    P[j].pBest=P[j].x;
                    P[j].value=temp;
                }
                if (temp<P[0].value)
                    gBest=P[j].pBest;
            }
            sort(P, P+N, cmp);
        }
        for (int i=0; i<dim; i++)
            printf("x%d = %.2f\n", i+1, gBest[i]);
    }
    // 测试
    int main() {
        Init();
        Run();
        return 0;
    }
    /**********输出**********/
    x1 = -4.00
    x2 = -4.00
    x3 = -4.00
    
    //实例测试2
    // f=sinx1+cosx2-cos(2*x3)，求fmin
    // -10<=x<=-4
    // 目标函数（适应度函数）
    double f(struct pp &P) {
        double vl=0;
        vl=sin(P.x[0])+cos(P.x[1])-cos(2*P.x[2]);
        return vl;
    }
    /**********输出**********/
    x1 = -7.85
    x2 = -9.42
    x3 = -6.28
    或
    x1 = -7.85
    x2 = -9.42
    x3 = -9.42
    对于cos(2*x3)来说周期为pie=3.14，因此两组输出等价

[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L0RYSDkyNA_size_16_color_FFFFFF_t_70]: https://img-blog.csdnimg.cn/20190315222311473.png?x-oss-process=image/watermark,type_ZmFuZ3poZW5naGVpdGk,shadow_10,text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L0RYSDkyNA==,size_16,color_FFFFFF,t_70
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L0RYSDkyNA_size_16_color_FFFFFF_t_70 1]: /images/20220307/0692cc4d3e8c4930b95afea5f665bc9e.png