main.c

用C编写的粒子群寻优的优化算法

 {
    //--------------------------------- If still no improvement, rebuild  tours
    if (move4==1) printf("  /rebuild, random");
    if (move4==2) printf("  /rebuild, arc_around");
    if (move4==3) printf("  /rebuild, blackboard");

    for (i=0;i<sw1.size;i++)
	{
         switch( move4)
         {
          case 1:
          norm_v=alea(1,G.N);
          sw1.p[i].v=alea_velocity(norm_v);    // random velocity
          sw1.p[i]=pos_plus_vel(sw1.p[i],sw1.p[i].v,0,0,0); //  => random new position
           break;

           case 2:
              new_p=arc_around(sw1,sw1.p[i]);
               sw1.p[i].v=coeff_pos_minus_pos(new_p,sw1.p[i],monotony,0,equiv_v,trace);
               sw1.p[i]=new_p;
           break;

           case 3:   
              new_p=BB_use(sw1.p[i]);
               sw1.p[i].v=coeff_pos_minus_pos(new_p,sw1.p[i],monotony,0,equiv_v,trace);
               sw1.p[i]=new_p;
               break;
         }
			if (sw1.p[i].best_f<best_best.best_f)/* Best of the best after the move */
			{
				best_best=sw1.p[i];
                     if (best_best.best_f<=target)  goto success;// Success
                    printf("\neval %.0f  value %f",eval_f, best_best.best_f);
                    goto moves;
			}
         	if (eval_f>=max_eval)  goto end_max_eval;
   }
}

// Primitive rehope
printf("\n Primitive rehope");
 for (i=0;i<sw1.size;i++)
 {
       sw1.p[i].v=alea_velocity(G.N);
       sw1.p[i]= pos_plus_vel(sw1.p[i],sw1.p[i].v,0,0,0);
 }

 goto moves;


//-----------------------------------
 end_move:

	if (save!=0) save_swarm(f_run,sw1);
      if (trace>1) display_swarm(sw1);
      printf("\n Eval %.0f",eval_f);
	if (trace>0) display_position(best_best,2);

      goto moves;

end_max_eval:
		printf("\n\n Stop at max eval %.0f",eval_f);
         goto end;

success:
n_success=n_success+1;
printf("\n SUCCESS");
            
end:
// Give the best position (sequence) found 

printf("\n Best solution found after %.0f evaluations",eval_f);
 if (best_best.best_x.s[0]!=0) best_best.best_x=rotate(best_best.best_x,0);
display_position(best_best,2);

 // Check if no error
 zzz=f(best_best,-1,-1);
 if (zzz!=best_best.best_f)
 {
     printf("\n ERROR. The true f value is in fact %f ",zzz);
 }
zzz=best_best.best_f;
if (zzz<eps_min) eps_min= zzz;
 if (zzz>eps_max) eps_max= zzz;
 eps_mean=eps_mean+zzz;
 mean_eval=mean_eval+eval_f;

if (save!=0) save_swarm(f_run,sw1);
save_swarm(f_trace,sw1);

if (trace>1) display_swarm(sw1);

run=run+1;
if (run<=run_max) goto loop_run;
else
{

printf("\n With  swarm size=%i, max eval=%.0f => %i successes/%i. Rate=%.3f ",swarm_size,max_eval,n_success,run_max,(float)n_success/run_max);

printf("\n Mean  %f [%f, %f]",eps_mean/run_max,eps_min,eps_max);
mean_eval=mean_eval/run_max;
printf("\n Mean eval %f",mean_eval);

goto the_end;
}

error_end:
printf("\n Sorry, I give up");

the_end:
  return EXIT_SUCCESS;
}
 /* ============================= Subroutines ================== */

/*------------------------------------------------------------------- 	CONVERGENCE_CASE */
// Useful only for "classical" non adaptive PSO. Will be removed when merged in TRIBES

