pso.h

基于粒子群优化算法

//PSO头文件PSO.hpp如下:
//Header __PSO_H
//Purpose Provide a class for particle swarm optimization
//Platform Visual C++
//Date 2008.05.05

//群微粒算法：本算法用群微粒算法求目标函数的最大值
//本算法的大致步骤和一些注意事项:
// (1) 派生自己的群微粒类，类中必须定义double GetFit(PARTICLE&)方法，用来计算每个微粒的适合度
// (2) 生成派生类实例，并在构造函数中指明微粒坐标维数和群体个数
// (3) 设置微粒坐标上界数组和下界数组，并用SetXup与SetXdown设置微粒坐标上下界
// (4) 用SetVmax方法设置微粒最大速度
// (5) 设置可选参数：C1，C2，W和通讯函数
// (6) 采用Run方法进行优化运算，优化后用GetBest方法获得最优个体适合度和坐标
// 基于粒子群的K均值聚类算法只需要将Run方法修改为K均值优化算法。

#ifndef __PSO_H
#define __PSO_H

//微粒类(class PARTICLE)定义
class PARTICLE
{
public:
double *X; //微粒的坐标数组
double *V; //微粒的速度数组
double *XBest; //微粒的最好位置数组
int Dim; //微粒的维数
double Fit; //微粒适合度
double FitBest; //微粒最好位置适合度

//构造函数
PARTICLE(); //空构造函数
PARTICLE(int n); //维数为参数的构造函数
//析构函数
~PARTICLE();

void SetDim(int d); //设置微粒的维数
};

//定义群粒子类
class PSO
{
	//保护成员
    protected:
       //变量
        PARTICLE *Particle; //微粒群数组
        int PNum; //微粒个数
        int GBestIndex; //最好微粒索引
        double W; //惯性权重
        double C1; //加速度系数1
        double C2; //加速度系数2
        double *Xup; //微粒坐标上界数组
        double *Xdown; //微粒坐标下界数组
        double *Vmax; //微粒最大速度数组
       //函数
        void Initialize(); //初始化群体
        void CalFit(); //计算全体适合度
        virtual void ParticleFly(); //微粒飞翔，产生新一代微粒
        //通讯函数，返回值为false时，系统停止优化
        bool (*Com)(double, //最优微粒适合度
        double*, //最优微粒坐标数组
        double**, //所有微粒坐标指针数组
        int); //当前最优微粒索引
    public:
        //构造函数
        PSO(); //空构造函数
        PSO(int dim, int num); //微粒个数,微粒维数为参数的构造函数
        //析构函数
        ~PSO();

        void SetXup(double*); //设置微粒坐标上界
        void SetXdown(double*); //设置微粒坐标下界
        void SetVmax(double*); //设置微粒最大速度，以数组为参数
        void SetVmax(double); //设置微粒最大速度，以上下界百分比为参数
        void SetW(double w){W=w;}; //设置权重
        void SetC1(double c){C1=c;};//设置加速度系数C1
        void SetC2(double c){C2=c;};//设置加速度系数C2
        void SetCom(void *p){Com=(bool(*)(double,double*,double**,int))p;};//设置通讯函数

        //计算特定微粒坐标所对应适合度，必须由派生的实际PSO类定义，以便计算适合度
        virtual double GetFit(PARTICLE&)=0;

        //运行类进行优化
        PARTICLE& Run(int max); //按最多次数限制运行PSO
        PARTICLE& Run(double fit); //按最佳适合度目标运行PSO
        double GetBest(double*); //获得最优微粒适合度和坐标
};

#endif
#include <stdlib.h>
#include <time.h>
#include <math.h>
//#include "PSO.H"

//微粒构造函数
PARTICLE::PARTICLE() //空构造函数
{
    X = 0; V = 0; XBest = 0; Dim = 0;
}
PARTICLE::PARTICLE(int n) //维数为参数的构造函数
{
    Dim = n;
    X = new double[Dim];
    V = new double[Dim];
    XBest = new double[Dim];
}

