pso.cpp

PSO算法的源代码

/*****************************************

文件名:  PSO.cpp
粒子群算法求函数优化。
console application, 行数: 300 lines。
modifyed by gaochong。 finished at 2008-04-07
PSO算法模板
默认是求最小值
具体情况要进行修改的是： *需要更改* ，可以搜索并进行修改 

******************************************/

#include <cstdio>   //include printf() 
#include <cstdlib>  //包含rand()和srand()
#include <cmath>    //include cos(),sin() 
#include <ctime>    //include time()

#include <limits>	//使用double 的最大值
#include <iostream> //include cout cin 

using std::numeric_limits;		//使用double 的最大值,numeric_limits<double>::max()
using std::cout;
using std::endl;

#define GANVARS 1	//参数的维数             *****需要更改******
#define PSO_popsize 30                        
#define PSO_maxgen 1000
#define Vmax 0.5 //速度极值，如果取值过大，不利于个体的优化,一般是upbound-lowbound

//生成[low,uper]的随机double值 
#define rnd(low,uper) ((rand()/(double)RAND_MAX)*((uper)-(low))+(low))
#define RUNTIMES 10   //程序运行次数 


double paramater_w; //历史系数， 由程序随机生成 
double paramater_c1=1.4; //认知系数，一般取2.0 
double paramater_c2=1.4; //社会系数, 一般取2.0 
double paramate_cy = 1.0e-10;

double lowBound[GANVARS],upperBound[GANVARS]; //种群中个体的范围，   *****需要更改****** 

//全局极值的坐标,注意由ParticleSwarm::getGBest()更新 
//粒子更新速度需要使用，用于全局通信用。 
double gbest_pos_global[GANVARS]; 


/************************************************************************/
/* 单个粒子类                                                           */
/************************************************************************/
class Particle {
public:
	double pos[GANVARS]; 		//粒子位置
	double v[GANVARS]; 			//粒子速度
	double pbest; 
	double pbest_pos[GANVARS];//个体最优解的坐标，对应每种证券组合比例的值
	double fitness; //当前算出的一个适应值
public:
	double calcFitness(double pos[]); 
	void updatePosition(); 
	void updatePBest();	 
};

//计算适应值的函数,评估个体，如果有约束，需要加罚函数呀！ ****需要更改*****
double Particle::calcFitness(double *p)
{
	
    int i,j;
    double serr=1,temp=0;
	for(i=1;i<=GANVARS;i++)
	{
		for(j=1;j<=5;j++)
		{
			temp=temp+j*cos(  (j+1) * p[i-1] + j  );
		}
	serr=serr*temp;
	}
return serr;
}

//更新位移和速度
void Particle::updatePosition()
{
	double r1,r2;
	int i,j,k;
	
	for(i=0;i<GANVARS;i++)
    {	
		paramater_w = rnd(0,1);
		
        //更新速度，利用粒子的pbest_pos，和全局的  gbest_pos_global[]
 	    v[i] = paramater_w  * v[i] +
               paramater_c1 * rnd(0,1) * (pbest_pos[i]        -  pos[i]) +
               paramater_c2 * rnd(0,1) * (gbest_pos_global[i] -  pos[i]); 

		//判断超出最大速度和最小速度。
		if (v[i]<-Vmax)
			v[i] = -Vmax; 
		if (v[i]>Vmax)
			v[i]=Vmax;

		//极值扰动吗?
		r1 = rnd(0,1);
		if (r1 < paramate_cy)
		{
			r2 = rnd(0,1);
			v[i]=r2*v[i];
		}

		//更新位移
		pos[i]+=v[i]; 

		//出界则拉回 
		if(pos[i]<lowBound[i])
			pos[i]=lowBound[i];
		if(pos[i]>upperBound[i])
			pos[i]=upperBound[i];
	}
}

//更新个体极值	
void Particle::updatePBest(){
	
	if(this->fitness<pbest)
	{
		pbest=this->fitness;
		for(int i=0;i<GANVARS;i++)
		{
			pbest_pos[i]=pos[i];	
		}
	}
}

/************************************************************************/
/*        粒子群类                                                      */
/************************************************************************/
class ParticleSwarm{
public:
	double gbest; //全局极值的适应值
	double gbest_pos[GANVARS]; //全局极值的坐标
    Particle PSO_pop[PSO_popsize];//单个粒子定义为粒子群类的属性
public: 
	void init();  //初始化种群
	void getGBest(); //获取全局极值
	void search(double *Array); //迭代,col参数是[0,RUNTIME-1]，指示是第几次运行 	  
};