struct coeff convergence_case(float kappa, float phi)
//*****************************General case
// Coefficients of the generalized iterative representation
// The matrix of the system is M
//       [alpha       beta*phi   ]
//       [-gamma  (delta-eta*phi)]       with phi=phi1+phi2
// NOTE: useful only for "classical" (non spherical) methods (see move_towards())
//********************************************
{
float			Aclass2,Bclass2;
struct coeff	coefft;
float			disc;
int				i;
float			khi;
float			maxne,maxneprim;
float			maxneclass2;
float			ne1class2,ne2class2;


casechoice:
printf("\nWhich convergence case ? (0,1,2 3,4,5) (suggestion: 2): ");
scanf("%i",&conv_case);

if (conv_case<0 || conv_case>5)
	{
	printf("\n Wrong value. Le'ts try again");
	goto casechoice;
	}

coefft.c[5]=phi;

switch(conv_case)
	{
	case 0:  // So that x(t+1) = x(t) + v(t+1)
		maxne=max_eig_basic(phi);
		khi=kappa/maxne;
		coefft.c[0] = khi;
		coefft.c[1] = khi;
		coefft.c[2] = khi;
		coefft.c[3] = (khi+1)*phi-1;
		coefft.c[4] = 1;
		break;



	case 1:  //(proved case)
		maxne=max_eig_basic(phi);
		khi=kappa/maxne; //********* to change according to the syggestion given
                     	//by the program
                     	// khi=kappa/maxne for sure convergence
		coefft.c[0] = khi;
		coefft.c[1] = khi;
		coefft.c[2] = khi;
		coefft.c[3] = khi;
		coefft.c[4] = khi;
		break;

	case 2: //(just for phi>=4) //(proved case)

		if (phi>=4)
			{
			disc=sqrt(phi*phi-4*phi);
			Aclass2=3/2-(3/4)*phi+(3/4)*disc-(1/2)*(phi*phi)+(1/4)*(phi*phi*phi)-(3/4)*(phi*phi)*disc;
			Bclass2=(1/4)*fabs(2-phi-2*(phi*phi)+phi*phi*phi+disc-3*phi*phi*disc);
			ne1class2=fabs(Aclass2+Bclass2);
			ne2class2=fabs(Aclass2-Bclass2);
			maxneclass2=ne2class2;
			if(ne1class2>ne2class2)
				maxneclass2=ne1class2;

			khi=kappa/maxneclass2;     //khi=kappa/maxneclass2 for sure convergence


			coefft.c[3]=khi*(2-phi+disc)/2;
			coefft.c[2]=khi*(2-phi+3*disc)/(4*phi);
			coefft.c[0]=2*coefft.c[3];
			coefft.c[1]=coefft.c[0];
			coefft.c[4]=2*coefft.c[2];
			}
		else
			{
			printf("\n For this case, you must have phi>=4. But here you have phi = %f",phi);
			goto casechoice;
			}
		break;

	case 3: //(experimental case)
		maxne=max_eig_basic(phi);
		khi=kappa/maxne;   //to begin, but not sure convergence
		printf("\n Initial khi: %f",khi);
  		 withkhi:
		if(khi>almostzero)
			{
			coefft.c[0]=1;
			coefft.c[1]=1;
			coefft.c[2]=khi;
			coefft.c[3]=khi*khi;
			coefft.c[4]=khi*khi;
			maxneprim=max_eig(coefft,phi);
			if (maxneprim>1)
				{
				khi=0.9*khi;
				goto withkhi;
				}
			}
		else
			{
			printf("\n No possible convergence");
			printf("\n You should change some parameters (phi1, phi2)");
			goto casechoice;
			}
		printf("\n    Final khi (which is used): %f",khi);
		break;

	case 4: //Called type 1' in the theoretical paper
		maxne=max_eig_basic(phi);
		khi=kappa/maxne;
		coefft.c[0] = khi;
		coefft.c[1] = khi;
		coefft.c[2] = 1;
		coefft.c[3] = 1;
		coefft.c[4] =1;
		break;

	case 5:  // So that x(t+1) = p + v(t+1)
		maxne=max_eig_basic(phi);
		khi=kappa/maxne;
		coefft.c[0] = khi;
		coefft.c[1] = khi;
		coefft.c[2] = khi;
		coefft.c[3] = (1-khi)*phi;
		//coefft.c[3] = (1-khi)*phi/2;  When phi is random ?
		coefft.c[4] = 1;
		break;

	default:
	break;
	}


	printf("\n Coeffs:");
	for (i=0;i<5;i++)
		printf(" %f", coefft.c[i]);

return coefft;
}
//============================================================== BB_UPDATE
void BB_update(struct particle p)
{
 /*
   Update the blackboard. Must be used each time a particle modify its  position
   BB is a global variable.
   G is a global variable.
   BB_option is a global variable
   ------------------------
   BB_option=0:
   For each arc (n1,n2) of x, the value f is added to BB[n1][n2][0]
   and 1 is added to BB[n1][n2][1].
   BB[n1][n2][0]/BB[n1][n2][1] is seen as the probable "cost" (penalty) of
   using the arc (n1,n2) in a tour.
   ----------------------  
   BB_option=1:
     For each arc (n1,n2) of x, the value f replace the previous one if it is smaller,
       and 1 is added to BB[n1][n2][1].
 */
 if (move[4]<3) return; // The blackboard won't be used. No need to update it!
 if (trace>1)
    printf("\n BB_update");
 int  i;
 int n1,n2;
 
  for (i=0;i<G.N;i++) // Warning. This suppose the last node is equal to the first one
  {
     n1=p.x.s[i];
     n2=p.x.s[i+1];
    if (BB_option==0)
    {
       BB[0].v[n1][n2]=BB[0].v[n1][n2]+p.f;
       BB[1].v[n1][n2]=BB[1].v[n1][n2]+1;
    }
    else
    {
       if ( BB[1].v[n1][n2]>0 && p.f<BB[0].v[n1][n2])
       {
          BB[1].v[n1][n2]=BB[1].v[n1][n2]+1;
          BB[0].v[n1][n2]=p.f;
       }
    }
  }

}

