粒子群python3实现

粒子群优化算法(详细易懂_很多例子)

在实现的过程中由于对其算法理解不够遇到了一些小问题。于是又找到了上面这篇文章加深了理解。

#实现过程

##1.初始化
（1）种群规模，即代码中的pop,为整数
（2）取值范围lower，upper
（3）初始化位置，x，求得对于适应值fitness,当前历史最优值pbest_x,pbest_y就设置为对应的x和fitness，当前全局最优gbest_x,gbest_y设置为pbest_x和pbest_y对应的最优位置。初始化速度v为随机值.
##2.根据fitness function，评价每个粒子的适应度。
##3.构建速度迭代
v[i] = w * v[i] + c1 * np.random.rand() * (pbest_x[i] - x[i]) + c2np.random.rand()(gbest_x - x[i])
##4.构建位置迭代x = x + v
##5.在新的位置和速度下寻找每个粒子新的当前历史最优位置和最优pbest_x,pbest_y。
##6.在新的位置和速度下寻找所有粒子新的全局最优位置和最优值gbest_x,gbest_y。
##7.设置最大迭代次数或者最佳适应度值的增量小于某个给定的阈值时算法停止，否则返回第二步。

#实现代码
这里设计了y=x1^2 + x2^2 + …+ xn^2（n=1,2,3…）作为fitness function。

import numpy as np
import matplotlib.pyplot as plt
#from mpl_toolkits.mplot3d import Axes3D

#初始化种群
def init_pop(pop=10, dim=3, lower=-300, upper=300):
    
    x = np.random.uniform(lower, upper, size=(pop,dim))
    #x = np.array([300,300]).reshape(2,1)
    y = fitness(x)
    v = np.random.uniform(lower, upper, size=(pop,dim))
    pbest_x = x
    pbest_y = y
    gbest_x = x[np.argmin(y)]
    gbest_y = y[np.argmin(y)]
    return  x,v,y,pbest_x, pbest_y, gbest_x, gbest_y

def fitness(x):
    value = []
    for i in x:
        value.append(np.dot(i.T, i)) 
    return  np.array(value).reshape(-1,1)


def update_pbest(x,pbest_x,pbest_y):
    new_fitness = fitness(x)
    for i in range(len(x)):
        if pbest_y[i] > new_fitness[i]:
            pbest_y[i] = new_fitness[i]
            pbest_x[i] = x[i]
    return pbest_x,pbest_y



def update_gbest(x,gbest_x, gbest_y):
    new_fitness = fitness(x)
    for i in new_fitness:
        if i < gbest_y:
            gbest_y = i
            gbest_x = x[np.where(new_fitness == gbest_y)]
    return gbest_x, gbest_y



def update_v(v, x, pbest_x, gbest_x, w=0.2, c1=2, c2=2, dim=3,v_max=30):
    for i in range(len(v)):
        v[i] = w * v[i] + c1 * np.random.rand() * (pbest_x[i] - x[i]) + c2*np.random.rand()*(gbest_x - x[i])
        for j in range(len(v[i])):
            if v[i][j] > v_max:
                v[i][j] = v_max
            elif v[i][j] < -v_max:
                v[i][j] = -v_max
    return v


def update_x(v,x):
    x = x + v
    return x


def run(max_iter=100):
    i = 0
    n_dim = 1
    x,v,y,pbest_x, pbest_y, gbest_x, gbest_y = init_pop(pop=1000, C1=2, C2=2, dim=n_dim, lower=-300, upper=300)
    y_best_list = []
    while i <= max_iter:
        plt.clf()
        scatter_x = x
        scatter_y = y
        x1 = gbest_x
        y1 = gbest_y
        plt.scatter(scatter_x, scatter_y, c='b')
        plt.scatter(x1, y1, c='r')
        plt.xlim(-300, 300)
        plt.savefig(str(i)+'.png')
        plt.pause(0.01)
        v = update_v(v,x, pbest_x, gbest_x,dim=n_dim)
        x = update_x(v,x)
        pbest_x,pbest_y = update_pbest(x, pbest_x, pbest_y)
        gbest_x, gbest_y = update_gbest(x, gbest_x, gbest_y)
        #print("更新第"+str(i)+"次,最佳位置在",gbest_x,"目前最优值为",gbest_y)
        i = i+1
        y_best_list.append(gbest_y)
    return None    

# =============================================================================
# =============================================================================
#  def run(max_iter=100):
#      i = 0
#      n_dim = 2
#      x,v,y,pbest_x, pbest_y, gbest_x, gbest_y = init_pop(pop=100, dim=n_dim, lower=-300, upper=300)
#      y_best_list = []
#      while i <= max_iter:
#          plt.clf()
#          scatter_x = x[:,0]
#          scatter_y = x[:,1]
#          z = y
#          ax = plt.axes(projection='3d')
#          ax.scatter3D(scatter_x,scatter_y,z,cmap='Blues')
#          plt.show()
#          v = update_v(v,x, pbest_x, gbest_x,dim=n_dim)
#          x = update_x(v,x)
#          pbest_x,pbest_y = update_pbest(x, pbest_x, pbest_y)
#          gbest_x, gbest_y = update_gbest(x, gbest_x, gbest_y)
#          print("更新第"+str(i)+"次,最佳位置在",gbest_x,"目前最优值为",gbest_y)
#          i = i+1
#          y_best_list.append(gbest_y)
#      return None
# 
# =============================================================================
y_best = run(max_iter=100)

下图是上述结果，可以看出所有的粒子都在不断地往（0，0）点靠近。