粒子群演算法解決非線性約束問題
阿新 • • 發佈：2021-08-09
粒子群演算法解決非線性問題
引入
上次我們介紹了粒子群演算法的各種改進，以及matlab軟體自帶的更強大的粒子群演算法，解決的問題都是連續的，無約束的；那麼我們能解決有約束的，非線性問題嗎？
當然可以，不過在此之前，我們需要搞清實現的思路。
解決非線性問題的兩種思路
- 直接在更新新的個體位置之前加入約束條件，滿足約束條件更新個體位置，不滿足約束條件，不更新個體位置
- 加入懲罰項，引入罰函式。在適應度函式上構造懲罰項，對違反約束的情況進行懲罰，將有約束的優化問題轉化為無約束的優化問題。
罰函式解決非線性約束問題理論
在解決非線性約束問題中最常見的手段是引入罰函式。
什麼是罰函式？
罰函式是指在求解最優化問題(無線性約束優化及非線性約束優化)時，在原有目標函式中加上一個障礙函式，而得到一個增廣目標函式，罰函式的功能是對非可行點或企圖穿越邊界而逃離可行域的點賦予一個極大的值，即將有約束最優化問題轉化為求解無約束最優化問題。
簡單來說，就是對對每一個粒子進行判斷，判斷該粒子的取值是否滿足約束條件，如果滿足，則進行正常的適應度計算，如果不滿足，則將該粒子的適應度直接賦值為一個很大的數。這個很大的數我們叫做懲罰。
約束條件一般由等式約束和不等式約束組成。我們一般將等式約束轉化為不等式約束處理。
加入懲罰項之後，我們新的目標函式變為：
其中\(σ\)是懲罰因子，一般取$σ=t√t $，P(x)是整體懲罰項，P(x)的計算方法如下：
1.對於不等式\(g_i(x)<=0\),懲罰項
即如果\(g_i(x)<=0\),則\(e_i(x)=0\),否則\(e_i(x)=-g_i(x)\)
2.對於等式約束,一個簡單的方法設定一個等式約束容忍度值\(ε\)
因此等式約束的懲罰項為:
整體懲罰性P(x)是各個約束懲罰項的和，即：
由於約束條件之間的差異，某些約束條件可能對個體違反程度起著決定性的作用，為了平衡這種差異，對懲罰項做歸一化處理:
最後得到的粒子xi的整體懲罰項P(x)的計算公式：
在粒子群演算法中，每一步迭代都會更新Pbest和Gbest，雖然可以將有約束問題轉換為無約束問題進行迭代求解，但是問題的解xi依然存在不滿足約束條件的情況，因此需要編制一些規則來比較兩個粒子的優劣，規則如下：
- 如果兩個粒子\(xi\)和\(xj\)都可行，則比較其適應度函式\(f(x_i)\)和\(f(x_j)\)，值小的粒子為優。
- 當兩個粒子\(xi\)和\(xj\)都不可行，則比較懲罰項\(P(x_i)\)和\(P(x_j）\)，違背約束程度小的粒子更優。
- 當粒子\(x_i\)可行而粒子\(x_j\)不可行，選可行解。
對於粒子的上下限約束 可以體現在位置更新函式裡，不必加懲罰項。 具體思路就是遍歷每一個粒子的位置，如果超除上下限，位置則更改為上下限中的任何一個位置.
罰函式解決非線性約束問題實現
例子
初始化引數
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 = [] # 記錄每次迭代過程中的種群適應度值變化
計算適應度函式和計算約束懲罰項函式
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
定義粒子群演算法的速度更新函式，位置更新函式
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
每個粒子歷史最優位置更優函式，以及整個群體歷史最優位置更新函式
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: 上一次迭代的全域性最優位置
: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
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] = pbesti
info[i, 0] = pbest_fitness
info[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()
本文來自：粒子群演算法求解帶約束優化問題 原始碼實現_總裁餘（餘登武）部落格-CSDN部落格_帶約束粒子群
可參考的其他相關部落格：