//============================================================== BB_USE
struct particle BB_use(struct particle p)
{  /*
 Build a tour, using the blackboard
  For each arc (n,*) to add, make a probabilistic choice,
  according to the line n of BB
 */

double   big_value;
int   candidate[Nmax];
double   candidate_penalty[Nmax];
int         candidate_nb;
double   c_max;
int      i;
int      j;
int      k;
int      m;
int      node1,node2;
double   penalty;
 struct particle pt;
 double  r;
 int  rank;
 int     used[Nmax]={0};

big_value=G.N*G.l_max;

 pt=p;
  node1=0; // The new tour begins in 0
  rank=1;
used[node1]=1;
pt.x.s[node1]=0;

//printf("\n");
 next_node:
candidate_nb=0;
  for (node2=0;node2<G.N;node2++) // For each possible new node ...
  {
     if (used[node2]>0) continue;
     if(rank==G.N-1) goto set_node; // Last node => no choice, anyway

     // Ask the blackboard
          if (BB[1].v[node1][node2]==0)  penalty=big_value;
          else
          {
             if (BB_option==0)
                penalty= BB[0].v[node1][node2]/BB[1].v[node1][node2] ;
             else
                penalty= BB[0].v[node1][node2];
          }
 
      candidate_penalty[candidate_nb]=penalty;
      candidate[candidate_nb]=node2;
      candidate_nb=candidate_nb+1;
   }   // next node2

   if (candidate_nb==1) goto set_node;  // Fi there is just one possible node, add it
   
   // Choose the new node at random, according to penalties
   c_max=0;
   for (i=0;i<candidate_nb;i++) c_max=c_max+ candidate_penalty[i];
   for (i=0;i<candidate_nb;i++) candidate_penalty[i]=1-candidate_penalty[i]/c_max;//Proba
   for (i=1;i<candidate_nb;i++) candidate_penalty[i]=candidate_penalty[i-1]+candidate_penalty[i];  // Cumulate
   r=alea_float(0,candidate_nb-1); // random choice

    //printf("\n %i: ",candidate_nb);for (i=0;i<candidate_nb;i++) printf(" %f",  candidate_penalty[i]);

