Python之粒子群演算法（含程式碼例項）

#f(x,y) = x^2 + y^2 + x

import numpy as np

import matplotlib.pyplot as plt

from matplotlib import cm

from mpl_toolkits.mplot3d import Axes3D# 匯入該函式是為了繪製 3D 圖

import matplotlib as mpl

#########

# 將資料繪圖出來

# 生成 X 和 Y 的資料

X = np.arange(-5,5,0.1) #-5 到 5 的等距陣列，距離為 0.1 ，注意區分開 range(), 它返回的是一個列表

Y = np.arange(-5,5,0.1)

X,Y = np.meshgrid(X,Y)# 該函式用來生成網格點座標矩陣。

# 目標函式

Z = X**2 + Y**2 + X

# 繪圖

fig = plt.figure()# 創立一個畫布

ax = Axes3D(fig)# 在這個畫布裡，生成一個三維的空間

surf = ax.plot_surface(X,Y,Z,cmap=cm.coolwarm)# 該函式是為了將資料在這三維空間裡視覺化出來。

plt.show()

###########

#######

# 計算適應度，這裡的適應度就是我們目標函式 Z 的值，因為我們要求 Z 的最小值。

# 這兩個函式，一般使用 mpl 畫圖的時候都會用到

mpl.rcParams['font.sans-serif'] = ['SimHei']# 指定預設字型

mpl.rcParams['axes.unicode_minus'] = False # 解決儲存影像是負號 '-' 無法顯示的問題

# 使用 matplotliblib 畫圖的時候經常會遇見中文或者是負號無法顯示的情況

#rcParams 函式里的引數可以修改預設的屬性，包括窗體大小、每英寸的點數、線條寬度、顏色、樣式、座標軸、座標和網路屬性、文字、字型等

def fitness_func(X):

     x = X[:,0]

     y = X[:,1]

     return x**2 + y**2 + x

##########

#############

# 更新速度，根據公式 V(t+1)=w*V(t)+c1*r1*(pbest_i-xi)+c1*r1*(gbest_xi)

def velocity_update(V,X,pbest,gbest,c1,c2,w,max_val):

     size = X.shape[0]# 返回矩陣 X 的行數

     r1 = np.random.random((size,1))# 該函式表示成 size 行 1 列的浮點數，浮點數都是從 0-1 中隨機。

     r2 = np.random.random((size,1))

     V = w*V + c1*r1*(pbest-X)+c2*r2*(gbest-X)# 注意這裡得到的是一個矩陣

     # 這裡是一個防止速度過大的處理，怕錯過最理想值

     V[V<-max_val] = -max_val

     V[V>-max_val] = max_val

     return V

#########

########

# 更新粒子位置，根據公式 X(t+1)=X(t)+V

def position_updata(X,V):

     return X+V

##########

#########

def pos():

     w = 1 # 設定慣性權重

     c1 = 2 # 設定個體學習係數

     c2 = 2 # 設定全域性學習係數

     r1 = None

     r2 = None

     dim = 2

     size = 20 # 這裡是初始化粒子群， 20 個

     iter_num = 1000 # 迭代 1000 次

     max_val = 0.5 # 限定最大速度為 0.5

     best_fitness = float(9e10) # 初始化適應度的值

     fitness_val_list = []

     # 初始化各個粒子的位置

     X = np.random.uniform(-5,5,size=(size,dim))

     # 初始化各個粒子的速度

     V = np.random.uniform(-0.5,0.5,size=(size,dim))

     p_fitness = fitness_func(X)# 得到各個個體的適應度值

     g_fitness =外匯跟單gendan5.com p_fitness.min()# 全域性最理想的適應度值

     fitness_val_list.append(g_fitness)

     pbest = X# 初始化個體的最優位置

     gbest = X[p_fitness.argmin()]# 初始化整個整體的最優位置

     # 迭代

     for i in range(1,iter_num):

         V = velocity_update(V, X, pbest, gbest, c1, c2, w, max_val)

         X = position_updata(X,V)

         p_fitness2 = fitness_func(X)

         g_fitness2 = p_fitness2.min()

         # 更新每個粒子的歷史的最優位置

         for j in range(size):

             if p_fitness[j] > p_fitness2[j]:

                 pbest[j] = X[j]

                 p_fitness[j] = p_fitness2[j]

             if g_fitness > g_fitness2:

                 gbest = X[p_fitness2.argmin()]

                 g_fitness = g_fitness2

             fitness_val_list.append(g_fitness)

             i += 1

     print(" 最優值是： %.5f" % fitness_val_list[-1])

     print(" 最優解是： x=%.5f,y=%.5f" % (gbest[0],gbest[1]))

     plt.plot(fitness_val_list,c='r')

     plt.title(' 迭代過程 ')

     plt.show()

#########

if __name__ == '__main__':

     pos()