算法原理

之前求解的無約束的問題。
粒子群算法求解無約束優化問題 源碼實現

算法原理如下

??????今天講解下求解約束優化的問題。該問題常用的方法是罰函數法。即如果一個解x不滿足約束條件，就對適應度值設置一個懲罰項。它的思想類似線性規劃內點法，都是通過增加罰函數，迫使模型在迭代計算的過程中始終在可行域內尋優。
假設有如下帶約束的優化問題

??????問題中既含有等式約束，也含有不等式約束，當一個解x不滿足約束條件時，則會對目標函數增加懲罰項，這樣就把帶約束優化問題變成無約束優化問題，新的目標函數如下：

其中σ是懲罰因子，一般取σ=t√t ，P(x)是整體懲罰項，P(x)的計算方法如下：

對于不等式gi(x)<=0,懲罰項

即如果gi(x)<=0,則ei(x)=0,否則ei(x)=-gi(x)

對于等式約束需要先轉換成不等式約束，一個簡單的方法設定一個等式約束容忍度值ε,新的不等式約束為

因此等式約束的懲罰項為

??????整體懲罰性P(x)是各個約束懲罰項的和，即：

??????由于約束條件之間的差異，某些約束條件可能對個體違反程度起著決定性的作用，為了平衡這種差異，對懲罰項做歸一化處理，下面的公式推導將不區分等式約束和不等式約束，在實際處理中做區分即可：

??????其中Lj表示每個約束條件的違背程度，m為約束條件的個數。公式中分子的意思是，對每個粒子xi計算違反第j個約束的懲罰項，然后求和；分母的意思：對每個粒子xi計算違背每個約束的懲罰項，然后求和，因此Lj也可以看成是第j個約束懲罰項的權重。
最后得到的粒子xi的整體懲罰項P(x)的計算公式：

?????? 在粒子群算法中，每一步迭代都會更新Pbest和Gbest，雖然可以將有約束問題轉換為無約束問題進行迭代求解，但是問題的解xi依然存在不滿足約束條件的情況，因此需要編制一些規則來比較兩個粒子的優劣，規則如下：

如果兩個粒子xi和xj都可行，則比較其適應度函數f(xi)和f(xj)，值小的粒子為優。

當兩個粒子xi和xj都不可行，則比較懲罰項P(xi)和P(xj），違背約束程度小的粒子更優。

當粒子xi可行而粒子xj不可行，選可行解。

?????? 對于粒子的上下限約束 可以體現在位置更新函數里，不必加懲罰項。 具體思路就是遍歷每一個粒子的位置，如果超除上下限，位置則更改為上下限中的任何一個位置

算例

語言python3.7
問題如下：

此函數圖像

?????? 為方便粒子群算法的迭代寫代碼，肯定要如下格式矩陣，用來存儲粒子群每個粒子xi的歷史最優位置Pbest、總的適應度值fitness、目標函數的值、約束1的懲罰項e1、約束2的懲罰項e2

粒子序號Pbest_fitnessPbest_efitnessfe1e2e

x1	l歷史最優位置對應的適應度	歷史最優位置對應的懲罰項	當前的適應度fitness=f+e	當前目標函數值	約束1的懲罰項e1	約束2的懲罰項e2	懲罰項的和e=L1e1+L2e2
x2
…

步驟1：初始化參數

import numpy as np import matplotlib.pyplot as plt import matplotlib as mplmpl.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 = 10pi = np.pix = 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] - 6return max(0, e)def calc_e2(X):"""計算第二個約束的懲罰項"""e = 3 * X[0] - 2 * X[1] - 5return max(0, e)def calc_Lj(e1, e2):"""根據每個粒子的約束懲罰項計算Lj權重值，e1, e2列向量，表示每個粒子的第1個第2個約束的懲罰項值"""# 注意防止分母為零的情況if (e1.sum() + e2.sum()) <= 0:return 0, 0else: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_velV[V > max_vel] = max_velreturn Vdef 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_eelse: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_eelse:return xi, xi_fitness, xi_edef update_gbest(gbest, gbest_fitness, gbest_e, pbest, pbest_fitness, pbest_e):"""更新全局最優位置:param gbest: 上一次迭代的全局最優位置: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_eelse:return pbesti, pbesti_fitness, pbesti_e# 有一個違反約束而另一個沒有違反約束if gbest_e <= 0.0000001 and pbesti_e > 0.0000001:return gbest, gbest_fitness, gbest_eif 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_eelse:return pbesti, pbesti_fitness, pbesti_e

