monkey-for-test-3.cpp

猴群算法

/*********************************************************************************************************/
/*  The algorithm on SPSA based on fuzzy simulation 
/*     see example 2 in Liu and Liu (Expected Value of Fuzzy Variable and Fuzzy Expected Value Models)   */
/*     IEEE TRASACTIONS ON FUZZY SYSTEMS,VOL.10,nO.4,AUGUST 2002                                         */
/*  Copyright     Tianjin  University                                                                    */
/*                 2004.5                                                                                */
/*********************************************************************************************************/
#include "math.h"
#include "stdlib.h"
#include "stdio.h"
#include "conio.h"
#include "UTlab.h"


#define O 1        // number of objectives
#define N 10000     //the number of decision variables


#define Bar 0.1   // the length of bar
#define Monk 5      //the number of Monkey

#define	epsilon 0.001         


#define	 TYPE -1    // 1=max; -1=min
#define	 Iter 2

static double OBJECTIVE[Monk+1][O+1];
static double Decision[Monk+1][N+1];
static double Optimal_Decision[N+1];
static double delta[Monk+1][N+1];
static double g[N+1];    //pesugradient
static int interation;

static int    evaluation(int iter); 
static void   Objectives(void);
static int    constraint_check(double x[]);
static void   Init_MC(void);
static void   Climb(void);
static void   Observations(void);
static void  Turn(void);
static void  Cooperation(void);




//以下求目标函数值 
static void Objectives(void)
{
	double x[N+1],y[N+1],z[N+1];
	int i,j;

	for(j=1;j<=Monk;j++)
	{
		OBJECTIVE[j][1]=0;
		y[0]=0;
		z[0]=0;
		for(i=1;i<=N;i++)
		{   
			x[i]=Decision[j][i];
			y[i]=y[i-1]+x[i]*x[i];
			z[i]=z[i-1]+cos(2*3.1415926*x[i]);				
		}
		OBJECTIVE[j][1]=-20*exp(-0.2*pow((1.0/N)*y[N],0.5))-exp((1.0/N)*z[N])+20.0+exp(1.0);
	}
	for(i=1;i<=Monk;i++)
	  OBJECTIVE[i][0]=OBJECTIVE[i][1];
}

//判断约束条件开始
static int constraint_check(double x[])
{
	int i;
	for(i=1;i<=N;i++)
	if(x[i]<-32||x[i]>32)  return 0;
	
    return 1;
}

static void Init_MC(void)
{
	double x[N+1];
	int i,j;

	for(j=1;j<=Monk;j++)
	{
		do
		{	
			for(i=1; i<=N; i++)
			x[i]=myu(-32,32);	
			
		}while( constraint_check(x) == 0 );
		for(i=1;i<=N;i++) 
			Decision[j][i] = x[i];
	}
	 
}


static void Climb(void)
{	
	int   i,j,k=1,p;  
	double xplus[Monk+1][N+1],xminus[Monk+1][N+1],y[N+1],temp[Monk+1][N+1],g[Monk+1][N+1];  
	double ak,ck;
	double yplus[Monk+1],yminus[Monk+1];
	

	while (k<=30 )
	{


		 
		 for(j=1;j<=Monk;j++)
		 {
			 for (i=1;i<=N;i++)
			 {
				 delta[j][i]=myu(0,1);
				 if(delta[j][i]<0.5)   delta[j][i]=-1;
				 else delta[j][i]=1;
			 }
		 }
		 for(j=1;j<=Monk;j++)
			 for (i=1;i<=N;i++)
			 {
		 		temp[j][i]=Decision[j][i];
			 }

		 for(j=1;j<=Monk;j++)
			 for (i=1;i<=N;i++)
			{
				 xplus[j][i]=Decision[j][i]+epsilon*delta[j][i];
				 xminus[j][i]=Decision[j][i]-epsilon*delta[j][i];
			}

		 

		 for(j=1;j<=Monk;j++)
		 {
			 for(i=1;i<=N;i++) 
				 Decision[j][i]=xplus[j][i];
			 Objectives();
			 yplus[j]=OBJECTIVE[j][1];
		 }

		 for(j=1;j<=Monk;j++)
		 {
			 for(i=1;i<=N;i++) 
				 Decision[j][i]=xminus[j][i];
			 Objectives();
			 yminus[j]=OBJECTIVE[j][1];
		 }
		
		 for(j=1;j<=Monk;j++)
		 {
			 for (i=1;i<=N;i++)
			 {
				 g[j][i]=(yplus[j]-yminus[j])*(delta[j][i]);
				 if(g[j][i]>0)  y[i]=temp[j][i]-epsilon;
				 else y[i]=temp[j][i]+epsilon;
				 if(y[i]>32) y[i]=32;
				 if(y[i]<-32) y[i]=-32;
			 }
			 if( constraint_check(y)==0)
			 {
				for (i=1;i<=N;i++)
				{
					if(g[j][i]>0)  y[i]=y[i]+epsilon;
					else y[i]=y[i]-epsilon;
				}
			 }
			 for (i=1;i<=N;i++) Decision[j][i]=y[i];
			 Objectives();
			 OBJECTIVE[j][1];
			 printf("OBJECTIVE[0][1]===%f,OBJECTIVE[Monk][1]=%f,%d\n",OBJECTIVE[0][1],OBJECTIVE[5][1],k);
		 }
		  
		 k=k+1;
	}
}



