pso8.c

C语言开发的微粒群优化算法源程序

/*

************** PSO 8.2: An Adaptive Particle Swarm Optimizer  ******************************

 [e-mail] 		2001-10-06      (ANSI C) 

What's new since?
(done) save data and post-mortem visualisation (cf Visu. Needs Leonardo free C compiler. See below)
(done) adaptive neighbourhood size for each particle 
(done) option "local queen"
(done) option "guided random generation"
(partly done) new estimation of the number of evaluations still needed 
(partly done) new metaphore similar to update_vel2, but for all kind of problems

Some future options
(in progress ) adaptive neighbourhood (each particule can choose its neighbours depending how much it trusts them)
(in progress) multicriteria optimization
(in progress) TSP

WARNING. There are probably still some bugs for energies (potential and kinetic)

This version does not _need_ any parameters at all, except, of course, the ones
defining the problem (file problems.txt):
- function
- dimension
- xmin,xmax, granularity (for the hyperparallepipedic search space)
- target
- admissible error
- granularity
- if all components of the solution has to be different (useful for some combinatorial problems)	

Nevertheless, for research purpose, you can modify the standard behaviour (see models.txt file).	

You can save some results (swarm positions, swarm energies) in text files, 
to import them later in a spreadsheet.
A 2 dimensional representation is also saved (file visu0.txt) for post-mortem animation 

(cf the program leonardo_visu.c)										

You can try to optimize some simple functions: 
1) find DD numbers that sum to S. 
2) find DD numbers whose product is S
3) Sphere. Find DD numbers  whose sum of squares is S
4) Rosenbrock (Banana)
5) Clerc's f1 (Alpine) (WARNING: new definition, different from the one in PSOx, x<=7)
6) Griewank
7) Rastrigin
8) Fifty-fifty Problem (in particular with granul=1)
9) Ackley
10) Foxholes  2D
11) Apple trees (see apple_trees.c)
12) Graph Coloring/Frequency Assignment (see color_kit.c)
13) PSO for PSO. Recursive use to find the best model to find a solution for a given problem (or set of problems)


15) El Farol


99) MyFunction (well, yours!) (see the source file myfunction.c)

(see also the file functions.txt).

Granularity. If you _know_  that the solution, on some dimensions,  is on an "integer" point, or even a rational one,
you can give this information
Fo example, if for a given dimension, the solution is an integer, then you can set granularity to 1 for this dimension:
the process may be faster (example: functions 4 and 6, or 1 and 2 if S is an integer)

Just for fun, using the function 2 (DD-product), you can then factorize an integer number.

IMPORTANT. This program is like a Swiss knife:
 it is not at all the best for any use, but you can use it as is for a lot of problems.

Values in the model.txt file are chosen so that the _first_ run shouldn't be too bad (engineer's point of view), 
even if you could have a better _mean_ result by trying a lot of times (mathematician point of view)

Have fun!

Equations
--------
For each particle, we have

v(t+1) = alpha*v(t) + beta*(pg(t) - x(t))
x(t+1) = x(t) + v(t+1)
with
v(t) := velocity at time t
x(t) := position at time t
pg(t):= previous best position found so far in the hood (at time t), including the particle itself

alpha is adaptive (theoreticalli. In practice, it is not necessary)
beta is depending on an adaptive b

Three options:

beta = b (the normal one)
beta = rand(b_min,b)  (for test purpose)
beta = rand(0, b) (for test purpose)

Note: unlike in the complete PSO version, 
the variable pi(t):= previous best position of the particle (at time t)
is NOT used


Three kinds of adaptations:

- adaptive swarm size
- adaptive coefficients alpha and beta
- adaptive neighbourhood size for each particle


Download
--------
Well, you have it, havn't you ?
For more information and last "public" version see
 http://www.mauriceclerc.net, Math stuff for PSO
And for even more information, see the Particle Swarm Central http://www.particleswarm.net

Compiling the program
-------------------
It is written in ANSI C.

If, by chance, you have a Mac, you should use the free compiler Leonardo
http://www.dis.uniroma1.it/~demetres/Leonardo/

  Note that a PC version might come.

Using the program
-----------------
See files problems.txt (and, if you want, models.txt)
"Problems" describe the problem(s) (function, dimension, search space etc.)
"Models" describe how PSO has to work (min max hood size,min max max swarm size etc.). 
You don't _need_ to modify it, but you can.
*/

#include <math.h>
#include <stdio.h>
#include <string.h>
#include <stdlib.h>
#include <time.h>

