main.c

这是Clerc最新的文章balanced PSO,包括论文和与其配套的源程序

/*
Balanced PSO
Research version
(based on Standard PSO 2007)
 Contact for remarks, suggestions etc.:
 [e-mail]

 Last update
2009-01-12 FEmaxSize parameter. Now, one can have far more quickly the curve
           "success probability" vs "number of fitness evaluations"
2009-01-09 Pressure vessel, compression spring
2009-01-05 CEC 2005 F2, F7, F9 (not rotated)
2008-12-27 Feasible space (FS), included in the search space (SS). Confinement inside the 
           search space, and no evaluation for positions in SS but not in FS
2008-12-22 Two strategies (only different w, for the moment). STILL TESTING
2008-12-14 Biased random initialisation
2008-12-01 Min potential initialisation
2008-11-24 Hammersley initialisation
2008-10-31 Gaussian local search
2008-10-13 Three CEC 2005 functions (F1, F6, F9)

Explanation
The "local area"  (i.e. "around the best known positions"
is rigorously defined.
With such a definition it appears that,
contrarily to usually said, Standard PSO does not perform any
exploitation at all, as soon as dimension is greater than 6.

(FYI the "exploitation rate" of each iteration is saved into the file
f_exploit.txt)

So, in this variant, a cheap local search is added.
I also replaced the (bad) pseudo random number generator of C by KISS


----------------------------------- Parameters/Options
clamping := true/false => whether to use clamping positions or not
init := different possible initialisations
randOrder:= true/false => whether to avoid the bias due to the loop
				on particles "for s = 1 to swarm_size ..." 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(c1)*(p(t)-x(t)) + R(c2)*(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))    
 */  

/* Codification
Each option has a "code". So, for example
PSO C1 P0 V5 L0 R0 r0 O1 S0
is for this PSO whith the following options:  
param.clamping =1;
param.init=0;
param.initVel=5;
param.localSearch=0;
param.R=0;
param.rand=0;
param.randOrder=1;
param.strat=0; 

i.e. Standard PSO 2007

By convention, we call it just SPSO07, and we add a code only when some options
are modified. For example in
PSO P2 V1
we use option 2 for the initialisation of the positions,
and option 1 for the initialisation of the velocities

*/

/* TO TRY
	- NoHope/ReHope 

*/
  
#include "stdio.h"
#include "math.h"
#include <stdlib.h>
#include <time.h>

#define ulong unsigned long // To generate pseudo-random numbers with KISS
#define RAND_MAX_KISS ((unsigned long) 4294967295)
  
#define	D_max 160		//160.  Max number of dimensions of the search space
#define R_max 5000	// 1000 Max number of runs
#define	S_max 58		//21 Max swarm size
#define stratMax 3 // Maximum number of different strategies
#define zero  1.0e-30			// 1.0e-30 // To avoid numerical instabilities with some compilers
#define infinity 1.0e30
#define NMax 1000		// Max number of points when the landscape is read on a file
#define ClassMax 20 // Max number of points for the curve "success rate vs FEmax"

#define ERROR(a) {printf("ERROR: %s\n", a);exit(-1);}
  // Structures
struct quantum 
{
  double q[D_max];
  int size; 
};

struct SS 
{ 
	int D;
	double max[D_max];  // Initialisation, and feasible space
	double min[D_max]; 
	struct quantum q;		// Quantisation step size. 0 => continuous problem
	double maxS[D_max]; // Search space
	double minS[D_max];
};

struct param 
{
  	double c1;		// Confidence coefficient
		double c2;
  	int clamping;	// Position clamping or not
	int K;			// Max number of particles informed by a given one
	double exploit1;
	double exploit2;
	int explOption;
	int init;
	int initVel;
	int localSearch; //
	double p;		// Probability threshold for random topology	
					// (is actually computed as p(S,K) )
	int rand; // Kind of randomness
	int randOrder;	// Random choice of particles or not
  	int S;			// Swarm size
  	int stop;		// Flag for stop criterion
  	double w[stratMax];		// Confidence coefficient
		double w0Rate; // Rate of particles using local strategy 0 (when param.strat=1)
  	int R; 			// What kind of random distribution
		double sigma; // Standard deviation if Gaussian distribution (R=1)
	int strat; // Strategy to use
	
		
};

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

struct velocity 
{  
	int size;  
	long double v[D_max]; 
};

struct problem 
{ 
	double epsilon; 	// Admissible error
	int evalMax; 		// Maximum number of fitness evaluations
	int function; 		// Function code
	int nb;// Number of points for the "repartition on nb points" problem (code -1)
	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
	struct position Pp[S_max]; // Previous previous best (for information)
	int lBest[S_max]; // Flag: lBest[s]=1 means Pp[s] was a local best
	int S; 						// Swarm size 
	struct velocity V[S_max];	// Velocities
  struct position X[S_max];	// Positions
	double w[S_max];
	double c[S_max];
	int R[S_max];

};

struct result 
{  
	double nEval; 		// Number of evaluations  
	struct swarm SW;	// Final swarm
	double error;		// Numerical result of the run
};