void ParticleSwarm::init(){
	
	gbest=numeric_limits<double>::max();
	srand((unsigned)time(NULL));

	//初始化Bound ****需要更改*****
	for(int t=0;t<GANVARS;t++)
	{
		lowBound[t]=-10;
		upperBound[t]=10;
	}

	//初始化种群,包括位置，速度，fitness和pbest 
	for(int i=0;i<PSO_popsize;i++)
	{   
        //x[] 保存位置。y[]保存速度 
		double x[GANVARS];
		double y[GANVARS];
	
		for(int j=0;j<GANVARS;j++)
		{		
			x[j] = rnd(lowBound[j],upperBound[j]);
            y[j] = rnd(-Vmax,Vmax);//[-Vmax,Vmax]之间的随机数 
		}
		
        //初始化位置和速度	
		for(int j=0;j<GANVARS;j++){
			PSO_pop[i].pos[j]=x[j];
			PSO_pop[i].v[j]=y[j];
		}
		
		//计算fitness 
		PSO_pop[i].fitness= PSO_pop[i].calcFitness(PSO_pop[i].pos);
		
		//初始化，将当前fitness赋给pbest
		PSO_pop[i].pbest=PSO_pop[i].fitness;  
		for(int m=0;m<GANVARS;m++)
        {
			PSO_pop[i].pbest_pos[m]=PSO_pop[i].pos[m];
		}
	}	
	
	getGBest();
	
    cout<<"gen: 0"<<endl;  //初始化算作第零代。 
	cout<<"The lowest value is" <<gbest<<endl;
	cout<<"The parameters are:"<<endl;
	for(int m=0;m<GANVARS;m++)
	{
		cout<<gbest_pos[m]<<" ";
	}
	cout<<endl;
}

void ParticleSwarm::getGBest()
{
	for(int i=0;i<PSO_popsize;i++)
	{
		if(PSO_pop[i].fitness<gbest)
		{
			gbest=PSO_pop[i].fitness;
			for(int j=0;j<GANVARS;j++)
			{
				gbest_pos[j]=PSO_pop[i].pos[j];	
				//也要更新全局的 gbest_pos_global,用于全体粒子通讯用。 
				gbest_pos_global[j]=PSO_pop[i].pos[j];		  
			}
			
		}
	}
}

//迭代
void ParticleSwarm::search(double *Array)
{
	int gen=0;
	
	//代数范围是:[0, PSO_maxgen-1],总共PSO_maxgen代。 
	while(gen<PSO_maxgen)
    	{
        Array[gen] = gbest;
        
        //每个粒子进行运动，求值，更行pbest 
		for(int k=0;k<PSO_popsize;k++)
		{
			PSO_pop[k].updatePosition();
			PSO_pop[k].fitness=PSO_pop[k].calcFitness(PSO_pop[k].pos);
			PSO_pop[k].updatePBest();
		}
		
		//更新gbest.注意包含更新用于通信的全局变量 
        	getGBest();

		gen++;		
		cout<<"gen:"<<gen<<endl;
		cout<<"The lowest value is" <<gbest<<endl;
		cout<<"The parameters are:"<<endl;
		for(int i=0;i<GANVARS;i++)
			cout<<gbest_pos[i]<<" ";
		cout<<endl;
	
	}
}


int main(int argc, char* argv[])
{
	ParticleSwarm p;
	int i,j;

	FILE *PSOData;
	double  *resultStore;
	double  outputdata = 0.0;
	
	//在堆上分配resultStore，保存每一代的gbest。是数组。每运行一次，保存一次。 
	resultStore = new double[PSO_maxgen];
	for (i = 0; i < PSO_maxgen; i++)
	{
			resultStore[i] = 0.0;
	}
	
	PSOData = fopen("PSO.dat","a+");//添加到文件末尾。 
	if (PSOData == NULL)
	{
		printf("The file is not exit!\n");
		exit(1);
	}
	else
	{		
		for (int i = 0; i < RUNTIMES; i++)
		{
			p.init();
			p.search(resultStore);	
			
			fprintf(PSOData,"Run number:  %d, Best values:  %f  \n",i,p.gbest);
			for (int j=0;j<GANVARS;j++)
			{
				fprintf(PSOData,"%f ",p.gbest_pos[j]);
			}
			fprintf(PSOData,"\n");
			fprintf(PSOData,"-----------------------------------------------------------------------------------------\n");
			
			//保存每一代的gbest,可以屏蔽掉下面两行。使文件简洁易读。 
			for (int row = 0; row <PSO_maxgen; row++)				
				fprintf(PSOData,"%d   %f\n",row,resultStore[row]);					
		}
	}
	
	//扫尾工作 
   	delete  []resultStore;
	fclose(PSOData);
	system("pause");
	return 0;
}