粒子群算法求解优化问题(c实现)

#include"stdio.h"
#include "math.h"
#include "stdlib.h"
#include "time.h"
#define N 100  //最大粒子数 
#define D 10  //最大问题维数 

double fitness(double x[],int n);
typedef double (*funType)(double [],int n);
double pso(funType fp,int n,double c1,double c2,double w,int m,int d,double xbest[]);

int main(){
 double x[D];
 double f;
 f = pso(fitness,40,1.5,1.5,0.6,2000,1,x);
 printf("最优的x为：");
 for(int i=0;i<1;i++)
  printf("%f ",x[i]); 
 printf("\n");
 printf("函数的最优值为：%f\n",f);
 return 0;

}
double fitness(double x[],int n){  //待优化函数 
 double f = 0;
 f = x[0]*x[0]-2*x[0]+1;  // f(x) = x.^2-2*x+1
 return f; 
}

/*
 fp: 待优化的目标函数
 n: 粒子数目
 c1,c2: 学习因子1，学习因子2
 w: 惯性权重
 M: 最大迭代次数
 d: 问题的维数
 x: 返回的最优解 
*/
double pso(funType fp,int n,double c1,double c2,double w,int m,int d,double xbest[]){
 int i,j,t;
 double x[N][D],v[N][D];
 double pg[D]; //全局最优解
 double y[N];  //粒子适应度
 double p[N][D];  
 //step1.初始化种群的个体
 srand((int)time(0)); 
 for(i=0;i<n;i++)           
  for(j=0;j<d;j++){
   x[i][j] = rand()*1.0/RAND_MAX;
   v[i][j] = rand()*1.0/RAND_MAX;
  }
 //step2. 计算粒子适应度，初始化p和pg 
 for(i=0;i<d;i++)
  pg[i] = x[n-1][i]; // pg为全局最优 
 for(i=0;i<n;i++){
  y[i] = fp(x[i],d);
  for(j=0;j<d;j++)
   p[i][j] = x[i][j];
 }
 for(i=0;i<n-1;i++){
  if(y[i]<y[n-1]){
   for(j=0;j<d;j++)
    pg[j] = x[i][j];
  }  
 }
 //step3. 进入主循环，按照公式依次迭代  
 for(t=0;t<m;t++){         
  for(i=0;i<n;i++){
   for(j=0;j<d;j++){
    v[i][j] = w*v[i][j] + c1*(rand()/RAND_MAX)*(p[i][j]-x[i][j]) + c2*(rand()/RAND_MAX)*(pg[j]-x[i][j]); //更新 
    x[i][j] = x[i][j] + v[i][j];
   }
   if(fp(x[i],d)<y[i]){
    y[i] = fp(x[i],d);
    for(j=0;j<d;j++)
     p[i][j] = x[i][j];
   }
   if(y[i]<fp(pg,d)){
    for(j=0;j<d;j++)
     pg[j] = p[i][j];
   }
  } 
  
 }
 
 for(i=0;i<d;i++) //最优解 
  xbest[i] = pg[i];  
 return fp(pg,d);
}