struct lBest {int rank; double x[D_max];}; // To sort the local bests

// Local exploitation areas
// rank[] = flags. rank[i]=1 means "particle i is a local best"
// ex[] = search space "around" each local best
struct exploit {int rank[S_max]; struct SS exI[S_max];};
struct posList {int size; struct position p[D_max+1];}; 

struct landscape {int N;double x[NMax];double fx[NMax];}; // 1D-landscape
  // Sub-programs
ulong	rand_kiss(); // For the pseudo-random number generator KISS
void	seed_rand_kiss(ulong seed); // "
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);
//double normL (struct velocity v,double L);
struct position hammersley(struct SS SS, int S, int deb, int s);
double perf (struct position x, struct problem pb);	// Fitness evaluation
struct position quantis (struct position x, struct SS SS);
struct problem problemDef(int functionCode, FILE *data);
struct result PSO ( struct param param, struct problem problem, int level);
int sign (double x);

// Subroutines for local search
struct exploit exploitArea(struct swarm SW, struct SS SS,double rho);
double explNb(struct swarm SW, struct exploit ex );
static int compareLBest (void const *a, void const *b); // For qsort
//struct result localSearch(struct result R,struct problem pb, 
//					struct param param,int g[],struct exploit ex,double rho);
  // Global variables
long double E;			// exp(1). Useful for some test functions
long double pi;			// Useful for some test functions
struct landscape funct; // When (x,f(x)) is read on a file
int dLBest; // For information
int randOption; // Kind of randomness to use
double distPow; // For decreasing random distribution coefficient  (R>=2)

int successFlag[R_max][ClassMax]; // For ClassMax different FEmax
int FEmax[ClassMax];
int FEmaxSize;
int FEclass;
int run;

  // File(s);
FILE * f_run;
FILE * f_synth;
FILE * f_expl;
FILE * fLandscape;
FILE * f_init;
FILE * f_init_save;
FILE * fCurve;

  // =================================================