static int evaluation(int iter)        
{
	float temp;
	int   i, j, k, label;

	Objectives();

	if(iter == 0)  
	{
		for(k=0; k<=O; k++) OBJECTIVE[0][k] = OBJECTIVE[1][k];
		for(j = 1; j <= N; j++) 
		{
			Decision[0][j] = Decision[1][j];
		}
	}

	for(i=0; i<= Monk; i++)
	{
		label = 0;  temp = OBJECTIVE[i][0];
		for(j=i+1; j<=Monk; j++)
			if((TYPE*temp)<(TYPE*OBJECTIVE[j][0]))
			{
				temp = OBJECTIVE[j][0];
				 label = j;
			}
		if(label != 0) 
		{
			 for(k=0; k<=O;k++) 
			 {
				  temp = OBJECTIVE[i][k];
				  OBJECTIVE[i][k] = OBJECTIVE[label][k];
				  OBJECTIVE[label][k] = temp;
			 }
			 for(j = 1; j <= N; j++) 
			 {
				  temp = Decision[i][j];
				  Decision[i][j] = Decision[label][j];
				  Decision[label][j] = temp;
			 }
		 }
	}
	return 1;
}


static void Observations(void)
{

	int i,j,l;
	double y[N+1], temp[Monk+1][N+1], temp1[Monk+1][O+1],Bar1;

	for(j=1;j<=Monk;j++)
	{
		l=0;
		Bar1=Bar;
		for(i=1;i<=N;i++) 
		{
			temp[j][i]=Decision[j][i];
			temp1[j][1]=OBJECTIVE[j][1];

		}
mark1:
		
		for(i=1;i<=N;i++) 
		{
			y[i]=myu(Decision[j][i]-Bar1,Decision[j][i]+Bar1);
			if(y[i]>32) y[i]=32;
				 if(y[i]<-32) y[i]=-32;

		}
		if  (constraint_check(y) == 0) goto mark1; 
	
		
		for(i=1;i<=N;i++) Decision[j][i]=y[i];
		Objectives();
		OBJECTIVE[j][1];			
		if(TYPE*OBJECTIVE[j][1]<TYPE*temp1[j][1]) 
		{   
			l++;
			if(l<5) 
			{
				goto mark1;
			}
			else if(5<=l && l<=10)
			{
				Bar1=Bar1+0.001;
				goto mark1;
			}
			else if(l>10)
			{
				for(i=1;i<=N;i++) 
					Decision[j][i]=temp[j][i];
			}
			OBJECTIVE[j][1]=temp1[j][1];					
		}
		

	}

}


static void  Turn(void)
{
	int   i, j, k;
	double x[N+1], y[N+1],beta;
	
	
	for(k = 1; k <= Monk; k++)
	{
		

		for(i = 1; i <= N; i++) 
		{
			x[i]=0;
			for(j=1;j<=Monk;j++)
			{
				x[i]=x[i]+Decision[j][i];
			}
			x[i]=x[i]/(Monk);			
		}
		do
		{
			for(i=1;i<=N;i++)
			{
				beta=myu(-1,1);
				y[i]=(x[i]-Decision[k][i]/(Monk)+beta*fabs(x[i]-Decision[k][i]/(Monk)-Decision[k][i]));
				if(y[i]>32) y[i]=32;
				 if(y[i]<-32) y[i]=-32;

			}
		}while(constraint_check(y) == 0);
		for(i=1;i<=N;i++) 
		{
			Decision[k][i]=y[i];	

		}
	}	
}

static void  Cooperation(void)
{
	int  i, j;
	double beta,x[N+1];

mark2:
	for(i=1;i<=N;i++) 
	{
		for(j = 1; j <= Monk-1; j++)
		{
			beta=myu(0,1);
			Decision[j][i]=beta*Decision[j][i]+(1-beta)*Decision[j+1][i];
		}
		Decision[Monk][i]=1.0/2.0*(Decision[Monk][i]+Decision[1][i]);
	}
	for(j = 1; j <= Monk; j++)
	{
		for(i=1;i<=N;i++) 
		{
			x[i]=Decision[j][i];
		}
		if(constraint_check(x) == 0) goto mark2;
		break;
		for(i=1;i<=N;i++) Decision[j][i]=x[i];
	}
	

	
}

int  main( void )
{
	int   i,j,l,n,z,k; 
	double mmm[21];
	FILE *fp;
//	fp=fopen("RESULT.txt","w");
	for(k=1;k<=20;k++){
	srand(90+10*k);
	
	Init_MC();
	evaluation(0);
	interation=1;
	do
	{ 
		for(n=1;n<=Iter;n++)
		{			
			Climb();	
			evaluation(n);
			Observations();	
		
		}
		Turn();
//		Cooperation();

		printf("\n Iter = %d\n ", interation);
		for(i=1;i<=N;i++) printf("%d\n  %3.4f ",k,Decision[0][i]);
		printf("YYY ");
		for(i=0;i<=O;i++) printf("%3.4f ", OBJECTIVE[0][i]);
		printf("\n\n");
	//	fprintf(fp,"%f  %f\n",interation/60.0,OBJECTIVE[0][i]/4.5);

//	    for(i=1;i<=N;i++) fprintf(fp,"%10.4f ",Decision[0][i]);

		interation++;
	}while(interation<601);
	mmm[k]=OBJECTIVE[0][i];
		
        fp=fopen("RESULT.txt","w");
	    for(z=1;z<=k;z++){
			fprintf(fp,"%d  %f\n",90+10*z,mmm[z]/4.5);
		}
		fclose(fp);
	}
	return 1;
	
}