standard_pso_2007.c

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

/*
Standard PSO 2007
 Contact for remarks, suggestions etc.:
 [e-mail]
 
 Last update
 2007-12-10 Warning about rotational invariance (valid here only on 2D)
 2007-11-22 stop criterion (option): distance to solution < epsilon
 			and log_progress evaluation
 2007-11-21 Ackley function
 
  -------------------------------- Contributors 
 The works and comments of the following persons have been taken
 into account while designing this standard.  Sometimes this is for 
 including a feature, and sometimes for leaving out one. 
 
 Auger, Anne
 Blackwell, Tim
 Bratton, Dan
 Clerc, Maurice
 Croussette, Sylvain 
 Dattasharma, Abhi
 Eberhart, Russel
 Hansen, Nikolaus
 Keko, Hrvoje
 Kennedy, James 
 Krohling, Renato
 Langdon, William
 Liu, Hongbo 
 Miranda, Vladimiro
 Poli, Riccardo
 Serra, Pablo
 Stickel, Manfred
 
 -------------------------------- Motivation
Quite often, researchers claim to compare their version of PSO 
with the "standard one", but the "standard one" itself seems to vary!
Thus, it is important to define a real standard that would stay 
unchanged for at least one year.
This PSO version does not intend to be the best one on the market
(in particular, there is no adaptation of the swarm size nor of the
coefficients). This is simply very near to the original version (1995),
with just a few improvements based on some recent works.
 --------------------------------- Metaphors
swarm: A team of communicating people (particles)
At each time step
    Each particle chooses a few informants at random, selects the best
	one from this set, and takes into account the information given by
	the chosen particle.
	If it finds no particle better than itself, then the "reasoning" is:
	"I am the best, so I just take my current velocity and my previous
	best position into account" 
----------------------------------- Parameters/Options
clamping := true/false => whether to use clamping positions or not
randOrder:= true/false => whether to avoid the bias due to the loop
				on particles "for s = 1 to swarm_size ..." or not
rotation := true/false => whether the algorithm is sensitive 
				to a rotation of the landscape or not 
You may also modify the following ones, although suggested values
are either hard coded or automatically computed:
S := swarm size
K := maximum number of particles _informed_ by a given one
w := first cognitive/confidence coefficient
c := second cognitive/confidence coefficient
 ----------------------------------- Equations
For each particle and each dimension
Equation 1:	v(t+1) = w*v(t) + R(c)*(p(t)-x(t)) + R(c)*(g(t)-x(t))
Equation 2:	x(t+1) = x(t) + v(t+1)
where
v(t) := velocity at time t
x(t) := position at time t
p(t) := best previous position of the particle
g(t) := best position amongst the best previous positions
		of the informants of the particle
R(c) := a number coming from a random distribution, which depends on c
In this standard, the distribution is uniform on [0,c]
Note 1:
When the particle has no informant better than itself,
it implies p(t) = g(t)
Therefore, Equation 1 gets modified to:
v(t+1) = w*v(t) + R(c)*(p(t)-x(t))
Note 2:
When the "non sensitivity to rotation" option is activated
(p(t)-x(t)) (and (g(t)-x(t))) are replaced by rotated vectors, 
so that the final DNPP (Distribution of the Next Possible Positions)
is not dependent on the system of co-ordinates.
 ----------------------------------- Information links topology 
A lot of work has been done about this topic. The main result is this: 
There is no "best" topology. Hence the random approach used here.  
 ----------------------------------- Initialisation
Initial positions are chosen at random inside the search space 
(which is supposed to be a hyperparallelepiped, and often even
a hypercube), according to a uniform distribution.
This is not the best way, but the one used in the original PSO.
Each initial velocity is simply defined as the half-difference of two
random positions. It is simple, and needs no additional parameter.
However, again, it is not the best approach. The resulting distribution
is not even uniform, as is the case for any method that uses a
uniform distribution independently for each component.
The mathematically correct approach needs to use a uniform
distribution inside a hypersphere. It is not very difficult,
and was indeed used in some PSO versions.  However, it is quite
different from the original one. 
Moreover, it may be meaningless for some heterogeneous problems,
when each dimension has a different "interpretation".
------------------------------------ From SPSO-06 to SPSO-07
The main differences are:
1. option "non sensitivity to rotation of the landscape"
	Note: although theoretically interesting, this option is quite
        computer time consuming, and the improvement in result may
		only be marginal. 
2. option "random permutation of the particles before each iteration"
	Note: same remark. Time consuming, no clear improvement
3. option "clamping position or not"
	Note: in a few rare cases, not clamping positions may induce an
	infinite run, if the stop criterion is the maximum number of 
	evaluations
		
4. probability p of a particular particle being an informant of another
	particle. In SPSO-06 it was implicit (by building the random infonetwork)
	Here, the default value is directly computed as a function of (S,K),
	so that the infonetwork is exactly the same as in SPSO-06.
	However, now it can be "manipulated" ( i.e. any value can be assigned)
	
5. The search space can be quantised (however this algorithm is _not_
   for combinatorial problems)
Also, the code is far more modular. It means it is slower, but easier
to translate into another language, and easier to modify.
 ----------------------------------- Use
 Define the problem (you may add your own one in problemDef() and perf())
 Choose your options
 Run and enjoy!
   
 ------------------------------------ Some results
 So that you can check your own implementation, here are some results.
 
(clamping, randOrder, rotation, stop) = (1,1,1,0) 
Function	Domain			error	Nb_of_eval	Nb_of_runs	Result
Parabola	[-100,100]^30	0.0001	6000		100			 52% (success rate)
	shifted		""			""		""			""			 7%	
	shifted		""			""		7000		""			 100% 
Griewank	[-100,100]^30	0.0001	9000		100			 51%
Rosenbrock	[-10,10]^30		0.0001	40000		50			 31
Rastrigin	[-10,10]^30		0.0001	40000		50			 53 (error mean)
Tripod		[-100,100]^2	0.0001	10000		50			 50%
(clamping, randOrder, rotation, stop) = (1,1,0,0) 
Parabola	[-100,100]^30	0.0001	6000		100			 0.69 (0%) 
	shifted		""			""		""			""			 2.16 (0%)			
Griewank	[-100,100]^30	0.0001	9000		100			 14%
Rosenbrock	[-10,10]^30		0.0001	40000		100			 38.7
Rastrigin	[-10,10]^30		0.0001	40000		100			 54.8 (error mean)
Tripod		[-100,100]^2	0.0001	10000		100			 47%
Because of randomness, and also depending on your own pseudo-random 
number generator, you just have to find similar values, 
not exactly these ones. 
 */  
  
