#给定初始化条件
c1 <- 2            #学习因子1
c2 <- 2            #学习因子2
w <- 0.7           #惯性权重
MaxDT <- 1000      #最大迭代次数
D <- 10            #搜索空间维数（未知数个数）
M <- 40            #初始化群体个体数目

#初始化种群的个体(可以在这里限定位置和速度的范围)
x <- matrix(rnorm(M*D),nrow = M,ncol = D,byrow = T)
v <- matrix(rnorm(M*D),nrow = M,ncol = D,byrow = T)

#先计算各个粒子的适应度，并初始化p[i]和gbest
fitness <- function(x)
{
  return(sum(x^2))
}

p <- apply(x,1,fitness)
y <- x
gbest <- x[which(p == min(p)),]            #gbest为全局最优

#进入主要循环，按照公式依次迭代，直到满足精度要求
for (t in 1:MaxDT)
{
  for (i in 1:M)
  {
    v[i, ] = w*v[i, ] + c1*runif(1)*(y[i, ] - x[i, ]) + c2*runif(1)*(gbest-x[i, ])
    x[i, ] = x[i, ] + v[i, ]
    if(fitness(x[i, ]) < p[i])
    {
      p[i] <- fitness(x[i, ])
      y[i, ] <- x[i, ]
    }
    if(p[i] < fitness(gbest))
    {
      gbest <- y[i, ]
    }
  }
}

gbest
fitness(gbest)