粒子群算法简单实现

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]上的随机数

算法流程：

实例测试：

// 测试实例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，因此两组输出等价