粒子群演算法求解帶約束最佳化問題 原始碼實現

步驟1 ：初始化引數

import numpy as np

import matplotlib.pyplot as plt

import matplotlib as mpl

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

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

# PSO 的引數

w = 1  # 慣性因子，一般取 1

c1 = 2  # 學習因子，一般取 2

c2 = 2  #

r1 = None  # 為兩個（ 0,1 ）之間的隨機數

r2 = None

dim = 2  # 維度的維度 # 對應 2 個引數 x,y

size = 100  # 種群大小，即種群中小鳥的個數

iter_num = 1000  # 演算法最大迭代次數

max_vel = 0.5  # 限制粒子的最大速度為 0.5

fitneess_value_list = []  # 記錄每次迭代過程中的種群適應度值變化

步驟2 ：這裡定義一些引數，分別是計算適應度函式和計算約束懲罰項函式

def calc_f(X):

    """ 計算粒子的的適應度值，也就是目標函式值， X 的維度是 size * 2 """

    A = 10

    pi = np.pi

    x = X[0]

    y = X[1]

    return 2 * A + x ** 2 - A * np.cos(2 * pi * x) + y ** 2 - A * np.cos(2 * pi * y)

def calc_e1(X):

    """ 計算第一個約束的懲罰項 """

    e = X[0] + X[1] - 6

    return max(0, e)

def calc_e2(X):

    """ 計算第二個約束的懲罰項 """

    e = 3 * X[0] - 2 * X[1] - 5

    return max(0, e)

def calc_Lj(e1, e2):

    """ 根據每個粒子的約束懲罰項計算 Lj 權重值， e1, e2 列向量，表示每個粒子的第 1 個第 2 個約束的懲罰項值 """

    # 注意防止分母為零的情況

    if (e1.sum() + e2.sum()) <= 0:

        return 0, 0

    else:

        L1 = e1.sum() / (e1.sum() + e2.sum())

        L2 = e2.sum() / (e1.sum() + e2.sum())

    return L1, L2

步驟3 ：定義粒子群演算法的速度更新函式，位置更新函式

def velocity_update(V, X, pbest, gbest):

    """

     根據速度更新公式更新每個粒子的速度

      種群size=20

    :param V: 粒子當前的速度矩陣， 20*2 的矩陣

    :param X: 粒子當前的位置矩陣， 20*2 的矩陣

    :param pbest: 每個粒子歷史最優位置， 20*2 的矩陣

    :param gbest: 種群歷史最優位置， 1*2 的矩陣

    """

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

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

    V = w * V + c1 * r1 * (pbest - X) + c2 * r2 * (gbest - X)  # 直接對照公式寫就好了

    # 防止越界處理

    V[V < -max_vel] = -max_vel

    V[V > max_vel] = max_vel

    return V

def position_update(X, V):

    """

     根據公式更新粒子的位置

    :param X: 粒子當前的位置矩陣，維度是 20*2

    :param V: 粒子當前的速度舉著，維度是 20*2

    """

    X=X+V# 更新位置

    size=np.shape(X)[0]# 種群大小

    for i in range(size):# 遍歷每一個例子

        if X[i][0]<=1 or X[i][0]>=2:#x 的上下限約束

            X[i][0]=np.random.uniform(1,2,1)[0]# 則在 1 到 2 隨機生成一個數

        if X[i][1] <= -1 or X[i][0] >= 0:#y 的上下限約束

            X[i][1] = np.random.uniform(-1, 0, 1)[0]  # 則在 -1 到 0 隨機生成一個數

    return X

步驟4 ：每個粒子歷史最優位置更優函式，以及整個群體歷史最優位置更新函式，和無約束約束最佳化程式碼類似，所不同的是新增了違反約束的處理過程

