psolib.h

pso算法编程的c语言实现

/***************************************************************************
                          psolib.h  -  description
                             -------------------
    begin                : Sun Jan 28 2001
    copyright            : (C) 2001 by Max Salazar
    email                : [e-mail]
 ***************************************************************************/

/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/

//Copia el arreglo s a d (De tipo double y longitud l)
void copy_array(double *s, double *d, unsigned int l) {
  unsigned int i;

  for(i = 0; i < l; i++)
    d[i] = s[i];
	
}

//Busca el valor maximo de un arreglo de tipo double
unsigned int search_max(double *f, unsigned int M) {
  unsigned int i, fmax = 0;

  for(i = 1; i < M; i++)
    if (f[fmax] < f[i])
      fmax = i;
	
  return fmax;
}

//Busca el valor maximo de un arreglo de tipo entero sin signo
unsigned int search_max(unsigned int *f, unsigned int M) {
  unsigned int i, fmax = 0;

  for(i = 1; i < M; i++)
    if (f[fmax] < f[i])
      fmax = i;
	
  return fmax;
}

//Busca el valor minimo de un arreglo de tipo double
unsigned int search_min(double *f, unsigned int M) {
  unsigned int i, fmin = 0;

  for(i = 1; i < M; i++)
    if (f[fmin] > f[i])
      fmin = i;
	
  return fmin;
}

//Busca el valor minimo de un arreglo de tipo entero sin signo
unsigned int search_min(unsigned int *f, unsigned int M) {
  unsigned int i, fmin = 0;

  for(i = 1; i < M; i++)
    if (f[fmin] > f[i])
      fmin = i;
	
  return fmin;
}
//*
void muta(double *pop, int M, unsigned int D,int Gen,int totGen,int fun){
  unsigned int i,j;
  int dimension=0;
  double linf[D],lsup[D],linftemp,lsuptemp,rango;
  double valtemp=0;
  ranges(fun,D,linf,lsup);
  for(i = 0;i < M;i++){
    if(flip(pow(1.0-(double)Gen/(totGen*pM),1.5))){
      dimension=RandomInt(0,D-1);
      rango=(lsup[dimension]-linf[dimension])*pow(1.0-(double)Gen/(totGen*pM),1.5)/2;//totGen
      valtemp=RandomDouble(rango,-rango);
      if(pop[dimension+(D*i)]-rango<linf[dimension]) linftemp=linf[dimension]; 
      else linftemp=pop[dimension+(D*i)]-rango;
      if(pop[dimension+(D*i)]+rango>lsup[dimension]) lsuptemp=lsup[dimension]; 
      else lsuptemp=pop[dimension+(D*i)]+rango;

      pop[dimension+(D*i)]=RandomDouble(linftemp,lsuptemp);//mutation_operator(pop[ngen+(D*i)],Gen);

    }
  }
}

//Se evaluan las particulas
void evaluation(double *pop, double *fitness, unsigned int M, unsigned int D, unsigned int fun, unsigned int NF) {
  unsigned int i, j, f;
		
  for(i = 0; i < M; i++)
    for(j = 0, f = fun; j < NF; j++, f++)
      fitness[(NF * i) + j] = chfun(&pop[D * i], f);
}

//Calcula la velocidad de las particulas
double velocity(double W, double Vi, double C1, double C2, double Pb, double Pg, double Xi) {
  return W * Vi + C1 * RandomDouble(0.0, 1.0) * (Pb - Xi)
    + C2 * RandomDouble(0.0, 1.0) * (Pg - Xi);
}

//Funcion que regresa 1 en caso de que el indiviudo (pop) a
//evaluar sea mejor que el establecido (pbests), de lo contrario
//la funcion devuelve 0.
unsigned int g_is_best(double *pop, double *pbests, unsigned int fun, unsigned int rest, double fitness, double fbests, unsigned int opt) {
  unsigned int i, j;
  double cont1, cont2, aux;
	
  i = 0;
  cont1 = 0;
  cont2 = 0;
 		
  for(j = 0; j < rest; j++)
    aux = chconst(pop, fun, j + 1);
  if (aux > 0.0)				//<<<<----Cuando las restricciones deben ser menores o iguales a 0 (>)
    cont1 += aux;
 		
  for(j = 0; j < rest; j++)
    aux = chconst(pbests, fun, j + 1);
  if (aux > 0.0)				//<<<<----Cuando las restricciones deben ser menores o iguales a 0 (>)
    cont2 += aux;
 					
  if ((cont1 == 0) && (cont2 == 0)) {
    if (((!(opt)) && (fitness < fbests)) || ((opt) && (fitness > fbests)))	// <<<------ min(<) - max(>)
      i = 1;
  }
  else {
    if ((cont1 == 0) || (cont2 == 0)) {
      if (cont1 == 0)
	i = 1;
    }
    else {
      if (cont1 < cont2)
	i = 1;
      else
	if ((cont1 == cont2) && (((!(opt)) && (fitness < fbests)) || ((opt) && (fitness > fbests))))	// <<<------ min(<) - max(>)
	  i = 1;
    }
  }
	
  return i;
}

