pso8.c

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

if (recursive<2)
{
	param=read_param(f_param,f_data,begin); // READ THE PROBLEM DESCRIPTION

	if (param.funct==13) display_model(models[1],1); // If PSO for PSO, display the second model
	
	// Check improvement every cycle_size time steps

	cycle_size=model.min_size; // Initial value
	
	// Just to avoid too long computation
	param.parapet=(float)(model.hood_min*model.max_size*model.hood_min);
	
	// Combinatorial problems are more difficult
	if (param.combin>0 || param.funct==12) param.parapet=param.parapet*(model.max_size-1)/2; 
	

		
	if (param.printlevel>0)
	{
		printf("\n=======================================\n Maximum swarm size: %i",model.max_size);
		printf("\n In any case, no more than %.0f evaluations without improvement",param.parapet);
	}

	param2=param; // Save for future use in recursive mode
	model2=model;
	display_param(param);	
}
else
{
	if (level>0)
	{
		param=param2;
	
		if (best_result[level-1].size>0)
		{
			model=pos_to_model(best_result[level-1],model2); // The new model at this second level is the best "position"
												// at the previous level
		}
		else
		{
			model=model2;
		}
	}											
}

if (begin==0) {batch_runs=param.batch_runs;begin=1;}

times=0; // you can run several times the same problem (parameter param.times)
ticks0=clock();
duration_tot=0;
volume_checked[level]=0; // Just for information and estimation of time still needed

min_eps=1/almostzero;

//========================================================= LOOP on problems to solve
times_loop: 
label[level]=0; // Label of the first particle. Just for later visualisation
nb_rehope=0; // Number of times Rehope process has been launched

failure_tot[level]=0; // Number of failures


if(param.printlevel>0)
{
	printf("\n\n--------------\n Level %i. %s %iD",level,functions[param.funct],param.DD);
	printf("  Run %i/%i",times+1,param.times);
}

ticks1=clock(); // Just for information. To evaluate processor time needed

chance[level]=0; // Equal to 1 if a solution is found
eval_f_tot[level] =0; // Number of objective function evaluations
eval_f_prev=0;
eval_f_max=(float)(1/almostzero); // At the very beginning, choose a very high max number of evaluations
eval_f[level]=0; // Number of evaluations counter
nb_no_progress=0;// If this number becomes > parapet, then => Failure
Iter=1;		// Number of iterations. Just for information
self_better[0]=0; // Just for stats
self_better[1]=0;
density_step=model.guided_rand;
dhood=1;// Increment for adaptive hood size
param.N=model.init_size;// Normally we begin with just a few particles, but you can change it (see models.txt)
H=1; // Just for Coloring problem. Projection will be done first on H1 (minimize constraints)
rank=-1;
if (model.freedom==6) rank=0; // Special case "indepnedent dimensions"

sw=init_swarm(param, model); 	// ******** SWARM INITIALIZATION *********

if (param.save==2) save_visu(sw,level,param.eps);
	    
param.N=sw.size; // Useful if the initial swarm has been in fact read on a file, at least partly 
Max_swarm=sw.size; // Just for information
Mean_swarm=(float)sw.size; // Just for information
Max_hood_size=0;// Just for information


if (chance[level]==1) goto end_by_chance; // Sometimes, you find a solution by chance at the very beginning
				
best_result[level]=sw.part[sw.best].prev_best; // Keep trace of the best result found so far
previous_eps=best_result[level].f;
currenteps=previous_eps;

for (i=0;i<Max_histo;i++) {best_eval[i][0][level]=eval_f_tot[level];best_eval[i][1][level]=currenteps;}

// Initial value for "enough improvement", to begin
enough_improv=(sw.part[sw.worst].prev_best.f-sw.part[sw.best].prev_best.f)/sw.part[sw.worst].prev_best.f;

if (param.printlevel>1) 
{
	printf("\n  Best result after %.0f evaluations ",eval_f_tot[level]);
	display_position (best_result[level]);
	printf("\n------ MOVING the swarm (size %i)",sw.size);
}						

