ps.cpp

标准粒子群优化算法源码

/*********************************************************
粒子群种群类
*********************************************************/
#include "stdafx.h"
#include "PS.h"
#include "time.h"

#define pi		3.1415926535897932

/*********************************************************
粒子评价函数：
由于适应度的评价与具体问题相关
因而需要一个外部函数评价粒子的适应度
*********************************************************/
extern double Evaluate(double *pw,int &err);
FILE *m_file;

//Rosenbrock function
double CPS::EvalFunc(double *x)
{
	int i;
	double value=0.;
	for(i=0;i<m_ParamNumber-1;i++)
	{
		value+=100*(x[i]*x[i]-x[i+1])*(x[i]*x[i]-x[i+1]);
		value+=(x[i]-1.)*(x[i]-1.);
	}
	return value;
}

/*
//Rastrigin function
double CPS::EvalFunc(double *x)
{
	int i;
	double value=0.;
	for(i=0;i<m_ParamNumber-1;i++)
	{
		value+=(x[i]*x[i]-10.*cos(2*pi*x[i])+10.);
	}
	return value;
}
*/
CPS::CPS()
{
	CPS(20,345,-.5,.5);
}

/*********************************************************
粒子群类构造函数：
需要传入：
pop_size－种群规模
param_number－要优化的参数个数
v_l－参数左边界
v_h－参数右边界
*********************************************************/
CPS::CPS(int pop_size,int param_number,double v_l,double v_h)
{
	this->m_PopSize=pop_size;
	this->m_ParamNumber=30;
	this->m_vl=-10.;
	this->m_vh=10.;
	this->m_vmax=0.5;
	this->m_sparrow=NULL;
	this->m_gbest=NULL;
	this->m_iter=0;
	m_file=fopen("proc_res.txt","w");
	if(m_file)
		fprintf(m_file,"Standard PSO Algorithm\n");
	InitPop();
}

CPS::~CPS()
{
	if(m_file)
		fclose(m_file);
	DestroyPop();
}

/*********************************************************
粒子群类初始化函数：
完成：
0 为个体、个体速度、个体最优解和全局最优解缓存分配内存
1 初始化个体值
2 初始化个体速度
*********************************************************/
void CPS::InitPop()
{
	int r;
	int i,j;
	double v;
	DestroyPop();
	srand( (unsigned)time( NULL ) );
	m_sparrow=new SparrowParticle[m_PopSize];
	this->m_gbest=new double[m_ParamNumber];
	this->m_fgbest=(double)0xffffffff;

	for(i=0;i<m_PopSize;i++)
	{
		m_sparrow[i].x=new double[m_ParamNumber];
		m_sparrow[i].v=new double[m_ParamNumber];
		m_sparrow[i].best_x=new double[m_ParamNumber];
		m_sparrow[i].best_f=(double)0xffffffff;
		for(j=0;j<m_ParamNumber;j++)
		{
			r=rand();
			m_sparrow[i].x[j]=(r*(m_vh-m_vl))/RAND_MAX+m_vl;
			r=rand();
			v=(r*m_vmax)/RAND_MAX;
			m_sparrow[i].v[j]=(r%2)?-v:v;
			m_sparrow[i].best_x[j]=m_sparrow[i].x[j];
		}
	}
}

/*********************************************************
粒子群类销毁函数：
完成：
0 销毁各指针，释放内存
*********************************************************/
void CPS::DestroyPop()
{
	int i;
	if(m_sparrow)
	{
		for(i=0;i<m_PopSize;i++)
		{
			delete []m_sparrow[i].x;
			delete []m_sparrow[i].v;
			delete []m_sparrow[i].best_x;
		}
		delete []m_sparrow;
		m_sparrow=NULL;
	}
	if(this->m_gbest)delete []m_gbest;
	this->m_gbest=NULL;
}

/*********************************************************
粒子群迭代进化函数：
w－惯性系数，随迭代次数增加呈线形递减
*********************************************************/
void CPS::Propagation(double w)
{
	int i,j;
	double f;

	for(i=0;i<m_PopSize;i++)
	{
		f=CountFit(m_sparrow[i].x);
		if(f<m_sparrow[i].best_f)				//若该解优于个体历史最优解，更新个体最优解
		{
			m_sparrow[i].best_f=f;
			memcpy(m_sparrow[i].best_x,m_sparrow[i].x,sizeof(double)*m_ParamNumber);
		}
		if(f<m_fgbest)							//若该解优于全局最优解，更新全局最优解
		{
			m_fgbest=f;
			memcpy(m_gbest,m_sparrow[i].x,sizeof(double)*m_ParamNumber);
		}
	}
	if(m_file)
		fprintf(m_file,"%04d:  %f\n",m_iter,m_fgbest);
	m_iter++;
	for(i=0;i<m_PopSize;i++)
	{
		for(j=0;j<m_ParamNumber;j++)			//更新粒子速度
		{
			m_sparrow[i].v[j]=m_sparrow[i].v[j]*w+2*rand()*(m_sparrow[i].best_x[j]-
				m_sparrow[i].x[j])/RAND_MAX+2*rand()*(m_gbest[j]-m_sparrow[i].x[j])/RAND_MAX;
			if(m_sparrow[i].v[j]>m_vmax)
				m_sparrow[i].v[j]=m_vmax;
			if(m_sparrow[i].v[j]<-m_vmax)
				m_sparrow[i].v[j]=-m_vmax;
			m_sparrow[i].x[j]+=m_sparrow[i].v[j];
		}
	}
}

void CPS::Iterate(int max_iter)
{
	int i,j;
	double w;
	double avg_f;
	int max_run=200;
	for(j=0,avg_f=0.;j<max_run;j++)
	{
		InitPop();
		for(i=0;i<max_iter;i++)
		{
//			cout<<i+1<<":   ";
			w=.5*(max_iter-i)/max_iter+.4;
			Propagation(w);
//			cout<<m_fgbest<<endl;
		}
		avg_f+=m_fgbest;
		printf("%02d:     %f\n",j,m_fgbest);
	}
	printf("Result:    %f\n",avg_f/max_run);
}

/*********************************************************
粒子适应度计算函数，调用外部评价函数
*********************************************************/
double CPS::CountFit(double *x)
{
	return EvalFunc(x);
}

int CPS::CountErr(double *x)
{
	int err;
	Evaluate(x,err);
	return err;
}