int main () 
{ 
	struct position bestBest;
	int d;			// Current dimension
	int dim;
	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 n;
	int nFailure;		// Number of unsuccessful runs
	double logProgressMean;
	struct param param;
	struct problem pb; 
	int  runMax; 
	struct result result; 
	double sum1,sum2;
	int SPSO07;
	double successRate;
	double variance;
	int wSize;
	double zz;
		
	f_run = fopen ("f_run.txt", "w");  
	f_synth = fopen ("f_synth.txt", "w"); 
	f_expl = fopen ("f_expl.txt", "w"); // Just for info (exploitation rate)
	fLandscape=fopen("fLandscape.txt","r"); // When data are on a file
	f_init_save=fopen("f_init_save.txt","w");
	fCurve=fopen("fCurve.txt","w");

	E = exp ((long double) 1); 
	pi = acos ((long double) -1);


FEmaxSize=20; // If >0, save info in order to plot the curve (success_rate vs FEmax)
for(n=0;n<FEmaxSize;n++) FEmax[n]=4000+1000*(n+1); 
	
param.stop = 1;	// Stop criterion
					// 0 => error < pb.epsilon  
					// 1 => eval = pb.evalMax	

SPSO07=0; // 1 => "simulate" Standard PSO 2007. 
					// All other parameters/options	are ignored
			
// ----------------------------------------------- PROBLEM
	functionCode = 8;	
			/* (see problemDef( ) for precise definitions)
				 0 Parabola (Sphere)
				 1 Griewank
				 2 Rosenbrock (Banana)
				 3 Rastrigin
				 2 Rosenbrock (Banana)
				 3 Rastrigin
				 4 Tripod (dimension 2)
				 5 Ackley
		  	 6 Center-bias test function
				 7 Pressure vessel
				 8 Compression spring

				CEC 2005 benchmark
				100 F1 (shifted Parabola/Sphere) 30D
				102 F6 (shifted Rosenbrock) 10D
				103 F9 (shifted Rastrigin) 10D 
				104 F2 Schwefel
				105 F7 Griewank (NOT rotated)
				106 F9 Ackley (NOT rotated)

				99 Test
				1000 Read (x,f(x)) on the file data.txt

				-1  Initialisation of N particles in a D-space => D*N variables
				-2  idem, but with another potential function
			 */ 
		
	runMax = 500; // Numbers of runs
	if (runMax > R_max) runMax = R_max;
			
				
	// -----------------------------------------------------
	// PARAMETERS
	// * = "suggested value"	

	if(SPSO07==1) goto runs;

			
	param.clamping =1; // code C (for "clamping")
			// 0 => no clamping AND no evaluation. WARNING: the program
			// 				may NEVER stop 
			// *1 => classical. Set to bounds, and velocity to zero (Standard PSO 2007)
//------
		param.init=2; // P.  Initialisation of the positions
		// -1 => read from file f_init (be sure it is consistent with the other data)
		// 0 => on each dimensionrandom position (uniform U(xmin,xmax) (Standard PSO 2007)
		// 1 => random position (centre biased)
		//* 2 => (improved) Hammersley
		// 3 => random (uniform) or Hammersley
		// 4 => on "discrete" regular points on each dimension
		// 5 =>  random position (centre biased), 2

		param.initVel=1; // V. Initialisation of the velocity for particle whose position is X
		// 0 => on each dimension U(xmin,xmax) -U(xmin,xmax)
		// * 1 => on each dimension U(xmin,xmax) - X(d)
		// 2 =>  Y - X, where Y is the initial position of another particle, chosen at random
		// 3 => set to zero
		// 4 => either 1 or 2
		// 5 => on each dimension 0.5*( U(xmin,xmax) -U(xmin,xmax)) 
                                  // ("Half-diff", i.e. Standard PSO 2007)
		
		// NOTE: Don't forget  to set param.strat (see below)

	param.localSearch=0; // L
			// 0 => no local search // SPSO 07
			// 1	=> random in some local areas
			// 2 => uniform in the best local area  
			// 3 => middle of some local areas
			// 4 => middle of the best local area *?
			// 5 => Gaussian in the best local area
		

//	param.exploit1; // Define the rate exploit/explor you want
												// Useful only if localSearch=1 or 2
//	param.exploit2;// Define the initial size of the exploitation area
												// Useful only if localSearch>0
//	param.explOption
	// 		0 => constant relative size of the exploitation area
	// 		1 => variable
	// 		2 => variable
	// *	3 => alternative {exploit2,1/S}
	//		4 => random uniform [0,param.exploit2]
	//		5 => random Gauss (1/S, exploit2)

// Suggested values

	switch (param.localSearch)
	{
		default:
			param.exploit1=0.5;
			param.exploit2=1./3;
		break;
		case 0:
		param.exploit2=1./2; // Just for information
												// and to check the exploitation rate
		param.explOption=0;
		break;
		case 1:
			param.exploit1=0.5;
			param.exploit2=1./3;
		break;
		case 2: // *
			param.exploit1=0.5;
			param.exploit2=1./2;
			param.explOption=3;
		break;
		case 3:
			param.exploit1=0.5;
			param.explOption=0;
		break;
		case 4:
			param.exploit1=0.5;
			param.explOption=0;
		break;
		case 5:
			param.exploit1=0.5;
			param.exploit2=1./2;
			param.explOption=3;
		break;
	}


//------	
	
  param.R=3; // R
  		// 0 => uniform distribution in [0,1]  // SPSO 07
  		// 1 => Gaussian distribution	
			//* 2 => decreasing function of the 1D distance (see also global variable distPow)
			//* 3 => either 0 or 2
	distPow=2;

	param.sigma=0.1; // Standard deviation for Gaussian distribution
										// Useful only if R>0

	param.rand=0; // r
								//* 0 => KISS pseudo-random number generator
								// 1 => the one of your compiler
  		
	param.randOrder=1; // O
						// 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  (SPSO07)

	param.strat=0; // S
			// General strategy to use  
			// 0 => classical (just one constant w) (SPSO 07)
			// 1 => constant w0 or random w in [w0,w1], according to w0Rate  (see  below)
			// 2 => idem, but non uniform w
			// 3 => random different w for each particle
			// 4 => forced "scout" particle(s)
			// 5 => variable c

	switch(param.strat)
	{
			default: case 3: wSize=1; break;
			case 1: case 2: case 4: wSize=2; break;
	}

	
runs: // =========================================================== 

	// RUNs
	randOption=param.rand;			
	// Initialize some objects

	pb=problemDef(functionCode,fLandscape);
	printf("\n Dimension %i",pb.SS.D);

	if(SPSO07==1)
	{
		#include "parameters.c"
		goto begin;
	}

	// 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;
	
	param.K=3; 											
	param.p=1-pow(1-(double)1/(param.S),param.K);
	 											
	// According to Clerc's Stagnation Analysis
	param.w[0] = 1 / (2 * log ((double) 2)); // 0.721
	param.c2 = 0.5 + log ((double) 2); // 1.193
	param.c1=param.c2; 
	param.w[1]=1;
//param.w[1]=0.9;
	param.w0Rate=1-1./param.S; //   In order to have about one "scout particle" 
	//param.w0Rate=1-1./(param.p*param.S);
					 // In order to have (in average) one "scout particle" in each neighbourhood

begin:
	// Some information
	printf ("\n Function %i ", functionCode);
	printf("\n Swarm size %i", param.S);
	printf("\n Mean size of the neighbourhood %f", param.p*param.S);

	printf("\n (clamping, R, randOrder, exploitation, stop_criterion) = (%i,%i, %i, %i, %i)",
		param.clamping, param.R, param.randOrder,param.explOption,  param.stop);
	printf("\n randomness %i",param.rand);
	printf("\n init %i",param.init);
	printf("\n initVel %i",param.initVel);

	printf("\n Strategy= %i",param.strat);
	printf("\n w ");
	for(d=0;d<wSize;d++)
			printf(" %f",param.w[d]);