#define E				2.718281828459045
#define	max_chaos		1
#define min_chaos		0
#define pi				3.141592653589793238462
#define Max_DD 			100 // Max number of dimensions
#define	Max_histo		10 // Max number of previous f_values stored
#define Max_hood		100 // Max size of the neighbourhood

#define Max_swarm_size 	100 // Theoretically infinite, but my computer does not like infinite...
#define Max_vel			100 // Max number of velocity components (combinatorial kit). Must be >=Max_DD
#define n_funct			15 //  ********** Modify this number if you add some functions. See also functions.txt
#define my_function		99

struct search_space	{int dim;float min[Max_DD];float max[Max_DD];float dx[Max_DD];float volume;};
struct position		{float label;int size;float x[Max_DD];float f;float ep;float prev_f[Max_histo];int adm;float X1;float X2;};
//(size, coordinates,f_value= abs(objective function value-target), potential energy, previous f_values)

struct velocity		{int size;float v[Max_vel];float ek;int v2[Max_vel];int combin;}; // (Note: v2 is useful only for combinatorial kit)
//(size, components, kinetic energy)

struct particle		{struct position pos;struct velocity vel;struct position prev_best;int H;float prev_prev_best_f;float hood_size;};
struct subswarm		{int size;int rank[Max_swarm_size];int best;int worst;};
struct swarm		{int size;struct particle part[Max_swarm_size];int best;int worst;}; 	
struct param 		{int funct;int DD;struct search_space H;float target;float eps;int batch_runs;int printlevel;float Max_Eval;int save;int N;int times;int mine;int all_diff;char data[20];float used_constr;int init_file;char init_swarm_r[20];char init_swarm_w[20];float max_fr;float parapet;int combin;};
struct vector		{int size;float x[Max_DD];};
struct vector_int	{int size;int val[Max_DD];};
struct list			{int size;float val[Max_swarm_size];};
struct matrix		{int size;int val[Max_DD][Max_DD];int nb_cont;float max_cont;}; // For square matrix of constraints
struct model		{float alpha_max;float alpha_min;float b_max;float b_min;int hood_min;int hood_max;int init_size;int max_size;int min_size;int freedom;int init_type;int H;int parallel;int local_search;int rehope;int selection;int confin;int gener_type;int move_type;int queen;int guided_rand;int rand_type;};	
struct dplus		{struct vector_int dplus[Max_DD]; int dmax; int dtot;};
struct color		{int size;int val[Max_DD];int max;int adm;};
struct color_list	{int size;int val[Max_DD];};

// Subroutines (declarations)

float			adapt_alpha(float m,struct model model);
float 			adapt_b(float m,float alpha,struct model model);
float 			adapt_hood(int sw_size, float hood_size, float dhood, int type,int hood_min,int hood_max);
struct swarm	add_particle(struct swarm sw,struct param param, struct model model,int rank);
int 			alea(int min, int max);
float			alea_density(int density_step, int d,float xmin,float xmax,int type);
int 			alea_diff(int min,int max, int num);
float 			alea_float(float min, float max);
struct position	all_different(struct position pos,struct param param);
float			apple_trees(struct position pos,int print_level); // For "apple trees" example
int				best(struct swarm sw);
int				better_than(struct particle part1,struct particle part2,float eps);			
struct particle complete_part(struct particle par, struct param param, struct model model);
struct particle constrain(struct particle part,struct param param,struct model model);
void			display_density(int dim,int density_step);
void			display_matrix(struct matrix M); // display a constraint matrix (Frequency Assignment problem
void			display_model(struct model model, int level);
void			display_param(struct param param);
void			display_position(struct position pos);
void			display_swarm(struct swarm sw, int type);
void			display_velocity(struct velocity vel);
float			distance(struct position pos1,struct position pos2);
float			distance_min(struct swarm sw, struct position pos);
float 			diverg(struct position pos);
struct position domain_pos(struct position pos,struct param param);
float			energy_p(struct particle part,float eps); // Potentiel energy of a particle
float			energy_p_swarm(struct swarm sw); // Potential energy of the swarm
float			energy_k(struct particle part,struct search_space H); // Kinetic energy of a particle
float			energy_k_swarm(struct swarm sw); // Kinetic energy of the swarm
float 			estim_eval(struct swarm sw, struct param param,struct model model);
float 			estim_eval2(struct swarm sw, struct param param,struct model model,int cycle_size);
float 			estim_eval3(int nb_past_eval,float eps);
float			estim_eval_particle(struct particle part,int nb_past_eval,float eps);
float 			eval(float total,float target);
struct position	homogen_to_carte(struct position pos);
struct particle	hood_queen(struct swarm sw,struct subswarm hood_s,struct param param,struct model model);
struct subswarm	hood_swarm(struct swarm sw,int part_rank,float eps);
struct particle	init_particle(struct param param, int type, struct model model,struct swarm sw,int rank_model);
struct swarm 	init_swarm(struct param param,struct model model);
float			improvement(struct position pos);
float 			lipschitz(struct swarm sw, int part_rank, struct subswarm hood);
float 			MAX(float a,float b);
float 			max_comp(struct position pos);
float 			MIN(float a,float b);
float 			min_comp(struct position pos);
struct particle	move_particle(struct particle part,struct particle partg,struct param param, struct model model,int rank);
struct swarm	move_swarm(struct swarm sw,struct param param, struct model model,int rank);
float			MyFunction(struct position pos,struct param param, struct model model);
float			particle_improv(struct particle part);
struct model	pos_to_model(struct position pos, struct model model);
float			PSO(int recursive,struct model model); // Useful for PSO_for_PSO (function 13)
void 			print_N_run(struct position best_result,struct param param);
struct param	read_param(FILE *f_param, FILE *f_data,int begin);
struct model	read_model(FILE *f_models);
float			regranul(float x,float granul);
struct swarm	rehope(struct swarm sw,struct param param, struct model model);
struct swarm	remove_particle(struct swarm sw,int rank, int Min_size);
void			save_energy(FILE *f_energy,struct swarm sw);
void			save_result(FILE *f_run,FILE *f_init_w,struct position best_result,struct param param,int begin, struct  model model);
void			save_swarm(FILE *f_swarm,struct swarm sw);
void			save_visu(struct swarm sw,int level,float eps);
float 			tot(struct position pos,struct param param, struct model model);
float			tot_fifty_fifty(struct position pos);
void			update_density(struct position pos,struct search_space H,int step);
struct velocity	vel_granul(struct velocity vel,struct search_space H);
int				worst(struct swarm sw);