cycle=0;  // Number of iterations since last adaptive update


do_Iter: //termination criterion is (currenteps<=eps or eval_f>=eval_f_max or eval_f_tot>stop)
						
	Iter=Iter+1; // number of time steps. Just for information
	move[0]=0;move[1]=0; // Number of moves, number of succesful moves
	if (param.save==1) 
	{
		save_swarm(f_swarm,sw);	
		save_energy(f_energy,sw);
	}

	eval_f_prev=eval_f_tot[level];
	sw=move_swarm(sw,param, model,rank); 					// *** HERE ! THE SWARM MOVES ! **********
	
	if (param.save==2) save_visu(sw,level,param.eps);
	if (param.save==3) fprintf(f_trace,"\n %.0f %.0f",move[0],move[1]);

	if (best_result[level].f<min_eps) min_eps=best_result[level].f;
	if(chance[level]==1) goto end_by_chance; // A solution has been found, maybe before all particles have moved
	
	if (rank>=0)
	{
		if (fabs (best_result[level].x[rank]-sw.part[sw.best].prev_best.x[rank])<param.eps)
		{
			rank=rank+1; if (rank==param.DD) rank=0; // Special case "independent dimensions"
		}
	}

	// Eventually update the best result
	if (sw.part[sw.best].prev_best.f<best_result[level].f)
	{
		best_result[level]=sw.part[sw.best].prev_best;// Abs(target-tot) for the best particle
		currenteps=best_result[level].f;	
	}
	
		
	if (param.printlevel>2) display_swarm(sw,1); 
	if (param.printlevel>1)
	{
		printf("\n  Best result after %.0f evaluations ",eval_f_tot[level]); 
		display_position(best_result[level]);
	}

	// For recursive PSO, systematically save intermediate results	
	if (param.funct==13)	 
	{	
		save_result(f_run,f_init_w,best_result[level],param,1, model);
	}

	if (param.printlevel>2)
	{
		printf("\n Energies (potential, kinetic):  %f  %f",E_p,E_k);
	}

	if (eval_f_tot[level]>stop) // Arbitrary stop criterion
	{
		failure_tot[level]=failure_tot[level]+1;
		printf("\n FAILURE eval_tot  %.0f > max= %.0f",eval_f_tot[level],stop);
		goto the_end;
	}

	if (eval_f[level]>eval_f_max) // Adaptive stop criterion
	{
		failure_tot[level]=failure_tot[level]+1;
			printf("\n FAILURE eval  %.0f > max=%.0f",eval_f[level],eval_f_max); 
		goto the_end;
	}

	                                // ***********  ADAPTATION(S)  (begin) *************
	cycle=cycle+1;
	cycle_size=sw.size;

	if (2*cycle<cycle_size) goto end_adapt; // Adaptation is not made at each iteration
										 // just from time to time (after "cycle_size" iterations) 
										// Note that cycle_size is not constant
									
		
		/* If improvement, "bad" particles are not useful => remove them
		   If deterioration, it is worthwhile to generate a particle
		*/
		
		// Store best evaluation  (just for info and prediction purpose		
		for (i=0;i<Max_histo-1;i++)
		{
			best_eval[Max_histo-i-1][0][level]=best_eval[Max_histo-i-2][0][level];
			best_eval[Max_histo-i-1][1][level]=best_eval[Max_histo-i-2][1][level];
		}
		best_eval[0][0][level]=eval_f_tot[level];
		best_eval[0][1][level]=currenteps;

		// Some inits for adaptation
		n_add=0; // Number of particles added. Just for information
		n_remove=0; // Number of particles removed. Just for information
		size0=sw.size;
		
		if (model.hood_min==0) // Option (cf models.txt): Neighbour = whole swarm
		{
			hood_size=sw.size;
		}
	
		for (i=0;i<size0;i++) //  Reinit add[] and remove[], status of particles to add or remove
		{
			add[i]=0;
			remove[i]=0;
		}
		