#include "stdio.h"
#include "math.h"
#include <stdlib.h>
#include <time.h>
  
#define	D_max 100		// Max number of dimensions of the search space
#define R_max 200		// Max number of runs
#define	S_max 100		// Max swarm size
#define zero  0			// 1.0e-30 // To avoid numerical instabilities
  // Structures
struct quantum 
{
  double q[D_max];
  int size; 
};

struct SS 
{ 
	int D;
	double max[D_max];
	double min[D_max]; 
	struct quantum q;		// Quantisation step size. 0 => continuous problem
};

struct param 
{
  	double c;		// Confidence coefficient
  	int clamping;	// Position clamping or not
	int K;			// Max number of particles informed by a given one
	double p;		// Probability threshold for random topology	
					// (is actually computed as p(S,K) )
	int randOrder;	// Random choice of particles or not
	int rotation;	// Sensitive to rotation or not
  	int S;			// Swarm size
  	int stop;		// Flag for stop criterion
  	double w;		// Confidence coefficient
};

struct position 
{ 
	double f;  
	int improved;   
	int size;  
	double x[D_max];
};

struct velocity 
{  
	int size;  
	double v[D_max]; 
};
struct problem 
{ 
	double epsilon; 	// Admissible error
	int evalMax; 		// Maximum number of fitness evaluations
	int function; 		// Function code
	double objective; 	// Objective value
						// Solution position (if known, just for tests)	
	struct position solution;
	struct SS SS;		// Search space
};

struct swarm 
{ 
	int best; 					// rank of the best particle
	struct position P[S_max];	// Previous best positions found by each particle
	int S; 						// Swarm size 
	struct velocity V[S_max];	// Velocities
  	struct position X[S_max];	// Positions 
};

struct result 
{  
	double nEval; 		// Number of evaluations  
	struct swarm SW;	// Final swarm
	double error;		// Numerical result of the run
};
struct matrix 	// Useful for "non rotation sensitive" option
{  
	int size;  
	double v[D_max][D_max]; 
};

  // Sub-programs
double alea (double a, double b);
int alea_integer (int a, int b);
double alea_normal (double mean, double stdev);
double distanceL(struct position x1, struct position x2, double L);
struct velocity aleaVector(int D,double coeff);
struct matrix matrixProduct(struct matrix M1,struct matrix M2);
struct matrix matrixRotation(struct velocity V);
struct velocity	matrixVectProduct(struct matrix M,struct velocity V);
double normL (struct velocity v,double L);
double perf (struct position x, int function,struct SS SS);	// Fitness evaluation
struct position quantis (struct position x, struct SS SS);
struct problem problemDef(int functionCode);
struct result PSO ( struct param param, struct problem problem);
int sign (double x);
  // Global variables
long double E;			// exp(1). Useful for some test functions
long double pi;			// Useful for some test functions
struct matrix reflex1;
long double sqrtD;

  // File(s);
FILE * f_run;
FILE * f_synth;

  // =================================================