// Global variables

double			almostzero=0.0000001;  //to avoid overflow by dividing by too small value;
int				AS=-3; // Arbitrary value. Means "Assigned"
int				batch_runs; 
int				begin;
float			best_eval[Max_histo][2][2];
struct position	best_result[2]; // Two levels, for eventual recursive mode (PSO_for_PSO)
int				chance[2]; // Flag, for solution found by chance
float			chaos=(float)0.6; // Init value for chaotic function
clock_t			clock_tick;
int				density[10][Max_DD]={0}; // For Guided Random Generation. Cf init_particle()
int				density_step; // Number of points on each dimension to define "cells" for Guided random generation 
struct dplus	dplus; // To store a graph (for each node, list of neighbours). Just for Graph Coloring
float			duration;
float			duration_tot;

float			E0;
float			E_k; // Kinetic energy of the swarm
float			E_p; // Potential energy of the swarm
float			eval_f[2],eval_f_tot[2];
int				failure_tot[2];
char			functions[30][30]; // Function names
int				H; // Flag for Coloring problem. Indicates what kind of projection to do
float			label[2]; // To labellize particles
int				level; // For recursive mode
struct matrix	M;
int				Max_hood_size; // Just for info
int				max_rand;
int				Max_swarm; // Just for info
float			Mean_swarm; // Just for info
struct model	models[2]; // Model level 0 and Initial model level 1 for recursive mode (PSO_for_PSO)
struct model	model2; // Modified model for recursive mode (PSO_for_PSO)
float			move[2]; // For statistical success estimations
int				n_add;
int				n_rand; // Type of pseudorandomness (see alea_float())
int				n_remove;
int				NA=-32000;// Arbitrary value. Means "non assigned"
int				NO=-2; // Arbitrary value. Means "no change"
struct param	param2; // For recursive mode (PSO_for_PSO)
float			rand_list[1000];
int				rand_type; // Type of randomness. Cf alea_float()
int				recursive=0;
int				self_better[2]; // Just for stat
struct	swarm	sw;
float			two_pi;
float			volume_checked[2];

// For Leonardo visualization

float			X1,X2;
FILE *f_visu0; // For visualization (X1,X2,value) (level 0)
FILE *f_visu1;

//--------
FILE *f_init_r; // (optional) to Read an initial swarm
FILE *f_init_w; // to Write an initial swarm for future use