/*
		for (i=0;i<size0;i++) // Save info about evolution of neighbourhood sizes
		{
			fprintf(f_run,"\n %.0f %f %f", sw.part[i].pos.label, eval_f_tot[level],sw.part[i].hood_size);
		}
*/
			
		i=0;
		next_hood: // Loop on neighbourhoods, to check them 
		{
			hood_s=hood_swarm(sw,i,param.eps); // Generate the neighbourhood
			best=hood_s.rank[hood_s.best];
			worst=hood_s.rank[hood_s.worst];
			
			// Evaluate local improvement
			
			improv=particle_improv(sw.part[best]);

			if (improv<enough_improv) // If not enough improvement
			{ 
				if(model.queen==1 && sw.size<Max_swarm_size)// Try the local queen
				{
					queen=hood_queen(sw,hood_s,param,model);
			
					if (queen.pos.f<=param.eps)
					{
						best_result[level]=queen.pos;
						chance[level]=1;
						printf("\n Solution thanks to a queen");
						goto end_by_chance;
					}

					if( queen.pos.f<sw.part[worst].prev_best.f) 
					{
						//printf("\n queen is better");								
						sw.part[sw.size]=queen;  // add the queen
						sw.size=sw.size+1;
						n_add=n_add+1; 
						goto next_part;
					}
				}

				add[best]=1;
				// Adapt hood size (increased because not enough improvement)
				sw.part[best].hood_size=adapt_hood(sw.size,sw.part[best].hood_size,dhood,1,model.hood_min,model.hood_max);
				goto next_part;
			}

			if (improv>=enough_improv) // If  enough improvement, remove the worst particle
			{
				remove[worst]=1;
				// Adapt hood size (decreased, because enough improvement)
				sw.part[best].hood_size=adapt_hood(sw.size,sw.part[best].hood_size,dhood,0,model.hood_min,model.hood_max);
				goto next_part;
			}


next_part:
			// Just once by complete neighbourhood
			if (model.hood_min>0) hood_size=(int)sw.part[i].hood_size; 
			i_step=(int)(hood_size+(float)hood_size/2 +1); 
			//i_step=1; // For each particle. For test purpose
			i=i+i_step;
			if (i<size0) goto next_hood;
		}


	size0=sw.size;
	
	// WARNING. Do add before to remove. Rank i is used in gener_particle_color
		for (i=0;i<size0;i++) // Adding
		{
			if (add[i]==1)
			{
				sw=add_particle(sw,param, model,i);  // add a particle 
				if(chance[level]==1) goto end_by_chance;

				// enough_improv is probably too "ambitious" => decreasing
				//enough_improv=enough_improv/(1+1/(float)sw.size);
				//enough_improv=enough_improv/(2-1/((float)sw.size)+1);	
				enough_improv=enough_improv/(2-1/exp((float)sw.size));
				
						
			}
		}
			
		for (i=0;i<size0;i++) // Removing
		{
			j=size0-i-1;
			if (remove[j]==1)
			{
				sw=remove_particle(sw,j, model.min_size);

				// enough_improv is reasonable. One try to increase it
				//enough_improv=enough_improv*(1+1/(float)sw.size);
				//enough_improv=enough_improv*(2-1/((float)sw.size)+1);
				enough_improv=enough_improv*(2-1/exp((float)sw.size));
					
			}
		}

		param.N=sw.size; // For eventual re-initialisation
		Max_swarm=(int)MAX((float)Max_swarm,(float)sw.size); // Just for information
		
		cycle=0;

		eval_f_max=estim_eval(sw,param,model); // **** HERE THE NUMBER OF EVALUATIONS TO DO IS RE-EVALUATED ***
//printf("\n estim1 %.0f",eval_f_max);		
		//eval_f_max=estim_eval2(sw,param,model,cycle_size);
//printf("\n  estim2 %.0f",eval_f_max);
	//eval_f_max=estim_eval3(3,param.eps);