   for (i=0;i<candidate_nb; i++)
   {
      if(r<candidate_penalty[i])
      {
         node2=candidate[i];
         goto set_node;
        }
    }
   set_node:
    //printf("\n r %f, i %i node2 %i",r,i,node2);
//printf(" %i",node2);
   pt.x.s[rank]=node2;  // Add the node
   used[node2]=1;   //  remember it is already used
   rank=rank+1;
   if (rank<G.N) goto next_node;

   pt.x.s[G.N]=pt.x.s[0]; // Complete the tour
   pt.f=f(pt,-1,-1); // evaluate
   if (pt.f<pt.best_f)  // eventually update the best previous position
   {
          pt.best_x=pt.x;
          pt.best_f=pt.f;
          BB_update(pt); // Update the blackboard
   }
  if (trace>1)
  {
     display_position(pt,1);
     display_position(pt,2);
  }
  BB_update(pt); // Update the blackboard
 return pt;
 }

//============================================================== ARC_AROUND
struct particle arc_around(struct swarm sw, struct particle p) 
{ // Build a tour, asking opinion of other particles for each arc
double   big_value;
int   candidate[Nmax];
double   candidate_penalty[Nmax];
int         candidate_nb;
double   c_max;
int      i;
int      j;
int      k;
int      m;
int      node1,node2;
double   penalty;
 struct particle pt;
 double  r;
 int  rank;
 int     used[Nmax]={0};

big_value=G.N*G.l_max;
 
 pt=p;
  node1=0; // The new tour begins in 0
  rank=1;
used[node1]=1;

 next_node:
candidate_nb=0;
  for (node2=0;node2<G.N;node2++) // For each possible new node ...
  {
     if (used[node2]>0) continue;
     if(rank==G.N-1) goto set_node; // Last node => no choice, anyway
     
       penalty=0; 
      for (j=0;j<sw.size;j++) // ... ask opinions
      {
         for (k=0;k<G.N-1;k++) // Find the arc in the current particle, if exists
          if (sw.p[k].best_x.s[k]==node1 && sw.p[k].best_x.s[k+1]==node2 )
          {
            penalty=sw.p[k].best_f;
             goto next_opinions;
          }
         // Can't find it in the current particle
         penalty=big_value;

      next_opinions:;
      }  
      // Mean value
      candidate_penalty[candidate_nb]=penalty/sw.size;
      candidate[candidate_nb]=node2;
      candidate_nb=candidate_nb+1;
   }   // next node2

   // Choose the new node at random, according to penalties
   c_max=0;
   for (i=0;i<candidate_nb;i++) c_max=c_max+ candidate_penalty[i];
   for (i=0;i<candidate_nb;i++) candidate_penalty[i]=1-candidate_penalty[i]/c_max;//Proba
   for (i=1;i<candidate_nb;i++) candidate_penalty[i]=candidate_penalty[i-1]+candidate_penalty[i];  // Cumulate
   r=alea_float(0,candidate_nb-1); // random choice
  
    //printf("\n %i: ",candidate_nb);for (i=0;i<candidate_nb;i++) printf(" %f",  candidate_penalty[i]);
 
   for (i=0;i<candidate_nb; i++)
   {
      if(r<candidate_penalty[i])
      {
         node2=candidate[i];
         goto set_node;
        }
    }
   set_node:
    // printf("\n r %f, i %i node2 %i",r,i,node2);
   pt.x.s[rank]=node2;  // Add the node
   used[node2]=1;   //  remember it is already used
   rank=rank+1;
   if (rank<G.N) goto next_node;
   
   pt.x.s[G.N]=pt.x.s[0]; // Complete the tour
   pt.f=f(pt,-1,-1); // evaluate
   if (pt.f<pt.best_f)  // eventually update the best previous position
   {
          pt.best_x=pt.x;
          pt.best_f=pt.f;
           pt.best_time=time;
   }
  if (trace>1)
  {
     display_position(pt,1);
     display_position(pt,2);
  }
  BB_update(pt); // Update the blackboard  
 return pt;
}