standard_pso_2007.c

用c语言实现标准的微粒群算法

	// fprintf( f_stag, "\nEND" );
	R.error = error;
	return R;  
}

  
    // ===========================================================
double alea (double a, double b) 
{				// random number (uniform distribution) in  [a b]
  double r;
  r = (double) rand ();
  r = r / RAND_MAX;
  return a + r * (b - a); 
}
  
// ===========================================================
int alea_integer (int a, int b) 
{				// Integer random number in [a b]
  int ir;
  double r;
    
	r = alea (0, 1);
  ir = (int) (a + r * (b + 1 - a));
    
	if (ir > b)	ir = b;
    
	return ir;  
}
  
    // ===========================================================
double alea_normal (double mean, double std_dev) 
{ 
  /*
  Use the polar form of the Box-Muller transformation to obtain a pseudo
	random number from a Gaussian distribution 
  */ 
  double x1, x2, w, y1;  
      // double y2;
      
  do  
  {
		x1 = 2.0 * alea (0, 1) - 1.0;
		x2 = 2.0 * alea (0, 1) - 1.0;
		w = x1 * x1 + x2 * x2;     
	}
  while (w >= 1.0);
    
	w = sqrt (-2.0 * log (w) / w);
  y1 = x1 * w;
  // y2 = x2 * w;
  if(alea(0,1)<0.5) y1=-y1; 
  y1 = y1 * std_dev + mean;
  return y1;  
}

	
// =============================================================
struct velocity aleaVector(int D,double coeff)
{
	struct velocity V;
	int d;
	int i;
	int K=2; // 1 => uniform distribution in a hypercube
			// 2 => "triangle" distribution
	double rnd;
	
	V.size=D;
	
	for (d=0;d<D;d++)
	{								
		rnd=0;
		for (i=1;i<=K;i++) rnd=rnd+alea(0,1);		
		V.v[d]=rnd*coeff/K;
	}
	
	return V;
}
    
      // ===========================================================
double normL (struct velocity v,double L) 
{   // L-norm of a vector
	int d;     
	double n;
      
	n = 0;
      
	for (d = 0; d < v.size; d++)
		n = n + pow(fabs(v.v[d]),L);
      
	n = pow (n, 1/L);     
	return n;  
}
    
// ===========================================================
double distanceL (struct position x1, struct position x2,double L) 
{  // Distance between two positions
	// L = 2 => Euclidean	
	int d;     
	double n;
      
	n = 0;
      
	for (d = 0; d < x1.size; d++)
		n = n + pow (fabs(x1.x[d] - x2.x[d]), L);
      
	n = pow (n, 1/L);
  return n;    
}

//============================================================
struct matrix matrixProduct(struct matrix M1,struct matrix M2)
{
	// Two square matrices of same size
	struct matrix Product;
	int D;	
	int i,j,k;
	double sum;
	D=M1.size;
	for (i=0;i<D;i++)
	{
		for (j=0;j<D;j++)
		{			
			sum=0;
			for (k=0;k<D;k++)
			{
				sum=sum+M1.v[i][k]*M2.v[k][j];	
			}
			Product.v[i][j]=sum;			
		}	
	}
	Product.size=D;
	return Product;
}	
//=========================================================
struct matrix matrixRotation(struct velocity V)
{
	/*
		Define the matrice of the rotation V' => V
	where V'=(1,1,...1)*normV/sqrt(D)  (i.e. norm(V') = norm(V) )
	
	*/
	struct velocity B;
	int i,j,d, D;
	double normB,normV,normV2;
	//struct matrix reflex1; // Global variable
	struct matrix reflex2;
	struct matrix rotateV;
	double temp;
	
	D=V.size;
	normV=normL(V,2); normV2=normV*normV;
	reflex2.size=D;
	
	// Reflection relatively to the vector V'=(1,1, ...1)/sqrt(D)	
	// norm(V')=1
	// Has been computed just once  (global matrix reflex1)	
	
	//Define the "bisectrix" B of (V',V) as an unit vector
	B.size=D;
	temp=normV/sqrtD;
					
	for (d=0;d<D;d++)
	{
		B.v[d]=V.v[d]+temp;
	}
	normB=normL(B,2);
			
	if(normB>0)
	{
		for (d=0;d<D;d++)
		{
			B.v[d]=B.v[d]/normB;
		}
	}
					
	// Reflection relatively to B
	for (i=0;i<D;i++)
	{
		for (j=0;j<D;j++)
		{
			reflex2.v[i][j]=-2*B.v[i]*B.v[j];
		}
	}
							