//printf("\n  estim3 %.0f",eval_f_max);
		eval_f[level]=0;
	 
			

		if (param.printlevel>0)
		{
			printf("\n");
			printf(" Eval: %.0f",eval_f_tot[level]);
		 	printf("  Estim.: %.0f",eval_f_max);
			printf(" Result: %.8f",currenteps);
			printf(" (part.  %.0f).",sw.part[sw.best].pos.label);
			printf(" Swarm size %i",sw.size);
			printf(" (+ %i, -%i)",n_add,n_remove);
		}
		
		
		end_adapt:;
		
		                     // ---------------------------  ADAPTATION(S)  (end)	

		Mean_swarm=Mean_swarm+sw.size;
		
		if (currenteps>=previous_eps)
		{
			if (currenteps>previous_eps) // Normally, it can't be the case, but there may be a bug...
				printf("\n WARNING. Non monotonous convergence. eps= %f (was %f)",currenteps,previous_eps);	
			nb_no_progress=nb_no_progress+1; // Number of iterations without any improvement

			if (model.rehope>0 && nb_no_progress>param.parapet) 
			{
				sw=rehope(sw,param, model); // REHOPE
				
				if(chance[level]==1) goto end_by_chance; // A solution has been found, maybe before all particles have moved 
				
					best_result[level]=sw.part[sw.best].prev_best;// Abs(target-tot) for the best particle
					currenteps=best_result[level].f;
					
				nb_rehope=nb_rehope+1;
				printf("%i/%i",nb_rehope, model.rehope);

				if (nb_rehope>model.rehope)
				{
					failure_tot[level]=failure_tot[level]+1;
					goto the_end;
				}

				eval_f_max=(float)(1/almostzero);  
				nb_no_progress=0;
			}
		}
		
		previous_eps=currenteps;			
		goto do_Iter; // Iterates till a criterion is met 
			
end_by_chance:

currenteps=best_result[level].f;
if(param.printlevel>0)
	printf("\n CHANCE !");

the_end: // ---------------------------- END

ticks2=clock();
clock_tick=ticks2-ticks1; // Processor time needed

if (param.funct==12) // For Coloring/ Frequency assignment, simplify best_result to display/save
	best_result[level]=color_shift(best_result[level]);

if (level==0)
{						
	save_result(f_run,f_init_w,best_result[level],param,1, model);
}

if(param.save==1) save_swarm(f_swarm,sw);
	
	duration = (float)(clock_tick/CLOCKS_PER_SEC);
	duration_tot=duration_tot+duration;
	eval_run[times]=eval_f_tot[level];
	
	if (level==0)
	{
		printf ("\n Run %i => eps %f, eval: %.0f, duration: %f",times+1,currenteps,eval_f_tot[level],duration);
		printf(" (by part. %.0f)",best_result[level].label);
		printf("\n Max swarm size = %i.",Max_swarm);
	

		if (times>0) printf(" Mean swarm size =  %.0f.",Mean_swarm/Iter);
		printf(" Max neighbourhood size = %i",Max_hood_size);
		
		if (model.move_type==2)  // Stats about  this option (self move), if used
		{
			printf(" Self better  %.2f ",(float)self_better[1]/(float)self_better[0]);
		}

		if (param.funct==12)  // Special case Coloring problem
		{
		// Better to verify ...
			color=convert_position_to_color(best_result[level]);
			i=check_color(color,M,dplus);
			printf(" Unsatisfied constraint(s): %i",i);
			i=diff_color(best_result[level]);
			printf(", %i colors",i);
		} 
	}

zzz=zzz+eval_f_tot[level];

if (eval_f_tot[level]>nb_eval_max) nb_eval_max=eval_f_tot[level];
if (eval_f_tot[level]<nb_eval_min) nb_eval_min=eval_f_tot[level];

if(param.printlevel>0)			
	print_N_run(best_result[level],param);
	
if (param.printlevel>1) display_swarm(sw,1);
			
	times=times+1;
	if (times<param.times) goto times_loop;
	
		ticks2=clock();
		mean_eval=zzz/param.times;
		Mean_swarm=Mean_swarm/Iter; // Just for information, mean swarm size