import matplotlib;
matplotlib.use('TkAgg')
from pylab import *
mpl.rcParams['font.sans-serif'] = ['SimHei'] # 指定默认字体
mpl.rcParams['axes.unicode_minus'] = False # 解决保存图像是负号'-'显示为方块的问题
import numpy as np
from tqdm import tqdm # 进度条设置
import random
# 更改了位置更新函数，
class PBO(object):
""""
M:种群最大容量
m：种群当前规模大小
N：目标函数解的维度
T：迭代次数
V：最大视野距离
P：初始化种群数量百分比
K：动态调整规模参数
B：变量限制范围
X：解集合
"""
def __init__(self):
self.Weight = [95, 75, 23, 73, 50, 51, 60, 57, 89, 98,
95, 75, 23, 73, 50, 51, 60, 57, 89, 98,
95, 75, 23, 73, 50, 51, 60, 57, 89, 98,
95, 75, 23, 73, 50, 51, 60, 57, 89, 98,
95, 75, 23, 73, 50, 51, 60, 57, 89, 98,
95, 75, 23, 73, 50, 51, 60, 57, 89, 98,
95, 75, 23, 73, 50, 51, 60, 57, 89, 98,
95, 75, 23, 73, 50, 51, 60, 57, 89, 98,
] # 物品体积
self.Value = [89, 59, 19, 43, 100, 101, 44, 16, 7, 64,
89, 59, 19, 43, 99, 101, 44, 16, 7, 64,
89, 59, 19, 43, 98, 101, 44, 16, 7, 64,
89, 59, 19, 43, 97, 101, 44, 16, 7, 64,
89, 59, 19, 43, 96, 101, 44, 16, 7, 64,
89, 59, 19, 43, 95, 101, 44, 16, 7, 64,
89, 59, 19, 43, 77, 101, 44, 16, 7, 64,
89, 59, 19, 43, 79, 101, 44, 16, 7, 64,
] # 物品价值
self.M = 100
self.N = len(self.Weight)
self.T = 1000
self.P = 0.75
self.K = 0.25
self.m = int(np.floor(self.M * self.P)) # 随后期的繁衍和死亡会变化
self.V = 0.3
# 二进制参数
self.a = 0.1
self.w = 0
self.b = 1 # 管控局部搜索
# 初始化种群
def init_x(self):
"""
:return:
"""
X = np.zeros((self.M, self.N))
X[:self.m] = np.random.choice([0,1], size=(self.m, self.N))
return X
def init_V(self,X):
V = np.zeros((self.M, self.N))
V[:self.m] = np.random.random((self.m, self.N))
V = V
return V
def fitness(self, X):
"""
有待更改
:param X:
:return:
"""
y = np.sum(X * self.Value)
w = np.sum(X * self.Weight)
if w > 407:
y = 0
return y
def up_data(self, index, angle, radius, action, X):
"""
局部搜索中的位置更新公式 公式看原版文献csdn有误
:param index:每一只熊的的编号
:param angle:angles [0,2pi]
:param X[index][0]:第index只熊的0维
:param radius:
:param action:
:param X:
:return:
"""
temp = np.zeros(self.N)
result = 0
if "add".__eq__(action):
x_temp = X[index][0] + radius * np.cos(angle[0])
else:
x_temp = X[index][0] - radius * np.cos(angle[0])
temp[0] = x_temp
# 因为从第2到n-1只熊由通用计算工时
for i in range(self.N - 2):
result += np.sin(angle[i])
if "add".__eq__(action):
x_temp = X[index][i + 1] + radius * (result + np.cos(angle[i + 1]))
else:
x_temp = X[index][i + 1] - radius * (result + np.cos(angle[i + 1]))
temp[i + 1] = x_temp
if "add".__eq__(action):
x_temp = X[index][self.N - 1] + radius * (result + np.sin(angle[self.N - 2]))
else:
x_temp = X[index][self.N - 1] - radius * (result + np.sin(angle[self.N - 2]))
temp[self.N - 1] = x_temp
return temp
def move_local(self, X, xbest, V):
"""
局部更新 （包含搜索半径 和局部更新）
:param X:
:return:
r=4acos(φ)sin(φ)
已调好
"""
# 对每个北极熊设计随机角度 根据公式计算出搜索半径 r=4acos(φ)sin(φ)
V1 = V.copy()
V2 = V.copy()
V3 = V.copy()
for i in range(self.m):
angles = np.random.random(self.N - 1) * 2 * np.pi # [0,2pi]用于局部搜索中的位置更新
sta = np.random.random() * self.V # [0,0.3]
angle0 = np.random.random() * np.pi / 2 # [0,pi/2]
# 搜索半径
r = 4 * sta * np.cos(angle0) * np.sin(angle0)
tmp = self.up_data(i, angles, r, 'add', X)
temp = tmp.copy()
V1[i] = self.a * (V[i] - temp)
temp = abs((np.exp(V1[i]) - np.exp(-1 * V1[i])) / (np.exp(V1[i]) + np.exp(-1 * V1[i])))
for j in range(self.N):
ran = np.random.random()
if ran >= temp[j]:
tmp[j] = X[i][j]
else:
tmp[j] = 1 - X[i][j]
if self.fitness(tmp) > self.fitness(X[i]):
X[i] = tmp
V3[i] = V1[i]
else:
tmp = self.up_data(i, angles, r, 'minus', X)
temp = tmp.copy()
V2[i] = self.a * (xbest - temp)
temp = abs((np.exp(V2[i]) - np.exp(-1 * V2[i])) / (np.exp(V2[i]) + np.exp(-1 * V2[i])))
for j in range(self.N):
ran = np.random.random()
# ran = 0.6
if ran >= temp[j]:
tmp[j] = X[i][j]
else:
tmp[j] = 1 - X[i][j]
if self.fitness(tmp) > self.fitness(X[i]):
X[i] = tmp
V3[i] = V2[i]
return X, V3
def move_globle(self, X, xbest, V):
"""
浮水漂移，全局搜索
X：M行N列的二维矩阵 ，每一个子矩阵代表一个多纬度的北极熊
全局位置更新
:return:
"""
V1 = V.copy()
for i in range(self.m):
# 每一只熊的W
w = np.sqrt(np.sum((xbest - X[i]) ** 2)) # w是n-1个维度的欧氏距离
alpha = np.random.random()
gamma = np.random.random() * w # 0-w之间的随机数
x_temp = X[i] + np.sign(w) * alpha + gamma
x_temp_copy = x_temp.copy()
V1[i] = self.a * (xbest - x_temp_copy)
x_temp_copy = abs((np.exp(V1[i]) - np.exp(-1 * V1[i])) / (np.exp(V1[i]) + np.exp(-1 * V1[i])))
for j in range(self.N):
ran = np.random.random()
if ran >= x_temp_copy[j]:
x_temp[j] = X[i][j]
else:
x_temp[j] = 1 - X[i][j]
if self.fitness(x_temp) > self.fitness(X[i]):
X[i] = x_temp
V[i] = V1[i]
return X, V
def fit_best(self, X, f_best, x_best):
fits = np.zeros(self.m)
for i in range(self.m):
# 找个体最优
fits[i] = self.fitness(X[i])
arg = np.argsort(fits[:self.m])
if fits[arg[-1]] >= f_best:
f_best = fits[arg[-1]]
x_best = X[arg[-1]]
return f_best, x_best, fits
def zhongqunfangyan(self, X, x_best, fits):
k = np.random.random()
arg = np.argsort(fits[:self.m])
# 种群繁衍
if self.m < self.M - 1 and k > 0.75:
# 从前10%中随机抽取一个 个体与最好个体交配 生成新的个体
idx = random.randint(1, np.floor(self.m * 0.1) + 1)
X[self.m] = (x_best + X[arg[self.m - 1 - idx]]) / 2
self.m += 1
# 最弱个体死亡
if self.m > self.M * 0.5 and k < 0.25:
# arg中排序最后的为最弱个体
# 死亡之后 后面的所有个体往前移动一个位置
X[arg[0]:self.m] = X[arg[0] + 1:self.m + 1]
self.m -= 1
return X
def main(self):
fitness_list = [] # 画图使用
x_best_list = []
X = self.init_x()
f_best = -1 * np.inf # 表示负无穷
x_best = 0
f_best, x_best, fits = self.fit_best(X, f_best, x_best)
f_bestV, V_best, V_fits = self.fit_best(X, f_best, x_best)
fitness_list.append(f_best)
V = self.init_V(X)
# 先局部后全局搜索
for i in tqdm(range(self.T)):
X, V = self.move_local(X, x_best, V)
X, V = self.move_globle(X, V_best, V)
# print(i)
# #计算当前适应度
f_best, x_best, fits = self.fit_best(X, f_best, x_best)
f_bestV, V_best, V_fits = self.fit_best(X, f_best, x_best)
fitness_list.append(f_best)
# print(x_best)
# 种群动态变化
X = self.zhongqunfangyan(X, x_best, fits)
V = self.zhongqunfangyan(V, V_best, fits)
X = np.trunc(X)
f_best, x_best, fits = self.fit_best(X, f_best, x_best)
print("最大适应度值:", f_best, '个体为：', x_best, "重量", np.sum(x_best * self.Weight))
plt.plot(fitness_list, "--")
plt.xlabel('迭代次数', size=15)
plt.ylabel('fitness', size=15)
plt.legend()
plt.show()
if __name__ == '__main__':
for i in range(20):
pbo = PBO()
pbo.main()
浙公网安备 33010602011771号