pso_basic.cpp

这是一个基本的微粒群算法

	Compare two position evaluations. In case of multiobjective
	it is the Pareto dominance
	Return 1 is fx better than fy, 0 else
	  */
	int	i;
	for(i=0;i<fx.taille;i++)
 {
    if (fx.fi[i]>fy.fi[i]) return 0;
 }
 // All fx values are <=
for(i=0;i<fx.taille;i++)
 {
      if (fx.fi[i]<fy.fi[i]) return 1;  // At least one is <
 }
return 0;    //You may also return 2, if this particular case
            // is interesting
}

//===========================================================
int	meilleure_info(int s,int option)
{
/*
Recherche le rang g de la meilleure informatrice de la particule s
option	=	0 => recherche classique
		=	1 => pseudo gradient
Return the rank of the best informer
option	= 0 => classical (really the best)
		= 1 => using pseudo-gradient
*/
double	delta1,delta2;
double	dist;
double	erreur_1,erreur_2;
int	g;
double	grad;
int	m;
struct f	min_f;
int	s2;

switch (option)
{

case 1: // Pseudo gradient
	g=s; // If nothing better, the best informer is
			// the particle itself
	grad=0;
	erreur_1=erreur_totale(P_m[s].f);

	for (s2=0;s2<S;s2++) 
	{
		
		if (LIENS[s2][s]==0) continue; // Pas de lien

		erreur_2=erreur_totale(P_m[s2].f);
		delta1=erreur_1-erreur_2;
		dist=distance(s,s2);
		if (dist<=0) continue;

		delta2=delta1/dist;
		if (delta2>grad)
		{
			g=s2; grad=delta2;
		}
	}
	// if (g==s) goto case0; // Possible variant
	break;


case 0: // Classical option: the "real" best
case0:
	for (g=0;g<S;g++) {if (LIENS[g][s]==0) continue; goto suite;}
suite:	
	min_f=P_m[g].f; 
	for (m=g+1;m<S;m++) 
	{
		if (LIENS[m][s]==0) continue;
			if (meilleur(P_m[m].f,min_f)==1) {g=m;min_f=P_m[m].f;}
	}
break;
}

return g;
}

//===========================================================
double min(double a, double b)
{
	if (a<b) return a; else return b;
}
//===========================================================
double perf(struct position x,int fonction)
{// Evaluate the function value at position x
double c;
int    D, d;
int		i,j,k;
double f,p,xd,x1,x2,x3,x4;
double	sum1,sum2;
double t0, tt, t1;

// For Coil compressing spring
static double	Fmax=1000.0;
static double	S=189000.0;
static double	lmax=14.0;
static double	dmin=0.2;
static double	Dmax=3.0;
static double	Fp=300;
static double	spm=6.0;
static double	sw=1.25;
static double	G=11500000;
double Cf,K,sp,lf;

// For Foxholes problem
static int a[2][25] =
{
   {-32, -16, 0, 16, 32, -32, -16, 0, 16, 32, -32, -16, 0, 16, 32,
      -32, -16, 0, 16, 32, -32, -16, 0, 16, 32  },
   {-32, -32, -32, -32, -32, -16, -16, -16, -16, -16,
      16, 16, 16, 16, 16, 32, 32, 32, 32, 32 }
};
// For polynomial fitting problem
   int const M=60;
   double py, y=-1, dx=(double)M;

nb_eval=nb_eval+1;

D=x.taille;

switch (fonction)
{
case 0: // Parabola (Sphere)
	f=0;
	p=0; // Shift
	for(d=0;d<D;d++) f=f+(x.x[d]-p)*(x.x[d]-p);
break;

case 1: // Parabola (Sphere)
		// (for multibojective test)
	f=0;
	p=2.1; // Shift
	for(d=0;d<D;d++) f=f+(x.x[d]-p)*(x.x[d]-p);
break;

case 2: // Griewank
	f=0;
	p=1;
	for (d=0;d<D;d++)
	{
		xd=x.x[d]-100;
		f=f+xd*xd;
		p=p*cos(xd/sqrt(d+1));
	}
	f=f/4000 -p +1;
	break;

case 3: // Rosenbrock	
	f=0;
   t0=x.x[0];
   for ( d=1; d < D; d++) 
   {
      t1 = x.x[d];
      tt = 1-t0;
      f += tt*tt;
      tt = t1-t0*t0;
      f += 100*tt*tt;     
      t0 = t1;
   } 
   break;

case 4: // Step
	f=0;
	for ( d=0; d < D; d++) f=f+(int)x.x[d];
	break;

case 5: // Stochastic
	f=0;
	for ( d=0; d < D; d++)
	{
		f=f+(d+1)*pow(x.x[d],4)+alea_normal(0,1);
	}
	break;

case 6: //Foxholes 2D
	f=0;
	for (j=0; j<25; j++)
   		{
		sum1=0;
      	for (d=0; d<2; d++)
      		{
         	sum1 =sum1+ pow (x.x[d] - a[d][j],6);
      		}
      	f=f+1/(j+1+sum1);
   		}
	f=1/(0.002+f);
	break;

case 7:	// polynomial fitting problem 
		// Cf. Differential Evolution Homepage, C code
		// on [-100 100]^9
	f=0;
   dx = 2/dx;
   for (i=0;i<=M;i++)
   {
      py = x.x[0];
      for (d=1;d<D;d++)
      {
	 py = y*py + x.x[d];
      }
      if (py<-1 || py>1) f+=(1-py)*(1-py);
      y+=dx;
   }
   py = x.x[0];
   for (d=1;d<D;d++) py=1.2*py+x.x[d];
   py = py-72.661;
   if (py<0) f+=py*py;
   py = x.x[0];
   for (d=1;d<D;d++) py=-1.2*py+x.x[d];
   py =py-72.661;
   if (py<0) f+=py*py;
   break;

case 8: // Clerc's f1, Alpine function, min 0 
      // Note that there are several global optima
	f=0;
	for( d=0;d<D;d++)
	{
		xd=x.x[d];
		f+=fabs(xd*sin(xd)+0.1*xd);
	}
	break;

 case 9:// Rastrigin. Minimum value 0. Solution (0,0 ...0)
	k=10;
	f=0;
	for (d=0;d<D;d++)
		{
		xd=x.x[ d];
		f+=xd*xd - k*cos(2*pi*xd);
		}
	f+=D*k;
   break;
   
 case 10: // Ackley
	sum1=0;
	sum2=0;
	for (d=0;d<D;d++)
		{
		xd=x.x[ d];
		sum1+=xd*xd;
		sum2+=cos(2*pi*xd);
		}
    y=D;
	f=(-20*exp(-0.2*sqrt(sum1/y))-exp(sum2/y)+20+E);
	break;

case 11: // see also case 12.  Multiobjective Lis-Eiben 1
x1=x.x[ 0];
x2=x.x[ 1];
x1=x1*x1+x2*x2;
f=pow(x1,1.0/8);
break;

case 12: // see also case 11.  Multiobjective Lis-Eiben 1
x1= x.x[ 0]-0.5;
x2= x.x[ 1]-0.5;
x2=x1*x1+x2*x2;
f=pow(x2,1.0/4);
break;

case 13: // 2D Tripod function (Louis Gacogne)
        // min 0 on (0  -50)
x1=x.x[0];
x2= x.x[1];
if(x2<0)
{
	f=fabs(x1)+fabs(x2+50);
}
else
{
	if(x1<0)
		f=1+fabs(x1+50)+fabs(x2-50);
	else
		f=2+fabs(x1-50)+fabs(x2-50);
}
break;

case 14: // Pressure vessel
	 	  // Ref. Onwubolu,G.,New Optimisation Techniques in Engineering p. 638
      // Min value 7197.73 (7019.03 if granularity=0)
/* D=4
  1.1 <= x1 <= 12.5        granularity 0.0625
  0.6 <= x2 <= 12.5         granularity 0.0625
  0 .0 <= x3 <= 240
  0.0 <= x4 <= 240
  constraints
  g1:= 0.0193*x3-x1 <=0
  g2 := 0.00954*x3-x2<=0
  g3:= 750*1728-pi*x3*x3*(x4-(4/3)*x3)  <=0
*/
	x1=x.x[0];
	x2=x.x[1];
	x3=x.x[2];
	x4=x.x[3];

   f=0.6224*x1*x3*x4 + 1.7781*x2*x3*x3 + 3.1611*x1*x1*x4 + 19.84*x1*x1*x3;

  // Constraints
  y=0.0193*x3-x1; if ( y>0) {c= 1+pow(10,10)*y;  f=f*c*c; }
  y=  0.00954*x3-x2; if (y>0) {c=1+y; f=f*c*c;  }
  y = 750*1728-pi*x3*x3*(x4+(4.0/3)*x3);  if (y>0) {c=1+y; f=f*c*c;  }
 break;

case 15: // Coil compression spring
	 	  // Ref. Onwubolu,G.,New Optimisation Techniques in Engineering p. 644
      // D=3
      // Min value 2.614 if granularity=0 (2.656 for the discrete problem,
      // but this basic PSO can't use a data list of possible values
      // (you may use Tribes)
      
  x1=x.x[0];
	x2=x.x[1];
	x3=x.x[2];

	f=pi*pi*x2*x3*x3*(x1+2)*0.25;

	// Constraints
	Cf=1+0.75*x3/(x2-x3) + 0.615*x3/x2;
	K=0.125*G*pow(x3,4)/(x1*x2*x2*x2);
	sp=Fp/K;
	lf=Fmax/K + 1.05*(x1+2)*x3;

	y=8*Cf*Fmax*x2/(pi*x3*x3*x3) -S;
	if (y>0) {c=1+y; f=f*c*c*c;}

	y=lf-lmax;
	if (y>0) {c=1+y; f=f*c*c*c;}

	y=sp-spm;
	if (y>0) {c=1+y; f=f*c*c*c;}
  
	y=sp-Fp/K;
	if (y>0) {c=1+pow(10,10)*y; f=f*c*c*c;}

	y=sw- (Fmax-Fp)/K;
	if (y>0) {c=1+pow(10,10)*y; f=f*c*c*c;}
	 break;

}

return f;
}
//===========================================================
//Pseudo-hasard d閠erministe KISS.
/*

 the idea is to use simple, fast, individually promising
 generators to get a composite that will be fast, easy to code
 have a very long period and pass all the tests put to it.
 The three components of KISS are
        x(n)=a*x(n-1)+1 mod 2^32
        y(n)=y(n-1)(I+L^13)(I+R^17)(I+L^5),
        z(n)=2*z(n-1)+z(n-2) +carry mod 2^32
 The y's are a shift register sequence on 32bit binary vectors
 period 2^32-1;
 The z's are a simple multiply-with-carry sequence with period
 2^63+2^32-1.  The period of KISS is thus
      2^32*(2^32-1)*(2^63+2^32-1) > 2^127
*/


static ulong kiss_x = 1;
static ulong kiss_y = 2;
static ulong kiss_z = 4;
static ulong kiss_w = 8;
static ulong kiss_carry = 0;
static ulong kiss_k;
static ulong kiss_m;


void seed_rand_kiss(ulong seed) 
{
    kiss_x = seed | 1;
    kiss_y = seed | 2;
    kiss_z = seed | 4;
    kiss_w = seed | 8;
    kiss_carry = 0;
}

ulong rand_kiss() 
{
    kiss_x = kiss_x * 69069 + 1;
    kiss_y ^= kiss_y << 13;
    kiss_y ^= kiss_y >> 17;
    kiss_y ^= kiss_y << 5;
    kiss_k = (kiss_z >> 2) + (kiss_w >> 3) + (kiss_carry >> 2);
    kiss_m = kiss_w + kiss_w + kiss_z + kiss_carry;
    kiss_z = kiss_w;
    kiss_w = kiss_m;
    kiss_carry = kiss_k >> 30;
    return kiss_x + kiss_y + kiss_w;
}