//微粒类析构函数
PARTICLE::~PARTICLE()
{
    if(Dim)
	{
        delete []X;
        delete []V;
        delete XBest;
	}
}

//设置微粒的维数
void PARTICLE::SetDim(int d)
{
   if(X) delete []X;
   if(V) delete []V;
   if(XBest) delete []XBest;
   Dim = d;
   X = new double[Dim];
   V = new double[Dim];
   XBest = new double[Dim];
}

//PSO构造函数
PSO::PSO()
{
    Particle = 0;
    PNum = 0;
    GBestIndex = 0;
    Xup = 0;
    Xdown = 0;
    W = 1;
    C1 = 2;
    C2 = 2;
    Com = 0;
}
//微粒维数、微粒个数为参数的构造函数
PSO::PSO(int dim, int num)
{
   Particle = new PARTICLE[num];
   for(int i=0; i<num; i++) 
	   Particle[i].SetDim(dim);
   PNum = num;
   GBestIndex = 0;
   Xup = new double[dim];
   Xdown = new double[dim];
   Vmax = new double[dim];
   W = 1;
   C1 = 2;
   C2 = 2;
   Com = 0;
}

//析构函数
PSO::~PSO()
{
    if(Particle) delete []Particle;
    if(Xup) delete []Xup;
    if(Xdown) delete []Xdown;
    if(Vmax) delete []Vmax;
}
//设置坐标上界
void PSO::SetXup(double *up)
{
    if(!Particle) return;
    for(int i=0; i<Particle[0].Dim; i++)
    Xup[i] = up[i];
}
//设置坐标下界
void PSO::SetXdown(double *d)
{
    if(!Particle) return;
    for(int i=0; i<Particle[0].Dim;i++)
    Xdown[i] = d[i];
}
//设置最大速度
void PSO::SetVmax(double *max)
{
    if(!Particle) return;
    for(int i=0; i<Particle[0].Dim;i++)
    Vmax[i] = max[i];
}
void PSO::SetVmax(double p)
{
    if(!Particle) return;
    for(int i=0; i<Particle[0].Dim; i++)
		Vmax[i] = (Xup[i]-Xdown[i])*p;
}

//初始化群体
void PSO::Initialize()
{
    if(!Particle) return;
    static int kk=(unsigned)time(NULL);  //
    srand((unsigned)time(NULL)+kk++);    //
    GBestIndex = 0;
    for(int i=0; i<PNum; i++)
	{
        for(int j=0; j<Particle[i].Dim; j++)
		{
            Particle[i].X[j] = rand()/(double)RAND_MAX*(Xup[j]-Xdown[j])+Xdown[j];//初始化坐标
            Particle[i].XBest[j] = Particle[i].X[j]; //初始化最好位置微粒的坐标
            Particle[i].V[j] = rand()/(double)RAND_MAX*Vmax[j]-Vmax[j]/2;//初始化速度
		}
    Particle[i].Fit = GetFit(Particle[i]); //计算该微粒适合度
    Particle[i].FitBest = Particle[i].Fit; //设最优适合度初值
    if(Particle[i].Fit>Particle[GBestIndex].Fit) GBestIndex = i;//查找群体最优微粒
	}
}

//计算群体各个微粒适合度
void PSO::CalFit()
{
    if(!Particle) return;
    for(int i=0; i<PNum; i++)
		Particle[i].Fit = GetFit(Particle[i]);
}