int main () 
{ 
	int d;			// Current dimension
	double error;			// Current error
	double errorMean;		// Average error
	double errorMin;		// Best result over all runs
	double errorMeanBest[R_max]; 
	double evalMean;		// Mean number of evaluations
	int functionCode;
	int i,j;
	int nFailure;		// Number of unsuccessful runs
	double logProgressMean;
	struct param param;
	struct problem pb; 
	int run, runMax; 
	struct result result; 
	double successRate;
	double variance;
		
	f_run = fopen ("f_run.txt", "w");  
	f_synth = fopen ("f_synth.txt", "w"); 
		
	E = exp ((long double) 1); 
	pi = acos ((long double) -1);
		
			
// ----------------------------------------------- PROBLEM
	functionCode = 0;	
			/* (see problemDef( ) for precise definitions)
				 0 Parabola (Sphere)
				 1 Griewank
				 2 Rosenbrock (Banana)
				 3 Rastrigin
				 4 Tripod (dimension 2)
				 5 Ackley
				100 Shifted Parabola 
				99 Test
			 */ 
		
	runMax = 30; // Numbers of runs
	if (runMax > R_max) runMax = R_max;
			
				
	// -----------------------------------------------------
	// PARAMETERS
	// * means "suggested value"		
			
	param.clamping =1;
			// 0 => no clamping AND no evaluation. WARNING: the program
			// 				may NEVER stop (in particular with option move 20 (jumps)) *
			// *1 => classical. Set to bounds, and velocity to zero
		
	param.randOrder=1; // 0 => at each iteration, particles are modified
						//     always according to the same order 0..S-1
						//*1 => at each iteration, particles numbers are
						//		randomly permutated
	param.rotation=0; // WARNING. Quite time consuming!
	// WARNING. Experimental code, completely valid only for dimension 2
	// 0 =>  sensitive to rotation of the system of coordinates
	//*1 => non sensitive (except side effects)
	// by using a rotated hypercube for the probability distribution
			
	param.stop = 0;	// Stop criterion
					// 0 => error < pb.epsilon
					// 1 => eval < pb.evalMax		
					// 2 => ||x-solution|| < pb.epsilon
					
		
	// -------------------------------------------------------
	// Some information
	printf ("\n Function %i ", functionCode);
	printf("\n (clamping, randOrder, rotation, stop_criterion) = (%i, %i, %i, %i)",
		param.clamping, param.randOrder, param.rotation, param.stop);
		
	// =========================================================== 
	// RUNs
			
	// Initialize some objects
	pb=problemDef(functionCode);
	
	// You may "manipulate" S, p, w and c
	// but here are the suggested values
	param.S = (int) (10 + 2 * sqrt(pb.SS.D));	// Swarm size
	if (param.S > S_max) param.S = S_max;
		
	printf("\n Swarm size %i", param.S);
	
	param.K=3; 											
	param.p=1-pow(1-(double)1/(param.S),param.K); 
	
	// According to Clerc's Stagnation Analysis
	param.w = 1 / (2 * log ((double) 2)); // 0.721
	param.c = 0.5 + log ((double) 2); // 1.193
		
	// According to Poli's Sampling Distribution of PSOs analysis	
	//param.w = ??; // in [0,1[
	//param.c =  
	//    smaller than 12*(param.w*param.w-1)/(5*param.w -7); 
	
	
	printf("\n c = %f,  w = %f",param.c, param.w);
	
	//---------------
	sqrtD=sqrt((long double) pb.SS.D);
		
	// Define just once the first reflexion matrix
	if(param.rotation>0)
	{					
		reflex1.size=pb.SS.D;	
		for (i=0;i<pb.SS.D;i++)
		{
			for (j=0;j<pb.SS.D;j++)
			{
				reflex1.v[i][j]=-2.0/pb.SS.D;
			}
		}
	
		for (d=0;d<pb.SS.D;d++)
		{
			reflex1.v[d][d]=1+reflex1.v[d][d];
		}
	}
								
	errorMean = 0;	    
	evalMean = 0;	    
	nFailure = 0;	
	//------------------------------------- RUNS	
	for (run = 0; run < runMax; run++)  
	{	
		//srand (clock () / 100);	// May improve pseudo-randomness            
		result = PSO (param, pb);
		error = result.error;
		
		if (error > pb.epsilon) // Failure
		{
			nFailure = nFailure + 1;
		}
	
		
			// Result display
		printf ("\nRun %i. Eval %f. Error %f ", run, result.nEval, error);
		printf("  x(0)= %f",result.SW.P[result.SW.best].x[0]);
			// Save result
			// fprintf( f_run, "\n%i %.0f %f ", run, result.nEval,  error );
			// for ( d = 0; d < SS.D; d++ ) fprintf( f_run, " %f",  bestResult.x[d] );
			
			// Compute/store some statistical information
		if (run == 0)
			errorMin = error;
		else if (error < errorMin)
			errorMin = error;
		evalMean = evalMean + result.nEval;	
		errorMean = errorMean + error;	
		errorMeanBest[run] = error;
		logProgressMean  = logProgressMean - log(error);		
	}		// End loop on "run"
			
	// ---------------------END 
	// Display some statistical information
	
	evalMean = evalMean / (double) runMax;   
	errorMean = errorMean / (double) runMax;
	logProgressMean = logProgressMean/(double) runMax;