def update_pbest(pbest, pbest_fitness, pbest_e, xi, xi_fitness, xi_e):

    """

     判斷是否需要更新粒子的歷史最優位置

    :param pbest: 歷史最優位置

    :param pbest_fitness: 歷史最優位置對應的適應度值

    :param pbest_e: 歷史最優位置對應的約束懲罰項

    :param xi: 當前位置

    :param xi_fitness: 當前位置的適應度函式值

    :param xi_e: 當前位置的約束懲罰項

    :return:

    """

    # 下面的 0.0000001 是考慮到計算機的數值精度位置，值等同於 0

    # 規則 1 ，如果 pbest 和 xi 都沒有違反約束，則取適應度小的

    if pbest_e <= 0.0000001 and xi_e <= 0.0000001:

        if pbest_fitness <= xi_fitness:

            return pbest, pbest_fitness, pbest_e

        else:

            return xi, xi_fitness, xi_e

    # 規則 2 ，如果當前位置違反約束而歷史最優沒有違反約束，則取歷史最優

    if pbest_e < 0.0000001 and xi_e >= 0.0000001:

        return pbest, pbest_fitness, pbest_e

    # 規則 3 ，如果歷史位置違反約束而當前位置沒有違反約束，則取當前位置

    if pbest_e >= 0.0000001 and xi_e < 0.0000001:

        return xi, xi_fitness, xi_e

    # 規則 4 ，如果兩個都違反約束，則取適應度值小的

    if pbest_fitness <= xi_fitness:

        return pbest, pbest_fitness, pbest_e

    else:

        return xi, xi_fitness, xi_e

def update_gbest(gbest, gbest_fitness, gbest_e, pbest, pbest_fitness, pbest_e):

    """

     更新全域性最優位置

    :param gbest: 外匯跟單gendan5.com 上一次迭代的全域性最優位置

    :param gbest_fitness: 上一次迭代的全域性最優位置的適應度值

    :param gbest_e: 上一次迭代的全域性最優位置的約束懲罰項

    :param pbest: 當前迭代種群的最優位置

    :param pbest_fitness: 當前迭代的種群的最優位置的適應度值

    :param pbest_e: 當前迭代的種群的最優位置的約束懲罰項

    :return:

    """

    # 先對種群，尋找約束懲罰項 =0 的最優個體，如果每個個體的約束懲罰項都大於 0 ，就找適應度最小的個體

    pbest2 = np.concatenate([pbest, pbest_fitness.reshape(-1, 1), pbest_e.reshape(-1, 1)], axis=1)  # 將幾個矩陣拼接成矩陣 ， 4 維矩陣（ x,y,fitness,e ）

    pbest2_1 = pbest2[pbest2[:, -1] <= 0.0000001]  # 找出沒有違反約束的個體

    if len(pbest2_1) > 0:

        pbest2_1 = pbest2_1[pbest2_1[:, 2].argsort()]  # 根據適應度值排序

    else:

        pbest2_1 = pbest2[pbest2[:, 2].argsort()]  # 如果所有個體都違反約束，直接找出適應度值最小的

    # 當前迭代的最優個體

    pbesti, pbesti_fitness, pbesti_e = pbest2_1[0, :2], pbest2_1[0, 2], pbest2_1[0, 3]

    # 當前最優和全域性最優比較

    # 如果兩者都沒有約束

    if gbest_e <= 0.0000001 and pbesti_e <= 0.0000001:

        if gbest_fitness < pbesti_fitness:

            return gbest, gbest_fitness, gbest_e

        else:

            return pbesti, pbesti_fitness, pbesti_e

    # 有一個違反約束而另一個沒有違反約束

    if gbest_e <= 0.0000001 and pbesti_e > 0.0000001:

        return gbest, gbest_fitness, gbest_e

    if gbest_e > 0.0000001 and pbesti_e <= 0.0000001:

        return pbesti, pbesti_fitness, pbesti_e

    # 如果都違反約束，直接取適應度小的

    if gbest_fitness < pbesti_fitness:

        return gbest, gbest_fitness, gbest_e

    else:

        return pbesti, pbesti_fitness, pbesti_e