//Intercambia los registros de los arreglos en
//caso de que se trate de una mejor particula.
void ifbest_interchange(double *fitness, double *pop, double *fbests, double *pbests, unsigned int M, unsigned int D, unsigned int fun, unsigned int rest, unsigned int opt) {
  unsigned int i, j;
	
  for(i = 0; i < M; i++) {
    if (g_is_best(&pop[D * i], &pbests[D * i], fun, rest, fitness[i], fbests[i], opt)) {
      fbests[i] = fitness[i];
      for(j = 0; j < D; j++)
	pbests[(D * i) + j] = pop[(D * i) + j];
    }
  }
	
	
}

//Busca a la mejor de las particulas, intentando
//que cumpla con todas sus restricciones
unsigned int search_opt_cons(double *f, unsigned int M, unsigned int D, unsigned int rest, unsigned int fun, double *pop, unsigned int opt) {
  unsigned int i, j, band = 0, cmin = 0, cmax = 0, it = 0, bestp;
  double aux, *cons;
  cons = new double [M];
	
  for(i = 0; i < M; i++) {
    cons[i] = 0;
    for(j = 0; j < rest; j++) {
      aux = chconst(&pop[D * i], fun, j + 1);
      if (aux > 0.0)				//<<<<----Cuando las restricciones deben ser menores o iguales a 0 (>)
	cons[i] += aux;
    }
  }
	
  cmax = search_max(cons, M);
  cmin = search_min(cons, M);
	
  bestp = cmin;
	
  if(cons[bestp] == 0.0)
    band = 1;
		
  do {
    if (band)
      if (cons[cmin] == 0.0) {
	if (((!(opt)) && (f[cmin] < f[bestp])) || ((opt) && (f[cmin] > f[bestp]))) {	// <<<------ min(<) - max(>)
	  bestp = cmin;
	}
	cons[cmin] = cons[cmax] + 1.0;
	cmin = search_min(cons, M);
      }
    band = 0;
				
    it++;
  } while ((band) && (it < M));
	
  delete [] cons;

  return bestp;	
}


///////////////////////////////////////////////////////
//MULTI-OBJETIVO
///////////////////////////////////////////////////////
/*
  unsigned int sumcompVec(double *noDomF, unsigned int NF, unsigned int h, unsigned int i, unsigned int opt) {
  unsigned int sum = 0, j;
	
  for(j = 0; j < NF; j++)
  if (((noDomF[(NF * h) + j] < noDomF[(NF * i) + j]) && (opt == 0)) || ((noDomF[(NF * h) + j] > noDomF[(NF * i) + j]) && (opt == 1)))
  //if (p->next == NULL)
  sum += 1;
  //else
  //	sum = sumcompVec(p->next, q->next, sum += 1, opt);
  //else
  //if (p->next != NULL)
  //	sum = sumcompVec(p->next, q->next, sum, opt);
				
  return sum;
  }
*/

unsigned int best_moo(double *fitness, double *fbests, unsigned int NF, unsigned int opt,unsigned int fun,double *pop) {
  unsigned int i, sum, ret;
	
  //  if(constraints(pop,fun) < constraints(fbests,fun))return 0;
  //if(constraints(pop,fun) > constraints(fbests,fun))return 1;
  //  cout<<"                        holas"<<endl;
  for(i = 0; i < NF; i++)
    if (((fitness[i] < fbests[i]) && (opt == 0)) || ((fitness[i] > fbests[i]) && (opt == 1)))
      sum += 1;
	
  if (sum == NF)
    ret = 0;
  else
    if (sum == 0)
      ret = 1;
    else
      ret = RandomInt(0, 1);
		
  return ret;
}

void ifbest_interchange_moo(double *fitness, double *pop, double *fbests, double *pbests, unsigned int M, unsigned int D, unsigned int NF, unsigned int opt,unsigned int fun) {
  unsigned int i, j;
		
  for(i = 0; i < M; i++) {
    if (best_moo(&fitness[NF * i], &fbests[NF * i], NF, opt,fun,pop+(NF*i))) {
      for(j = 0; j < NF; j++)
	fbests[(NF * i) + j] = fitness[(NF * i) + j];
      for(j = 0; j < D; j++)
	pbests[(D * i) + j] = pop[(D * i) + j];
    }
  }
}