//微粒飞翔,产生新一代微粒
void PSO::ParticleFly()
{
	static double FitBak[100];
	if(!Particle) return;
	static int tt=(unsigned)time(NULL);
	srand((unsigned)time(NULL)+tt++);
	//整个群体飞向新的位置
	for(int i=0; i<PNum; i++)
	{
		for(int j=0; j<Particle[i].Dim; j++)
			Particle[i].V[j] = W*Particle[i].V[j]+// (FitBak[i]-Particle[i].Fit)+//修改速度
				rand()/(double)RAND_MAX*C1*(Particle[i].XBest[j]-Particle[i].X[j])+
				rand()/(double)RAND_MAX*C2*(Particle[GBestIndex].XBest[j]-Particle[i].X[j]);
		for(j=0; j<Particle[i].Dim; j++) //检查速度最大值
		{
			if(Particle[i].V[j]>Vmax[j]) Particle[i].V[j] = Vmax[j];
			//若微粒速度大于最大速度则等于最大速度
			if(Particle[i].V[j]<-Vmax[j]) Particle[i].V[j] = -Vmax[j];
			//若微粒速度小于最小速度则等于最小速度
		}
		for(j=0; j<Particle[i].Dim; j++)
		{
			Particle[i].X[j] += Particle[i].V[j]; //修改坐标
			if(Particle[i].X[j]>Xup[j]) Particle[i].X[j]=Xup[j];
			//若速度大于最大速度则等于最大速度,起到保护作用
			if(Particle[i].X[j]<Xdown[j]) Particle[i].X[j]=Xdown[j];
		}
	}
	//计算各微粒适合度
	CalFit();  //calculate fit value
	for(i=0; i<PNum; i++) FitBak[i] = Particle[i].Fit;
	//设置新的个体最好位置
	for(i=0; i<PNum; i++)
	   if(Particle[i].Fit>=Particle[i].FitBest)
	   {
		   Particle[i].FitBest = Particle[i].Fit;
		   for(int j=0; j<Particle[i].Dim; j++)
		      Particle[i].XBest[j] = Particle[i].X[j];
	   }
	//设置新的最优个体
	GBestIndex = 0;
	for(i=0; i<PNum; i++)
	   if(Particle[i].FitBest>=Particle[GBestIndex].FitBest && i!=GBestIndex) GBestIndex = i;
}
///*
//运行群粒算法
PARTICLE& PSO::Run(int n)
{
   Initialize();
   double *opt_p = new double[Particle[0].Dim]; //通讯用数组，最优点坐标
   double **opt_a = new double*[PNum]; //通讯用数组，所有点坐标
   for(int i=0; i<n; i++)
   {
      ParticleFly();
      if(Com) //通讯函数存在，完成通讯
	  {
		  //拷贝最优点坐标
          for(int k=0; k<Particle[0].Dim; k++) opt_p[k]=Particle[GBestIndex].XBest[k];
          for(k=0; k<PNum; k++) opt_a[k] = Particle[k].X; //指向所有点坐标
          if(!Com(Particle[GBestIndex].FitBest,opt_p,opt_a,GBestIndex)) break;
	  }
   }
   delete []opt_p;
   delete []opt_a;
   return Particle[GBestIndex];
}
PARTICLE& PSO::Run(double fit)
{
   double *opt_p = new double[Particle[0].Dim]; //通讯用数组，最优点坐标
   double **opt_a = new double*[PNum]; //通讯用数组，所有点坐标
   Initialize();
   do
   {
       ParticleFly();
       if(Com) //通讯函数存在，完成通讯
	   {
           for(int k=0; k<Particle[0].Dim; k++) 
			   opt_p[k]=Particle[GBestIndex].XBest[k];//拷贝最优点坐标
           for(k=0; k<PNum; k++) 
			   opt_a[k] = Particle[k].X; //指向所有点坐标
           if(!Com(Particle[GBestIndex].FitBest,opt_p,opt_a,GBestIndex)) break;
	   }
   }while(Particle[GBestIndex].FitBest<fit);
   delete []opt_p;
   delete []opt_a;
   return Particle[GBestIndex];
}
//*/
//返回最佳个体
double PSO::GetBest(double *r)
{
    for(int i=0; i<Particle[GBestIndex].Dim; i++)
        r[i] = Particle[GBestIndex].XBest[i];
    return Particle[GBestIndex].FitBest;
}