/*該程序用PSO算法求解約束性優化問題， */ 
 
//庫文件 
#include"stdio.h" 
#include"stdlib.h" 
#include"math.h" 

//隨機數定義 ,r1 ~U(0,1)，r2 ~U(0,1)為兩個相互獨立的隨機函數
#define rdint(i) (rand()%(int)(i)) 
#define rdft() (float)((double)rdint(16384)/(16383.0)) 
#define rnd(a,b) (rdint((int)(b)-(int)(a)+1)+(int)(a)) 
 
//宏定義 
#define POPSIZE 20 //粒子群規模
#define DIMENSION 2 //微粒的多維向量的維數為3
 
//全局變量定義 
float SHK=0.73;//收縮因子
float W=0.9; //權重
float C1=1.8; //c1、c2為加速常數，通常在0~2間取值。
float C2=1.8; 
float VMAX=2; /*為了減少在進化過程中，微粒離開搜索空間的可能性，vij通常限定于一定
			  范圍內，即v∈[-vmax ,vmax]。如果問題的搜索空間限定在[-xmax ，xmax ]內，則可設定v =k·x ，
              0.1≤k≤1.0。*/ 
float XMIN=-2.0; //搜索空間限定,-2<x1,x2<2
float XMAX=2.0; 
float P[DIMENSION]; //全局最好位置，對應該Goldstein-Price函數最小值的x1,x2取值
float PBEST; //全局最佳適應值，對應該Goldstein-Price函數的函數最小值
 
struct individual //定義粒子屬性
{ 
 float number[DIMENSION]; //個體當前位置
 float best[DIMENSION]; //個體最好位置
 float bestfitness; //個體最佳適應值
 float fitness; //個體當前適應值
 float speed[DIMENSION]; //個體當前速度
}individual[POPSIZE]; 
 
void initiate(void); 
void calculation(int number); 
void globalbest(int number); 
void localbest(int number); 
 
//程序初始化定義 
void initiate() 
{ 
 int i,j; 
 for(i=0;i<POPSIZE;i++) //對每個粒子的位置進行初始化
  for(j=0;j<DIMENSION;j++) 
   individual[i].number[j]=rdft()*(XMAX-XMIN)+XMIN; 
 for(i=0;i<POPSIZE;i++) //對每個粒子的速度進行初始化
  for(j=0;j<DIMENSION;j++) 
   individual[i].speed[j]=VMAX*rdft(); 
 for(i=0;i<POPSIZE;i++) //粒子初始化時，其初始化位置即設為個體最佳位置
  for(j=0;j<DIMENSION;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; 
 if(individual[number].bestfitness>individual[number].fitness) 
  for(i=0;i<DIMENSION;i++) 
   individual[number].best[i]=individual[number].number[i]; 
 individual[number].bestfitness=individual[number].fitness; 
} 
 
//種群歷史最優位置修改程序 
void globalbest(int number) 
{ 
 int i,j; 
 float 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<DIMENSION;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<DIMENSION;j++) 
     P[j]=individual[i].best[j]; 
    PBEST=individual[i].bestfitness; 
   } 
 } 
} 
 
//適應值函數計算程序 
void calculation(int num) 
{ //計算編號為“nmu”的粒子的當前適應值
 int i; 
 float s=0.0,h=0.0; 
 s=pow(individual[num].number[0]+individual[num].number[1]+1.0,2.0)*(19.0-14.0*individual[num].number[0]+3.0*pow(individual[num].number[0],2.0)-14.0*individual[num].number[1]+6.0*individual[num].number[0]*individual[num].number[1]+3.0*pow(individual[num].number[1],2.0))+1.0; 
 h=30.0+pow(2.0*individual[num].number[0]-3.0*individual[num].number[1],2.0)*(18.0-32.0*individual[num].number[0]+12.0*pow(individual[num].number[0],2.0)+48.0*individual[num].number[1]-36.0*individual[num].number[0]*individual[num].number[1]+27.0*pow(individual[num].number[1],2.0)
); 
 individual[num].fitness=s*h; 
} 
 
//主程序 
main() 
{ 
 int i,j,k,t,total=0; 
 float sum=0,Total_fit=0; 
 for(j=0;j<100;j++) //進行100次的試驗
 { 
  initiate(); //將粒子群初始化
      for(i=0;i<1000;i++) //一次優化過程中最大的迭代次數是1000
  { 
   W=0.9-i*0.5/1000; /*設置權重，將W設定為隨著進化而線性減少，這樣使粒子在前期有較高的探索能力以得到合適的種子，而在后期有較高的開發能力以加快收斂速度*/
   if((PBEST-3.0)<0.0001) 
   { 
    total++; //統計試驗成功（即達到迭代終止條件）的次數，total<=100.
    sum=sum+i; //統計每次試驗（共100次）之后所進行的總的迭代次數。
	Total_fit=Total_fit+PBEST; 
    break; //當全局最佳適應值PBEST達到迭代終止條件時跳出循環。
   } 
   else 
   { 
    for(k=0;k<POPSIZE;k++) 
    { 
     for(t=0;t<DIMENSION;t++) 
     {  //一次迭代過程，根據進化方程進行優化
		 individual[k].speed[t]=SHK*(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]=2*XMAX-individual[k].number[t]; 
	 } 
	 calculation(k); //一個粒子一次進化結束后，計算該粒子的適應值individual[num].fitness
     localbest(k); //一個粒子一次進化結束后，計算該粒子個體最佳適應值individual[number].bestfitness
	} 
	globalbest(1); //所以粒子進行一次進化后，計算全局最佳適應值PBEST
   } 
  } 
	  
		  printf("迭代次數：%d,   函數最佳適應值：%f\n",i,PBEST); //一次試驗結束后后，計算迭代次數i和最小值PBEST
 } 
  printf("\n試驗成功概率：%d(共100次)\n",total); 
  printf("總的迭代次數：%f\n",sum); 
  printf("平均迭代次數=總的迭代次數/試驗成功次數=%f\n",sum/total); 
  printf("總的最佳適應值之和：%f\n",Total_fit);
  printf("平均最佳適應值=總的適應值之和/試驗成功次數=%f\n",Total_fit/total);
}