unsigned int compVec(double *noDomF, unsigned int NF, unsigned int h, unsigned int i, unsigned int opt,unsigned int fun,double *noDomP,unsigned int D) {
  unsigned int sum = 0, ret, j;
	
  //sum = sumcompVec(noDomF, NF, h, i, opt);
    if(constraints(noDomF+(NF*h),fun)>0&&constraints(noDomF+(NF*i),fun)>0){
   return 4;
   }
  else if(constraints(noDomF+(NF*h),fun)>0)return 2;
  else if(constraints(noDomF+(NF*i),fun)>0)return 1;
  //if(constraints(noDomP+(D*h),fun) > constraints(pop+(D*i),fun))return 1;
  //if(constraints(noDomP+(D*h),fun) < constraints(pop+(D*i),fun))return 2;

  //  cout<<"                                          rakam"<<endl;

  for(j = 0; j < NF; j++)
    if (((noDomF[(NF * h) + j] < noDomF[(NF * i) + j]) && (opt == 0)) || ((noDomF[(NF * h) + j] > noDomF[(NF * i) + j]) && (opt == 1)))
      sum += 1;
	
  if (sum == NF)
    ret = 1;
  else
    if (sum == 0)
      ret = 2;
    else
      ret = 3;
			
  return ret;		
}
	
void delPartDom(double * noDomP, double *noDomF, unsigned int *lastPart, unsigned int D, unsigned int NF, unsigned int opt, unsigned int fun) {
  unsigned int h, i, j;// k;
  h = 0;
  i = h + 1;
  do {
    do {
      switch (compVec(noDomF, NF, h, i, opt,fun,noDomP,D)) {
      case 1: for(j = 0; j < D; j++)
	noDomP[(D * i) + j] = noDomP[(D * (*lastPart - 1)) + j];
      for(j = 0; j < NF; j++)
	noDomF[(NF * i) + j] = noDomF[(NF * (*lastPart - 1)) + j];/*
								    for(k = i; k < (*lastPart -1); k++) {
								    for(j = 0; j < D; j++)
								    noDomP[(D * k) + j] = noDomP[(D * (k + 1)) + j];
								    for(j = 0; j < NF; j++)
								    noDomF[(NF * k) + j] = noDomF[(NF * (k + 1)) + j];
								    }*/
      *lastPart -= 1;
      break;
      case 2:	 for(j = 0; j < D; j++)
	noDomP[(D * h) + j] = noDomP[(D * (*lastPart - 1)) + j];
      for(j = 0; j < NF; j++)
	noDomF[(NF * h) + j] = noDomF[(NF * (*lastPart - 1)) + j];/*
								    for(k = h; k < (*lastPart -1); k++) {
								    for(j = 0; j < D; j++)
								    noDomP[(D * k) + j] = noDomP[(D * (k + 1)) + j];
								    for(j = 0; j < NF; j++)
								    noDomF[(NF * k) + j] = noDomF[(NF * (k + 1)) + j];
								    }*/
      *lastPart -= 1;

      i = h + 1;
      break;
      case 3:	i += 1;
	break;
      case 4:
	noDomP[(D * i) + j] = noDomP[(D * (*lastPart - 1)) + j];
      for(j = 0; j < NF; j++)
	noDomF[(NF * i) + j] = noDomF[(NF * (*lastPart - 1)) + j];/*
								    for(k = i; k < (*lastPart -1); k++) {
								    for(j = 0; j < D; j++)
								    noDomP[(D * k) + j] = noDomP[(D * (k + 1)) + j];
								    for(j = 0; j < NF; j++)
								    noDomF[(NF * k) + j] = noDomF[(NF * (k + 1)) + j];
								    }*/
      *lastPart -= 1;

	noDomP[(D * h) + j] = noDomP[(D * (*lastPart - 1)) + j];
      for(j = 0; j < NF; j++)
	noDomF[(NF * h) + j] = noDomF[(NF * (*lastPart - 1)) + j];/*
								    for(k = h; k < (*lastPart -1); k++) {
								    for(j = 0; j < D; j++)
								    noDomP[(D * k) + j] = noDomP[(D * (k + 1)) + j];
								    for(j = 0; j < NF; j++)
								    noDomF[(NF * k) + j] = noDomF[(NF * (k + 1)) + j];
								    }*/
      *lastPart -= 1;

      i = h + 1;

	break;

      }
    } while (i < *lastPart);