	for (d=0;d<D;d++)
	{
		reflex2.v[d][d]=1+reflex2.v[d][d];
	}
						
	// Multiply the two reflections
	// => rotation				
	rotateV=matrixProduct(reflex2,reflex1);
	return rotateV;
	
}
//==========================================================
struct velocity	matrixVectProduct(struct matrix M,struct velocity V)
{
	struct velocity Vp;
	int d,j;
	int Dim;
	double sum;
	Dim=V.size;
	for (d=0;d<Dim;d++)
	{
		sum=0;
		for (j=0;j<Dim;j++)
		{
			sum=sum+M.v[d][j]*V.v[j];	
		}
		Vp.v[d]=sum;
	}	
	Vp.size=Dim;
	return Vp;
}
// ===========================================================
int sign (double x) 
{     
	if (x == 0)	return 0;
  	if (x < 0)	return -1;    
			return 1;   
}
    
// ===========================================================
struct position quantis (struct position x, struct SS SS) 
{     
	/*
	  Quantisatition of a position
	Only values like x+k*q (k integer) are admissible 
	 */ 
  int d;
  double qd;
  struct position quantx;

	quantx = x;     
	for (d = 0; d < x.size; d++)
	{
	    qd = SS.q.q[d];	
	  
		if (qd > zero)	// Note that qd can't be < 0
	  {     
			qd = qd * (SS.max[d] - SS.min[d]) / 2;	      
			quantx.x[d] = qd * floor (0.5 + x.x[d] / qd);	    
		}
	}
	return quantx;    
}
// ===========================================================
double perf (struct position x, int function, struct SS SS) 
{				// Evaluate the fitness value for the particle of rank s   
	int d;
	double DD;
	int  k;
	double f, p, xd, x1, x2;
	double s11, s12, s21, s22;
	double sum1,sum2;
	double t0, tt, t1;
	struct position xs; 
	// Shifted Parabola/Sphere (CEC 2005 benchmark)		
	static double offset_0[30] =
{ 
-3.9311900e+001, 5.8899900e+001, -4.6322400e+001, -7.4651500e+001, -1.6799700e+001,
-8.0544100e+001, -1.0593500e+001, 2.4969400e+001, 8.9838400e+001, 9.1119000e+000, 
-1.0744300e+001, -2.7855800e+001, -1.2580600e+001, 7.5930000e+000, 7.4812700e+001,
 6.8495900e+001, -5.3429300e+001, 7.8854400e+001, -6.8595700e+001, 6.3743200e+001, 
 3.1347000e+001, -3.7501600e+001, 3.3892900e+001, -8.8804500e+001, -7.8771900e+001, 
-6.6494400e+001, 4.4197200e+001, 1.8383600e+001, 2.6521200e+001, 8.4472300e+001
};	
	xs = x;

	switch (function)
	{
		
		case 100:
		for (d = 0; d < xs.size; d++) 
	  	{
			xs.x[d]=xs.x[d]-offset_0[d];
		}
		
		case 0:		// Parabola (Sphere)
	  	f = 0;

		for (d = 0; d < xs.size; d++) 
	  	{    
			xd = xs.x[d];      
			f = f + xd * xd;    
		}	  
		break;
	
		case 1:		// Griewank
	  f = 0; 
		p = 1;
	  
		for (d = 0; d < xs.size; d++)
	  {      
			xd = xs.x[d];	      
			f = f + xd * xd;	      
			p = p * cos (xd / sqrt ((double) (d + 1)));	    
		} 
		f = f / 4000 - p + 1;	  
		break;
	
		case 2:		// Rosenbrock
		
	  f = 0;  
		t0 = xs.x[0]  + 1;	// Solution on (0,...0) when
	  															// offset=0
	  for (d = 1; d < xs.size; d++)
	  {     
			
			t1 = xs.x[d]  + 1;	      
			tt = 1 - t0;	      
			f += tt * tt;      
			tt = t1 - t0 * t0;      
			f += 100 * tt * tt;	      
			t0 = t1;    
		}  
		break;
	
		case 3:		// Rastrigin
	  k = 10;  
		f = 0;
	  
		for (d = 0; d < xs.size; d++)    
	  {     
			xd = xs.x[d];	      
			f =f+ xd * xd - k * cos (2 * pi * xd);	    
		}	  
		f =f+ xs.size * k;  
		break;
	
		case 4:		// 2D Tripod function
	  // Note that there is a big discontinuity right on the solution
	  // point. 
		
	x1 = xs.x[0] ; 
	x2 = xs.x[1];  
	s11 = (1.0 - sign (x1)) / 2;
	s12 = (1.0 + sign (x1)) / 2; 
	s21 = (1.0 - sign (x2)) / 2;
	s22 = (1.0 + sign (x2)) / 2;
	  
	//f = s21 * (fabs (x1) - x2); // Solution on (0,0)
	f = s21 * (fabs (x1) +fabs(x2+50)); // Solution on (0,-50)  
	f = f + s22 * (s11 * (1 + fabs (x1 + 50) +
			fabs (x2 - 50)) + s12 * (2 +
			fabs (x1 - 50) +
			fabs (x2 - 50)));	  
	break;
	case 5:  // Ackley
	sum1=0;
	sum2=0;
	DD=x.size;
	pi=acos(-1);
	for (d=0;d<x.size;d++)
		{
		xd=xs.x[d];
		sum1=sum1+xd*xd;
		sum2=sum2+cos(2*pi*xd);
		}
	f=-20*exp(-0.2*sqrt(  sum1/DD  ))-exp(sum2/DD)+20+exp(1);
	
	break;
		
	case 99: // Test
	xd=xs.x[0];
	//if (xd<9) f=10-xd; else f=10*xd-90;	
	if (xd<=1) f=10*xd; else f=11-xd;	
	break;
	
	}
      
return f;    
}