步驟5：PSO

流程圖如圖：

約束體現在更新粒子最優 ，和全局最優上。

# 初始化一個矩陣 info, 記錄： # 0、種群每個粒子的歷史最優位置對應的適應度， # 1、歷史最優位置對應的懲罰項， # 2、當前適應度， # 3、當前目標函數值， # 4、約束1懲罰項， # 5、約束2懲罰項， # 6、懲罰項的和 # 所以列的維度是7 info = np.zeros((size, 7))# 初始化種群的各個粒子的位置 # 用一個 20*2 的矩陣表示種群，每行表示一個粒子 X = np.random.uniform(-5, 5, size=(size, dim))# 初始化種群的各個粒子的速度 V = np.random.uniform(-0.5, 0.5, size=(size, dim))# 初始化粒子歷史最優位置為當當前位置 pbest = X # 計算每個粒子的適應度 for i in range(size):info[i, 3] = calc_f(X[i]) # 目標函數值info[i, 4] = calc_e1(X[i]) # 第一個約束的懲罰項info[i, 5] = calc_e2(X[i]) # 第二個約束的懲罰項# 計算懲罰項的權重，及適應度值 L1, L2 = calc_Lj(info[i, 4], info[i, 5]) info[:, 2] = info[:, 3] + L1 * info[:, 4] + L2 * info[:, 5] # 適應度值 info[:, 6] = L1 * info[:, 4] + L2 * info[:, 5] # 懲罰項的加權求和# 歷史最優 info[:, 0] = info[:, 2] # 粒子的歷史最優位置對應的適應度值 info[:, 1] = info[:, 6] # 粒子的歷史最優位置對應的懲罰項值# 全局最優 gbest_i = info[:, 0].argmin() # 全局最優對應的粒子編號 gbest = X[gbest_i] # 全局最優粒子的位置 gbest_fitness = info[gbest_i, 0] # 全局最優位置對應的適應度值 gbest_e = info[gbest_i, 1] # 全局最優位置對應的懲罰項# 記錄迭代過程的最優適應度值 fitneess_value_list.append(gbest_fitness) # 接下來開始迭代 for j in range(iter_num):# 更新速度V = velocity_update(V, X, pbest=pbest, gbest=gbest)# 更新位置X = position_update(X, V)# 計算每個粒子的目標函數和約束懲罰項for i in range(size):info[i, 3] = calc_f(X[i]) # 目標函數值info[i, 4] = calc_e1(X[i]) # 第一個約束的懲罰項info[i, 5] = calc_e2(X[i]) # 第二個約束的懲罰項# 計算懲罰項的權重，及適應度值L1, L2 = calc_Lj(info[i, 4], info[i, 5])info[:, 2] = info[:, 3] + L1 * info[:, 4] + L2 * info[:, 5] # 適應度值info[:, 6] = L1 * info[:, 4] + L2 * info[:, 5] # 懲罰項的加權求和# 更新歷史最優位置for i in range(size):pbesti = pbest[i]pbest_fitness = info[i, 0]pbest_e = info[i, 1]xi = X[i]xi_fitness = info[i, 2]xi_e = info[i, 6]# 計算更新個體歷史最優pbesti, pbest_fitness, pbest_e = \update_pbest(pbesti, pbest_fitness, pbest_e, xi, xi_fitness, xi_e)pbest[i] = pbestiinfo[i, 0] = pbest_fitnessinfo[i, 1] = pbest_e# 更新全局最優pbest_fitness = info[:, 2]pbest_e = info[:, 6]gbest, gbest_fitness, gbest_e = \update_gbest(gbest, gbest_fitness, gbest_e, pbest, pbest_fitness, pbest_e)# 記錄當前迭代全局之硬度fitneess_value_list.append(gbest_fitness)# 最后繪制適應度值曲線 print('迭代最優結果是：%.5f' % calc_f(gbest)) print('迭代最優變量是：x=%.5f, y=%.5f' % (gbest[0], gbest[1])) print('迭代約束懲罰項是：', gbest_e)#迭代最優結果是：1.00491 #迭代最優變量是：x=1.00167, y=-0.00226 #迭代約束懲罰項是： 0.0 # 從結果看，有多個不同的解的目標函數值是相同的，多測試幾次就發現了# 繪圖 plt.plot(fitneess_value_list[: 30], color='r') plt.title('迭代過程') plt.show()

作者：電氣 余登武

總結

以上是生活随笔為你收集整理的粒子群算法求解带约束优化问题 源码实现的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。