FILE *f_data; // For additional data, defining the problem (see Coloring/Frequency Assignment)
FILE *f_energy; // To save swarm's energy
FILE *f_functs; // Function names
FILE *f_models; // PSO models (level 0, and 1 if recursive)
FILE *f_move_color; // For Coloring problem. Save particle positions in a 2D space (colors, constraints)
FILE *f_param; 	// Parameter file
FILE *f_rand;   // List of random numbers in [0,1]
FILE *f_run; 	// To save the run
FILE *f_swarm;  // To save swarm
FILE *f_synth; // 
FILE *f_trace;

// External subroutines

#include "color_kit.h" // For Graph coloring problems
#include "combinatorial_kit.h" // For combinatorial problem
#include "color_kit.c"
#include "combinatorial_kit.c"
#include "myfunction.c" // Write your own objective function in this file.
#include "apple_trees.c" // For "apple trees" example
#include "hyper_vol_1_300.txt" // List of unit hypersphere volumes for dimensions 1,2,...300

/* =============================================================================== */
void main()
{
int				i;
float			Mean_eval;
char 			mr[1];

f_functs=fopen("functions.txt","r");
f_models=fopen("models.txt","r");

two_pi=(float)(2*pi);

/* Initialize pseudorandom numbers  (see alea_float() ), just in case
   you use this option (see models.txt, model.rand_type=1)
*/ 

f_rand=fopen("rand_list.txt","r");
fscanf(f_rand,"%i",&max_rand);

for (i=0;i<max_rand;i++)
{
	fscanf(f_rand,"%f",&rand_list[i]);
}
fclose (f_rand);

// Initialize function names. Just for display

for (i=1;i<=n_funct;i++)
	{
	fscanf(f_functs,"%s\n",&functions[i]);
	printf("\n %i %s",i,functions[i]);
	}
strcpy(functions[99],"My function ");
printf("\n-----------------------------------\n");

fclose (f_functs);

// Model(s)
printf("\n Read Model 0");
models[0]=read_model(f_models); // Level 0 model
printf("\n Read Model 1 (in case you use recursive PSO)");
models[1]=read_model(f_models); // Level 1 model (useful only for recursive use  (see param.funct=13 in tot() )
fclose (f_models);
display_model(models[0],0);

//===================  Other files  (they must be open during the whole process)
f_visu0=fopen("visu0.txt","w"); // For post mortem visualisation
f_visu1=fopen("visu1.txt","w"); // The same, in case of recursive PSO

f_energy=fopen("energy.txt","w");
f_move_color=fopen("f_move_color.txt","w");
f_param=fopen("problems.txt","r");

f_run=fopen("run.txt","w");
f_swarm=fopen("swarm.txt","w");
f_synth=fopen("synth.txt","w");
f_trace=fopen("trace.txt","w");
//-----------------


level=0;
best_result[0].size=0; // Still no "best result" at all


Mean_eval=PSO(0,models[0]); // ******* CORE OF THE PROGRAM *******


printf("\n Any letter and Return: ");
scanf("%s",&mr); 
return;
}

/*============================ S U B R O U T I N E S =====================================*/
/* Except the first one, they are in alphabetical order                                   */

float PSO(int recursive,struct model model)  // Main routine. Separated for optional recursive use
{
/* Main subroutine */

int				add[Max_swarm_size];
int				best;
struct color	color;
float 			currenteps;
int				cycle,cycle_size;
float			dhood;
float			enough_improv;
float			eval_f_max;
float			eval_f_prev;
float			eval_f_sup=0; // Just for particular tests
float			eval_run[100]; // Just for information
struct subswarm	hood_s;
int				hood_size;
int				i;
int				i_step;
int				j;
float			improv;
int 			Iter;
float			mean_eval;

float			min_eps;
float			nb_eval_max;
float			nb_eval_min;
float			nb_no_progress;
int				nb_rehope;
struct param 	param={0};
float			previous_eps;
struct particle	queen;
int				rank;
int				remove[Max_swarm_size];
int				size0;
float			stop;
clock_t 		ticks0,ticks1, ticks2;
int				times;
int				worst;
float			zzz;

if (param.printlevel>0)
	printf("\n Level %i. PSO", level);

rand_type=model.rand_type; // Global variable. Type of randomness (see alea_float())	
n_rand=0; // Number of times alea_float() is called
stop=40000; /* Arbitrary stop criterion (objective function evaluations)
              You may remove it, for there is also an adaptive automatic stop criterion
 			*/
zzz=0;
nb_eval_max=0;
nb_eval_min=(float)-NA; 

begin=0;

if (level==0)
	{
	save_result(f_run,f_init_w, best_result[level],param,0, model); // Titles
	fprintf(f_energy,"eval_f swarm_size potential_energy kinetic_energy\n");
	}

// Parameters of the problem(s)
read_param: