pso8.c

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

float diverg(struct position pos)
{
// Compute how much the direction of the position is far from the "balanced" one (1,1, ...1)
// (for 2D visualization)


int		d;
float	r;
float	rx;
float 	x;

r=0;
for (d=0;d<pos.size;d++)
{
	r=r+pos.x[d]*pos.x[d];
}	
r=(float)sqrt (r);

rx=(float)pos.size;
rx=(float)(r/sqrt (rx));
x=0;
for (d=0;d<pos.size;d++)
{
	x=x+(pos.x[d]-rx)*(pos.x[d]-rx);
}	
x=(float)sqrt(x);

return x;

}



// ----------------------------------------------------------------------------- DOMAIN_POS
struct position domain_pos(struct position pos,struct param param)
{
// Modify the position according to the "granularity" of the search space
int				d;
struct position post;

post=pos;
for (d=0;d<pos.size;d++)
{
	post.x[d]=regranul(pos.x[d],param.H.dx[d]);
}
	
return post;	
}		

// ----------------------------------------------------------------------------- ENERGY_K
float energy_k(struct particle part,struct search_space H) // Kinetic energy  ( Sum of V^2/2 )
{ // Kinetic energy of a particle
int		d;
float	dx;
float 	e;
float	v;

e=0;

for (d=0;d<part.vel.size;d++)
{
	v=part.vel.v[d];
	dx=H.max[d]-H.min[d];
	if (dx>almostzero)
		e=e+v*v/(dx*dx);
}

return e/2;
}

// ----------------------------------------------------------------------------- ENERGY_K_SWARM
float energy_k_swarm(struct swarm sw) // Kinetic energy  ( Sum of V^2/2 )
{ // Kinetic energy of a swarm
int		i;
float 	e;

e=0;

for (i=0;i<sw.size;i++)
{
	e=e+sw.part[i].vel.ek;
}
return e;
}

// ----------------------------------------------------------------------------- ENERGY_P
float energy_p(struct particle part, float eps)
{ // Potential energy of a particle
float e;
if (eps>almostzero) {e=part.pos.f/eps;} else {e=part.pos.f;};

return e;

}
// ----------------------------------------------------------------------------- ENERGY_P_SWARM
float energy_p_swarm(struct swarm sw)  // Potential energy (Sum of abs(f(X)-target) )
{ // Potential energy of the swarm
int		i;
float 	e;

e=0;

for (i=0;i<sw.size;i++)
{
	e=e+sw.part[i].pos.ep;
}
return e;
}

// ----------------------------------------------------------------------------- ESTIM_EVAL
float estim_eval(struct swarm sw, struct param param,struct model model)
{/* Re-estimate how many evaluations are probably still needed to reach a solution
	globally, for the whole swarm
*/


float	estim;
int		hood;
int		i;
float	nb_estim;

if (param.combin==1) return (float)(1/almostzero); // Combinatorial problem

estim=0;
nb_estim=0;
hood=model.hood_min;

for (i=0;i<sw.size;i++) // Estimate for each particle, and take the max
{
	if(sw.part[i].pos.prev_f[0]<sw.part[i].pos.prev_f[hood-1])
	{
		estim=MAX(estim,estim_eval_particle(sw.part[i],MIN(hood,Max_histo),param.eps));
		nb_estim=nb_estim+1;
	}
}

if (nb_estim==0)
{
	if (param.printlevel>1)
		printf("\n WARNING: impossible to estimate the number of evaluations needed");
	estim=(float)(1/almostzero);
}
else
{
	estim=sw.size*(estim+hood);// Add cycle_size and multiply by the number of particles

	estim=10*estim; // Arbitrary pessimistic point of view
	estim=MIN(estim,(float)(1+1/almostzero));
}

if (model.move_type==2) // Each move needs two evaluations
{
	estim=2*estim;
}

if (model.local_search==1)
{
	estim=estim+param.DD;
}

if (model.local_search>1)
{
	estim=estim+param.DD*param.DD;
}

if (param.funct==12) // Coloring problem is more difficult ...
{
	estim=estim*dplus.dmax;
}
return estim; 
}

// ----------------------------------------------------------------------------- ESTIM_EVAL2
float estim_eval2(struct swarm sw, struct param param,struct model model, int cycle_size)
{/* Re-estimate how many evaluations are probably still needed to reach a solution
	globally, for the whole swarm

Note: hyper_vol[] is a static table, containing unit hypersphere volumes for dimension 1,2,...300
Note: level and volume_checked[level] are global variables
*/
int				d;
float			df;
float			df_mean;
float			dim;
float			dmax;
float			ds;
float			dx;
float			dx_mean;
float			estim;
int				i;
int				j;
int				k;
struct subswarm	hood;
int				ncount;
float			r;
float			search_power_hood;
float			search_power_swarm;
float			swarm_volume;

if (param.combin==1) return (float)(1/almostzero); // Combinatorial problem. The following heuristic does'nt work

dim=(float) param.H.dim;
search_power_swarm=0;

/*
swarm_volume=1;
for (d=0;d<dim;d++) // Volume of the "best" swarm
{
	dmax=0;
	for (i=0;i<sw.size;i++) // For each particle ...
	{
		for (j=0;j<i;j++)
		{
			dx=fabs(sw.part[i].prev_best.x[d]-sw.part[j].prev_best.x[d]);
			if (dx>dmax) dmax=dx;
		}
	}
	swarm_volume=swarm_volume*dx;
}
*/

for (i=0;i<sw.size;i++) // For each particle ...
{
	hood=hood_swarm(sw,i,param.eps); // ... find the neighbourhood

	dx_mean=0;
	df_mean=0;
	ncount=0;
	for (j=0;j<hood.size;j++) // For each neighbour (including the particle itself ...
	{
		k=hood.rank[j];
		dx=distance (sw.part[i].pos,sw.part[k].prev_best);
		if (dx<almostzero) continue;
			df=fabs (sw.part[i].pos.f-sw.part[k].prev_best.f);
			dx_mean=dx_mean+dx;
			df_mean=df_mean+df/dx;	
			ncount=ncount+1;
	}	
	df=df_mean/(float)ncount;
	
	r=param.eps/df;
		
	if (dim>1)
	{
		ds=hyper_vol[(int)(dim-2)]*pow(r,dim-1);// small hypersphere
	}
	else
	{
		ds= 1;			
	}
			
	// "volume" of the cylinder which could be checked
	// (Note: this formula assumes that alpha is constant)														
	search_power_hood=(model.b_max-model.b_min)*ds*dx_mean;  
	
	search_power_swarm=search_power_swarm+search_power_hood;
}

if (search_power_swarm<almostzero) return 1/almostzero;
printf("\n volume_checked[level] %f, search_power_swarm %f, H.volume %f",volume_checked[level],search_power_swarm,param.H.volume);

volume_checked[level]=volume_checked[level]+search_power_swarm*(float)cycle_size;
estim=((param.H.volume-volume_checked[level])/search_power_swarm)*(float)sw.size; // Estimation of the number of evaluations

return estim; 
}

// ----------------------------------------------------------------------------- ESTIM_EVAL3
float estim_eval3(int nb_past_eval,float eps)
{
 /* Re-estimate how many evaluations are probably still to do
    We consider the points (u,v) := (nb_of_eval, best_f_value)
	and approximate the curve with v=lambda/u (option: lambda/u^2)  Then, we solve v<eps. 
Note: best_eval[][2] is a global variable
Note: level is a global variable
 */

int		i;
double	delta;
double	lambda;
double	mu;
int		nb;
int		option;
double 	s1,s2;
double	u;
double	u0;
double	u1;
double	uc;
double	v1;
double	vc;
option=1;

// model lambda/(u-u0)^mu
nb=nb_past_eval;
if (nb>Max_histo-1) nb=Max_histo-1;
s1=0;s2=0;
u0=best_eval[nb][0][level];
u1=best_eval[nb-1][0][level];
v1=best_eval[nb-1][1][level];
uc=best_eval[0][0][level];
vc=best_eval[0][1][level];

if (v1<almostzero) return 1/almostzero;
if (v1<=vc) return 1/almostzero;
if (u1<=u0) return 1/almostzero;

mu=log(vc/v1)/log((u1-u0)/(uc-u0));
lambda=pow(u1-u0,mu)*v1;

//printf("\n log(vc/v1) %f",log(vc/v1));
//printf(" log(...) %f",log((u1-u0)/(uc-u0)));
//printf("\n v1,vc: %f %f.  lambda, mu: %f %f",v1,vc,lambda,mu);

u=pow(2*lambda/eps,1/mu) + u0-uc;
return u;

/*
for (i=0;i<nb;i++)
{
	u=best_eval[i][0][level]-u0;
//printf("\n u,v: %f %f",u, best_eval[i][1][level]);
	if (u<almostzero) return 1/almostzero;
	
	if (option==2)
	{
		s1=s1+best_eval[i][1][level]/(u*u);
		s2=s2+1/(u*u*u*u);
	}
	else
	{
		s1=s1+best_eval[i][1][level]/(u);
		s2=s2+1/(u*u);
	}
	
}
//if (s2<almostzero) return 1/almostzero;

lambda=s1/s2; 

//printf("\n nb %i, s1, s2, %f %f, lambda %f, u0, u %f %f",nb,s1,s2,lambda,u0,uc);
//printf("\n");
//for (i=0;i<=nb;i++) printf(" %.0f",best_eval[i][0][level]);

if (option==2)
{
	delta=vc/(lambda*(uc-u0)*(uc-u0));
	 u=sqrt(lambda*delta/eps);
	} 
else 
{
	delta=vc/(lambda*(uc-u0));
	u=lambda*delta/eps;
}

u=u+u0-uc;

return u;
*/

}


// ----------------------------------------------------------------------------- ESTIM_EVAL_PARTICLE
float estim_eval_particle(struct particle part,int nb_past_eval, float eps)
{
 /* Re-estimate how many evaluations are probably still to do
    using previous evaluations for a particle which have improved its position
	For this particle, we consider the points (u,v) := (evaluation rank, f value)
	and approximate the curve with v=lambda/u. Then, we solve v<eps. 
 */


int		i;
float 	estim;
float 	s1,s2;
float	x;

s1=0;s2=0;
for (i=0;i<nb_past_eval;i++)
{
	x=(float)(nb_past_eval-i);
	s1=s1+part.pos.prev_f[i]/x;
	s2=s2+1/(x*x);
}
estim=s1/s2; 

return estim/eps;
}

// ----------------------------------------------------------------------------- EVAL
float eval(float total,float target)
{ // How far from the target is the value "total"
float v;
v=(float)fabs(target-total);
return v;
}

// ----------------------------------------------------------------------------- HOMOGEN_TO_CARTE
struct position	homogen_to_carte(struct position pos)
{
/* Convert a position given by its homogeneous coordinates into 
the same position given by its cartesian coordinates. The homogeneous coordinates
are supposed to be given relatively to the unit points:
(0, ... 0), (1,0,...0) (0,1,0...,0) ... (0,...,0,1)

See, for example, Apple Trees problem

*/
int				d;
struct position post;
float			s;

//post=pos;
post.size=pos.size-1;

s=0;
for (d=0;d<pos.size;d++)
{
		s=s+pos.x[d];
	}
if (fabs (s)<almostzero)
{
	printf("\n WARNING. Incorrect homogeneous coordinates. Return null point");
	for (d=0;d<post.size;d++)
		post.x[d]=0;
	return post;
}

for (d=0;d<post.size;d++)
	post.x[d]=pos.x[d+1]/s;

return post;
}

// ----------------------------------------------------------------------------- HOOD_QUEEN
struct particle	hood_queen(struct swarm sw,struct subswarm hood_s,struct param param,struct model model)
{
/*
 Weighted position (queen) of the particles in a given neighbourhood,
 using a second order approximation of the objective function
*/
float			c;
float			c_tot;
int				d;
int				dim;
int				i;
struct particle queen={0};
int				rank;

dim=sw.part[0].pos.size;

c_tot=0; 
for (i=0;i<hood_s.size;i++) // Total weight
{