pso_basic.cpp

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

/* Basic PSO
[e-mail]
Last update : 2004-08-09

Constriction coefficient
D-rectangular DPNP (distributions of possible next positions)
Random information links
Interval confinement
Granularity confinement
Pseudo-gradient option (to find the best informer)
Multiobjective (results are saved on a file. You will have then to
                "clean up" them by using Pareto dominance)
Parameter values and function choice are hard coded.
A good pseudo random number generator is included (KISS), but you may replace it.

Note: the original version was in French. Hence the variable names ;-)
 */
#include "stdio.h"
#include "math.h"	
#include <stdlib.h>

#define	D_max 35  // Max number of dimensions of the search space
#define	S_max 100 // Max swarm size
#define F_max 5 // (for multiobjective) Max number of functions
#define ulong unsigned long 	// For KISS pseudo random numbers generator
#define RAND_MAX_KISS ((unsigned long) 4294967295) // idem
					// You don't need KISS if you have a good rand() function

// Structures
struct f	{int taille;double fi[F_max];};
struct	vecteur	{int taille;double v[D_max];};
struct	position {int taille;double x[D_max]; struct f f;};

// Sub-programs 
double	alea(double a,double b); 
int		alea_entier(int a,int b);
double alea_normal(double mean, double std_dev);
void	confin_granul(int s);
void	confin_interv(int s);
double	distance(int s1,int s2);
double	erreur_totale(struct f f);
double	max(double a, double b);
int		meilleur(struct f fx,struct f fy);
int		meilleure_info(int s,int option);
double	min(double a, double b);
double	perf(struct position x,int fonction);
ulong	rand_kiss();
void	seed_rand_kiss(ulong seed);

// Global variables
int	D;				// Search space dimension
double  E; // Useful for some test functions
struct f	fmin;	// Objective(s) to reach
double	granul[D_max];	// Granularity for each dimension (0 = continue, 1 = integer)
int		nb_eval; // Total number of evaluations
int	LIENS[S_max][S_max]; // Information links

struct	position meilleure;	// Best of the best position 
struct	position P[D_max];	// Positions
struct	position P_m[D_max];// Best positions found by each particle
int	S;				// Swarm size
double pi; // Useful for some test functions
struct	vecteur V[D_max];// Velocities
double	xmin[D_max],xmax[D_max];	// Intervals defining the search space

// File(s)
FILE	*f_multi;	// In case of multiobjective, save results, in order to later find
					// the Pareto front
// =================================================
void main()
{
double	c1,cmax;	// Confidence coefficient
int	d;				// Current dimension
double	eps;		// Wanted accuracy
double	eps_moyen;	// Average accuracy 
double	erreur;		// Total error for a given position
double	erreur_prec; // Total error after previous iteration
int	eval_max;		// Max number of evaluations

double	eval_moyen;	// Mean number of evaluations
int	fonction[F_max];
int	g;				// Rank of the best informer
int	i;
int	init_liens;		// Flag to (re)init or not the information links
int	k;
int	K;				// Max number of informed particles by a given one
int	m;
double min; // Best result along several runs
double	phi;		// Control exploration/exploitation
int		n_f;		// (multiobjective) number of functions to optimize
int	n_echec;		// Number of failuers
int	n_exec,n_exec_max; // Nbs of executions
int	option_meilleure; // Option for the best informer choice
						// 0 => really the best (classical)
						// 1 => using pseudo-gradient
int	s;				// Rank of the current particle


seed_rand_kiss(1); // Initialize KISS
E=exp(1);
pi=acos(-1);

// Control parameters
phi=4.14; //4.14 is usually a good compromise exploration/ exploitation
// Confidence coefficients are computed according to
// a theoretical constriction coefficient
c1=2/(phi-2+sqrt(phi*phi-4*phi));
cmax=c1*phi/2;
S=20;	// Swarm size. Typically 20, but may be far smaller or easy problems
K=3;	// Max number of informed particles by a given one
      // Either 3 or S (thumb rule)
// Information links topology is partly random
// Proximity distributions are classical D-rectangles

//---------------------------------------------------------------
n_f=1;			// Nb of functions to optimize simultaneously
fonction[0]=15; //Function code (see perf())
fonction[1]=12; // Test multiobjective for n_f=2 and fonction[0]=11
/*
 0 Parabola (Sphere)
 1 Parabola (Sphere) with shift
 2 Griewank
 3 Rosenbrock
 4 Step
 5 Stochastic
 6 Foxholes 2D
 7 Polynomial fitting problem 9D
 8 Alpine
 9 Rastrigin
10 Ackley
11 Multiobjective Lis-Eiben 1, function 1
12 Multiobjective Lis-Eiben 1, function 2
13 Tripod 2D
14 Pressure Vessel
15 Coil Compression Spring 
*/

D=3;	// Search space dimension 

/* // D-cube data
for (d=0;d<D;d++) // Intervals and granularities
{xmin[d]=-100; xmax[d]=100; granul[d]=0;}
*/

/* // Data for Pressure Vessel
xmin[0]=1.1; xmax[0]=12.5; granul[0]=0.0625;
xmin[1]=0.6; xmax[1]=12.5; granul[1]=0.0625;
xmin[2]=0.001; xmax[2]=240; granul[2]=0;
xmin[3]=0.001; xmax[3]=240; granul[3]=0;
*/

 // Data for (continuous) Coil Compression Spring
xmin[0]=1; xmax[0]=70; granul[0]=0;
xmin[1]=0.6; xmax[1]=3; granul[1]=0;
xmin[2]=0.2; xmax[2]=1; granul[2]=0;


// Save if multiobjective
if (n_f>1)
	f_multi=fopen("f_multi.txt","w"); 

eps=0.00001; 		// Accuracy
fmin.taille=n_f;
for (i=0;i<n_f;i++)
	fmin.fi[i]=0.0;	// Objective(s) to reach
n_exec_max=5; 		// Numbers of runs
eval_max=94000; 	// Max number of evaluations
option_meilleure=0;	// 1 or 0.Using pseudo-gradient or not

// Initialisation of information variables
n_exec=0; eval_moyen=0; eps_moyen=0; n_echec=0;

init: // Initialisations of positions and velocities
n_exec=n_exec+1;
for (s=0;s<S;s++) // Loop on particles
{
	P[s].taille=D; 
	for (d=0;d<D;d++) P[s].x[d]=alea(xmin[d],xmax[d]);	
	V[s].taille=D; 
	for (d=0;d<D;d++) 
		V[s].v[d]=alea((xmin[d]-xmax[d])/2,(xmax[d]-xmin[d])/2);
	confin_granul(s);
}

// First evaluations
nb_eval=0;

for (s=0;s<S;s++)
{	 
	P[s].f.taille=n_f;
	for (i=0;i<n_f;i++)
		P[s].f.fi[i]=fabs(perf(P[s],fonction[i])-fmin.fi[i]); 
	P_m[s]=P[s]; // Best position = current one
}

// Save the best of the bests
meilleure=P_m[0];
for (s=0;s<S;s++) if (meilleur(P_m[s].f,meilleure.f)==1) meilleure=P_m[s];

// Current min error
erreur_prec=erreur_totale(meilleure.f);
init_liens=1; // So that information links will be initialized

boucle:
if(init_liens==1)
{
	// Who informs who, at random
	for (s=0;s<S;s++) // 
	{
		for (m=0;m<S;m=m++) LIENS[m][s]=0; 
		LIENS[s][s]=1; // However, each particle informs itself
	}

	for (m=0;m<S;m=m++) // Other links
	{
		for (k=0;k<K;k++){s=alea_entier(0,S-1);LIENS[m][s]=1;}
	}
}

// The swarm MOVES
for (s=0;s<S;s++)  // For each particle ...
{
	g=meilleure_info(s,option_meilleure); // .. find the best informer
	// ... compute new velocity
	for (d=0;d<D;d++) 
	{
		V[s].v[d]=c1*V[s].v[d]+alea(0,cmax)*(P_m[s].x[d]-P[s].x[d]);
		V[s].v[d]=V[s].v[d]+alea(0,cmax)*(P_m[g].x[d]-P[s].x[d]);
	}
	// ... and move
	for (d=0;d<D;d++) P[s].x[d]=P[s].x[d]+V[s].v[d];

	// ... interval confinement
	confin_interv(s);

	// ... granularity confinement
	confin_granul(s);
	
	// ... evaluate the new position
	for (i=0;i<n_f;i++) P[s].f.fi[i]=fabs(perf(P[s],fonction[i])-fmin.fi[i]);

	// ... update the best position
	if (meilleur(P[s].f,P_m[s].f)==1) P_m[s]=P[s];
	
	// ... update the best of the bests
	if (meilleur(P_m[s].f,meilleure.f)) meilleure=P_m[s];

}

// Check if finished
erreur=erreur_totale(meilleure.f);
// If no improvement, information links will be reinitialized
if(erreur>=erreur_prec) init_liens=1; else init_liens=0;
erreur_prec=erreur;

if (erreur>eps && nb_eval<eval_max) goto boucle;
if (erreur>eps) n_echec=n_echec+1;

// Result display
printf("\n\nExec %i Eval= %i. Value(s) ",n_exec,nb_eval);
for (i=0;i<n_f;i++) printf(" %f",meilleure.f.fi[i]);
printf("\n Position :\n");
for (d=0;d<D;d++) printf(" %f",meilleure.x[d]);

if (n_exec==1) min=erreur; else if(erreur<min) min=erreur;

// Multiobjective result save
if (n_f>1)
{
	for (i=0;i<n_f;i++)
		fprintf(f_multi,"%f ",meilleure.f.fi[i]);
	for (d=0;d<D;d++) fprintf(f_multi," %f",meilleure.x[d]);
	fprintf(f_multi,"\n");
}

// Compute and display some statistical information
eval_moyen=eval_moyen+nb_eval;
eps_moyen=eps_moyen+erreur;

if (n_exec<n_exec_max) goto init;

eval_moyen=eval_moyen/(double)n_exec;
eps_moyen=eps_moyen/(double)n_exec;
printf("\n\n Eval. (mean)= %f",eval_moyen);
printf("\n Error (mean) = %f",eps_moyen);
printf("\n Success rate = %.2f%%",100*(1- n_echec/(double)n_exec));
if (n_exec>1) printf("\n Best min value = %f",min);
}

//===========================================================
double	alea(double a,double b)
{	// rand number (uniform distribution) in [a b]
double 	r;
 r=(double)rand_kiss()/RAND_MAX_KISS;
// r=(double)rand(); r=r/RAND_MAX;
return a+r*(b-a);
}
//===========================================================
int	alea_entier(int a,int b)
{// Integer rand 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;
          y1=y1*std_dev+mean;
         return y1;
  }
  
//===========================================================
void confin_granul(int s)
{ // Granularity confinement
	int d;
  double xd;

for (d=0;d<D;d++)
	{
    if (granul[d]>0)
    {
    xd= P[s].x[d];
    V[s].v[d]=V[s].v[d]-xd;
    if (xd>=0) 			P[s].x[d]=granul[d]*floor (xd/granul[d]+0.5);
            //"1/granul" is the minimum distance between two points

    else  P[s].x[d]=-granul[d]*floor(-xd/granul[d]+0.5);
	V[s].v[d]=V[s].v[d]+P[s].x[d];
    }
 }
}

//===========================================================
void confin_interv(int s)
{ // Interval confinement (keep the particle "in the box")
	int d;
	for (d=0;d<D;d++)
	{
		if (P[s].x[d]<xmin[d]) {P[s].x[d]=xmin[d];V[s].v[d]=0;}
		if (P[s].x[d]>xmax[d]) {P[s].x[d]=xmax[d];V[s].v[d]=0;}
	}
}

//===========================================================
double	distance(int s1,int s2)
{
// Euclidean distante between the best memorized positions
// of particles s1 and s2
int	d;
double	dist;
double	z;

dist=0;

for (d=0;d<P_m[s1].taille;d++)
{ 	
	z=P_m[s1].x[d]-P_m[s2].x[d]; 
	dist=dist+z*z;}
	return sqrt(dist);
}

//===========================================================
double	erreur_totale(struct f f)
{
	double	err=0;
	int		i;
	double	z;

	for (i=0;i<f.taille;i++)
	{
		z=f.fi[i]; err=err+z*z;
	}
	return sqrt(err);
}

//===========================================================
double max(double a, double b)
{
	if (a>b) return a; else return b;
}

//===========================================================
int	meilleur(struct f fx,struct f fy)
{
	/*