//===================================================
struct problem problemDef(int functionCode)
{
	int d;
	struct problem pb;
	
	pb.function=functionCode;	
	pb.epsilon = 0.0001;	//0.0001 Acceptable error
	pb.objective = 0;       // Objective value
	
	// Define the solution point, for test
	// NEEDED when param.stop = 2 
	// i.e. when stop criterion is distance_to_solution < epsilon
	for (d=0; d<30;d++)
	{
		pb.solution.x[d]=0;
	}
	
	
	// ------------------ Search space
	switch (pb.function)
	{   
		case 0:			// Parabola
		case 100:
		pb.SS.D =30;// 30 	// Dimension																								// values
		
		for (d = 0; d < pb.SS.D; d++)
		{   
			pb.SS.min[d] = -100; // -100
			pb.SS.max[d] = 100;	// 100
			pb.SS.q.q[d] = 0;	// Relative quantisation, in [0,1].   
		 }
				
		pb.evalMax = 6000;// 6000	// Max number of evaluations for each run
		break;
				
		case 1:		// Griewank
		pb.SS.D = 30;	// Dimension   
	
		 // Boundaries
		 for (d = 0; d < pb.SS.D; d++) 
		 {	
				pb.SS.min[d] = -100;
				pb.SS.max[d] = 100;
				pb.SS.q.q[d] = 0;
			}
				
		pb.evalMax = 9000;	    
		break;
				
		case 2:		// Rosenbrock
		pb.SS.D = 30;	// Dimension
				
					// Boundaries
		for (d = 0; d < pb.SS.D; d++)
		{	
			pb.SS.min[d] = -10; pb.SS.max[d] = 10;			
			pb.SS.q.q[d] = 0;	      
		}
				
		pb.evalMax =40000; // 40000	    
		break;
				
	
		case 3:		// Rastrigin
		pb.SS.D = 30;	// Dimension
				
					// Boundaries
		for (d = 0; d < pb.SS.D; d++)
		{
			pb.SS.min[d] =-10; 
			pb.SS.max[d] =10; 	  
			pb.SS.q.q[d] = 0;	
		 }
				
		pb.evalMax = 40000; 	    
		break;
				
		case 4:		// Tripod
		pb.SS.D = 2;	// Dimension
	
					// Boundaries
		for (d = 0; d < pb.SS.D; d++)
		{
			pb.SS.min[d] = -100;
			pb.SS.max[d] = 100;
			pb.SS.q.q[d] = 0;	
		 }
				
		pb.evalMax = 10000; 	    
		break;
		case 5: // Ackley
		pb.SS.D = 10;	// Dimension
		// Boundaries
		for (d = 0; d < pb.SS.D; d++)
		{
			pb.SS.min[d] = -100; // -32
			pb.SS.max[d] = 100; // 32
			pb.SS.q.q[d] = 0;	
		 }
		pb.evalMax = 6000;		 
		break;
		
		 
	// case n: // Add here your own problem
		 
	case 99: // Test
	pb.SS.D = 1;	// Dimension
		// Boundaries
		for (d = 0; d < pb.SS.D; d++)
		{
			pb.SS.min[d] = 0;
			pb.SS.max[d] = 10;
			pb.SS.q.q[d] = 0;	
		 }
		pb.evalMax = 10000;		 
		break;	
			
	}
		
	pb.SS.q.size = pb.SS.D;
	return pb;
}