粒子群算法测试程序(PSO)

// PSO.cpp : Defines the entry point for the console application.
//

#include "stdafx.h"

#include"stdio.h"
#include"stdlib.h"
#include"time.h"
#include"math.h"

#define rdint(i) (rand()%(int)(i))/*随机数的定义*/
#define rdft() (double)((double)rdint(16384)/(16383.0))
#define rnd(a,b) (rdint((int)(b)-(int)(a)+1)+(int)(a))

#define POPSIZE 20
#define DIMENSION1 2   /*  目标函数中所含变量x的个数  */
#define DIMENSION2 6   //约束条件的个数

double W=1.0;
double C1=1.8;
double C2=1.8;
double VMAX=3.0;
double XMIN=2.0;
double XMAX=4.0;
double P[DIMENSION1];
double PBEST;
double H,K;
double g[DIMENSION2];//约束条件的值

struct indi
{
 double number[DIMENSION1];
 double best[DIMENSION1];
 double bestfitness;
 double fitness;
 double speed[DIMENSION1];
}individual[POPSIZE];

void initiate(void);
void calculation(int number);
void globalbest(int number);
void localbest(int number);
void judge(int num);

void judge(int num)
{

 /*     6个约束条件     */
 g[0]=4*individual[num].number[0]-individual[num].number[1];
 g[1]=individual[num].number[1]-10*individual[num].number[0];
 g[2]=-2*K*individual[num].number[1]+3*3.14*pow(individual[num].number[0],3);
 g[3]=-50*K*individual[num].number[1]+75*3.14*pow(individual[num].number[0],3);  
 g[4]=-H+3.98*individual[num].number[1];
 g[5]=3*pow(individual[num].number[1],3)-2000*pow(individual[num].number[0],4);
}
void initiate()
{
 int i,j;

 for(i=0;i<POPSIZE;i++)
 {
    do
    {
     for(j=0;j<DIMENSION1;j++)
    individual[i].number[j]=rdft()*(XMAX-XMIN)+XMIN;
   judge(i);

  }while(g[0]<=0 && g[1]<=0 && g[2]<=0 && g[3]<=0 && g[4]<=0 && g[5]<=0);  
 }//如果不在可行域就对这个粒子重新进行初始化
 
 for(i=0;i<POPSIZE;i++)
  for(j=0;j<DIMENSION1;j++)
   individual[i].speed[j]=VMAX*rdft();
 for(i=0;i<POPSIZE;i++)
  for(j=0;j<DIMENSION1;j++)
   individual[i].best[j]=individual[i].number[j];
 for(i=0;i<POPSIZE;i++)
  calculation(i);
 for(i=0;i<POPSIZE;i++)
  individual[i].bestfitness=individual[i].fitness;
 globalbest(0);
}

void localbest(int number)
{
 int i;
 judge(number);//进行判断，如果不在可行域内就不对粒子的历史最优值进行更新
 bool bTrue=g[0]<=0 && g[1]<=0 && g[2]<=0 && g[3]<=0 && g[4]<=0 && g[5]<=0;
 if(individual[number].bestfitness>individual[number].fitness&&bTrue)
 {
  for(i=0;i<DIMENSION1;i++)
   individual[number].best[i]=individual[number].number[i];
  individual[number].bestfitness=individual[number].fitness;
 }
}

void globalbest(int number)
{
 int i,j;
 double s=0;
 int flag=0;
 if(number==0)
 {
  s=individual[0].fitness;
  flag=0;
  for(i=1;i<POPSIZE;i++)
   if(individual[i].fitness<s)
   {
    s=individual[i].fitness;
    flag=i;
   }
   for(i=0;i<DIMENSION1;i++)
    P[i]=individual[flag].number[i];
   PBEST=individual[flag].fitness;
 }
 else
 {
  for(i=0;i<POPSIZE;i++)
  {
   
   if(individual[i].bestfitness<PBEST)
   {
    for(j=0;j<DIMENSION1;j++)
     P[j]=individual[i].best[j];
    PBEST=individual[i].bestfitness;
   }
  }
 }
}

void calculation(int num)
{
 K=(4*individual[num].number[1]-individual[num].number[0])/4*(individual[num].number[1]-individual[num].number[0])
  +0.615*individual[num].number[0]/individual[num].number[1];
 H=(2000*pow(individual[num].number[0],5)/pow(individual[num].number[1],3))+1.5*individual[num].number[0]
 +75*3.14*pow(individual[num].number[0],3)/(4*K);
 //L=2000*3.14*pow(individual[num].number[0],4)/pow(individual[num].number[1],2)+1.5*3.14*individual[num].number[1];
 individual[num].fitness=0.25*3.14*pow(individual[num].number[1],2)*H;  
}

void main(void)
{
 int i,j,k,t,total=0;
 double sum=0;
 // srand(time(NULL));
 for(j=0;j<10;j++)
 {
  initiate();
  for(i=0;i<5000;i++)
  {
   W=1.0-i*0.6/4999;
   if((PBEST-67997.0)<0.1)
   {
    total++;
    sum=sum+i;
    break;
   }
   else
   {
    for(k=0;k<POPSIZE;k++)
    {
     for(t=0;t<DIMENSION1;t++)
     { 
      individual[k].speed[t]=W*individual[k].speed[t]+C1*rdft()*(individual[k].best[t]-individual[k].number[t])+C2*rdft()*(P[t]-individual[k].number[t]);
      if(individual[k].speed[t]>VMAX)
       individual[k].speed[t]=VMAX;
      individual[k].number[t]=individual[k].number[t]+individual[k].speed[t];
      if(individual[k].number[t]<XMIN)
       individual[k].number[t]=2*XMIN-individual[k].number[t];
      if(individual[k].number[t]>XMAX)      
       individual[k].number[t]=fmod(individual[k].number[t],(XMAX-XMIN))+XMIN;
     }

     calculation(k);
     localbest(k);
    }
    globalbest(1);
   }
  }
  //  printf("%3d,%f/n",i,PBEST);
 }
 printf("Total number is  %d/n",total);
 printf("Sum is %f/n",sum);
 printf("Average is %f/n",sum/total);
 printf("%f/n",PBEST);
 getchar();
}