simplepso.cpp

粒子群算法的一个简单的C实现

//PSO
//completed at 13:47 2006.3.21

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>

#define MYTYPE double
#define KMAX 50

// declare functions
/*------------------------------------------------------------------*/
MYTYPE fun(MYTYPE);  //待求最优解问题
MYTYPE fun1(MYTYPE);
MYTYPE getmnrand(MYTYPE ,MYTYPE );  //随机数产生函数
//MYTYPE gaussrand(void);
MYTYPE pso(int , MYTYPE , MYTYPE ,MYTYPE (*function)(MYTYPE ),int (*cmpbetter)(MYTYPE, MYTYPE));  //算法主函数
int itok(int);  //退出计算条件判断函数
int comparebt(MYTYPE,MYTYPE);  //比较函数
MYTYPE abbs(MYTYPE i);  //绝对值函数
double get_rand(double low_end, double high_end);
void prin(int i,MYTYPE * p,char *s); //打印函数
/*------------------------------------------------------------------*/

int main(void)
{
	MYTYPE it;
	srand( (unsigned)time( NULL ) );
	it=pso(20, 0.0, 100.0, fun, comparebt );
	pso(20, 0.0, 100.0, fun1, comparebt );
	return 0;
}

// My PSO Algorithm
/*------------------------------------------------------------------*/
MYTYPE pso(int pinit, MYTYPE a, MYTYPE b, MYTYPE (*function)(MYTYPE ),int (*cmpbetter)(MYTYPE, MYTYPE))
{
	if(b<=a)
	{exit(0);return 0.0;}

	int i,FLAG=0;
	long num=0;  //叠代计算次数
	
	MYTYPE *particle=(MYTYPE *)calloc(sizeof(MYTYPE), (unsigned)pinit);
	for(i=0;i<pinit;i++)
		*(particle+i)=getmnrand(a,b);  //粒子位置初始化

	prin(pinit,particle,"particle");   //打印初始化结果
	
	MYTYPE *velocity=(MYTYPE *)calloc(sizeof(MYTYPE), (unsigned)pinit);
	for(i=0;i<pinit;i++)
		*(velocity+i)=0.0;  //粒子速度初始化
	
	MYTYPE gbest,gbestvalue;  //

	MYTYPE *pbest=(MYTYPE *)calloc(sizeof(MYTYPE), (unsigned)pinit);  //建立pBest数组
	MYTYPE *pbestvalue=(MYTYPE *)calloc(sizeof(MYTYPE), (unsigned)pinit);//建立pBestValue数组
	MYTYPE *pnowvalue=(MYTYPE *)calloc(sizeof(MYTYPE), (unsigned)pinit);//建立pNowValue数组
	
	//进入计算主循环
	do
	{
		for(i=0;i<pinit;i++)//分粒子计算
		{
			*(particle+i)=*(particle+i)+*(velocity+i);//利用速度和前一次位置计算当前位置
			
			*(pnowvalue+i)=(* function)(*(particle+i));//计算当前位置待求函数值

			if(num==0)
			{
				*(pbest+i)=*(particle+i);
				*(pbestvalue+i)= (* function)(*(particle+i));
			}//初始化当前粒子最优值
						
			if((*cmpbetter)(*(pnowvalue+i),*(pbestvalue+i)))
			{
				*(pbest+i)=*(particle+i);
				*(pbestvalue+i)=*(pnowvalue+i);
			}//当前函数值优于该粒子最优值
		}

		if(FLAG==0)
		{
			gbest=*(particle);
			gbestvalue=fun(gbest);
			FLAG=1;
		}//初始化全局最优值

        for(i=0;i<pinit;i++)
		{
			if((*cmpbetter)(*(pbestvalue+i),gbestvalue))
			{
				gbest=*(pbest+i);
				gbestvalue=*(pbestvalue+i);
			}
		}

		for(i=0;i<pinit;i++)
		{
			*(velocity+i)=get_rand(0.0,1.0)*(gbest-(*(particle+i)))-0.5*get_rand(0.0,1.0)*(*(pbest+i)-(*(particle+i)));
		}//Calculate the velocity

		if(num<10)
		{
			printf("\n the %dth time pbest:\n",num);
			prin(pinit,pbest,"pbest");
		}
	
		num++;

	}while(!itok(num));

	printf("\nrot times:%d",num);
	printf("\ngbest=%f\n",gbest);
	prin(pinit,pbest,"pbest");

	free(particle);
	free(velocity);
	free(pbest);
	free(pbestvalue);
	free(pnowvalue);

	return gbest;
}
/*------------------------------------------------------------------*/

MYTYPE getmnrand(MYTYPE i,MYTYPE j)
{
	return (MYTYPE)get_rand((MYTYPE)i,(MYTYPE)j);
}

double get_rand(double low_end, double high_end) 
{ 
    /* Assume rand() has been seeded already with srand */ 
    int num; 
    double retval; 

    num = rand(); /* Gives a number between 0 and RAND_MAX */ 

    retval = ((double) num / (double) (RAND_MAX-1)) * (high_end - 
              low_end) + low_end; 

    return retval; 
} 

/*------------------------------------------------------------------*/
int itok(int n)
{
	if (n>KMAX)
		return 1;
	else 
		return 0;
}

/*------------------------------------------------------------------*/
int comparebt(MYTYPE z1,MYTYPE z2)
{
	if(z1<z2)	
		return 1;
	else
		return 0;
}

/*------------------------------------------------------------------*/
MYTYPE fun(MYTYPE x)
{
	return (x+1)*cos(0.1*x);
}

MYTYPE fun1(MYTYPE x)
{
	return (x-50.0)*(x-60.0)+20.0;
}

/*------------------------------------------------------------------*/
double abbs(MYTYPE i)
{
	if(i<0.0)
		return -i;
	else
		return i;
}

/*------------------------------------------------------------------*/
void prin(int i,MYTYPE * p, char *s)
{
	int j;
	for(j=0;j<i;j++)
	{
		printf("%s[%d]=%f\n",s,j+1,*(p+j));
	} 
}

/*------------------------------------------------------------------*/