简介
粒子群算法(Particle Swarm optimization,简称PSO)是由Eberhart博士和kennedy博士发明的一种启发式算法,是通过模拟鸟群觅食行为而发展起来的一种基于群体协作的随机搜索算法。通常认为它是群集智能 (Swarm intelligence, SI) 的一种。
算法思想
通过模拟鸟群的捕食过程,把每只鸟看成PSO算法中的一个粒子,也就是我们需要求解问题的可行解。整个种群中的鸟在搜寻食物的过程中,会根据条件不断地改变自己的位置和速度。这个条件就是根据个体(可以看成一只叫“我”的鸟)的历史最优p_best
和整个种群中全局最优g_best
来进行改变下一时刻的位置。可以想象成,每只鸟会一定程度上追随该时刻找食物最厉害的那只鸟所在的位置,这里的位置就是指目标函数的极值点位置。
核心算子
1 2 3 |
# 伪代码 v[i](t+1) = w * v[i](t) + c1 * rand() * (pbest[i] - present[i](t)) + c2 * rand() * (gbest - present[i](t)) present[i](t+1) = present[i](t) + v[i](t+1) |
v[i](t+1)
表示一个i
粒子下一时刻的速度,它是通过计算该粒子历史最优值pbest[i]
和全部粒子全局最优值gbest
与该时刻当前位置present[i](t)
的差值得到的。其中w
是惯性权值,c1
和c2
表示学习参数。有了下一时刻的速度后,就能得到下一时刻的位置present[i](t+1)
。
示例
我们还是以一个简单的例子来介绍算法的流程。假设我们需要寻找函数y=x12+x22+x33+x44
在[1,30]
之间的最大值。我们很容易就知道,当x1=x2=x3=x4=30
时,该函数能取到最大值。
构建一个叫PSO
的类,这个类包括算法的所有操作方法:
1 2 3 4 5 6 7 8 9 10 11 12 |
class PSO: def __init__(self, parameters): pass def fitness(self, ind_var): pass def update_operator(self, pop_size): pass def main(self): pass |
使用__init__()
方法初始化参数,包括自变量可取的最大值,最小值,种群大小,迭代代数;使用fitness()
方法作为适应度函数评估该个体的函数值,在这里就是函数y
的值;使用update_operator()
方法根据位置和速度更新粒子下一时刻的位置,挑选出当前代该粒子和种群中所有粒子的最好位置作为历史记录;使用main()
方法实现整个算法的循环。
__init__()方法
__init__()
方法的代码如下:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 |
def __init__(self, parameters): # 初始化 self.NGEN = parameters[0] # 迭代的代数 self.pop_size = parameters[1] # 种群大小 self.var_num = len(parameters[2]) # 变量个数 self.bound = [] # 变量的约束范围 self.bound.append(parameters[2]) self.bound.append(parameters[3]) self.pop_x = np.zeros((self.pop_size, self.var_num)) # 所有粒子的位置 self.pop_v = np.zeros((self.pop_size, self.var_num)) # 所有粒子的速度 self.p_best = np.zeros((self.pop_size, self.var_num)) # 每个粒子最优的位置 self.g_best = np.zeros((1, self.var_num)) # 全局最优的位置 # 初始化第0代初始全局最优解 temp = -1 for i in range(self.pop_size): for j in range(self.var_num): self.pop_x[i][j] = random.uniform(self.bound[0][j], self.bound[1][j]) self.pop_v[i][j] = random.uniform(0, 1) self.p_best[i] = self.pop_x[i] # 储存最优的个体 fit = self.fitness(self.p_best[i]) if fit > temp: self.g_best = self.p_best[i] temp = fit |
初始化方法接受传入的参数,包括最大值,最小值,种群大小和迭代代数。通过这些参数初始化产生一个粒子种群,并记录个体最优和全局最优位置。
fitness()方法
fitness()
方法用来评估该位置的函数值,代码如下:
1 2 3 4 5 6 7 |
def fitness(self, ind_var): x1 = ind_var[0] x2 = ind_var[1] x3 = ind_var[2] x4 = ind_var[3] y = x1**2 + x2**2 + x3**3 + x4**4 return y |
update_operator()方法
update_operator()
方法里包含核心算子,根据速度和位置更新下一时刻的速度,代码如下:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 |
def update_operator(self, pop_size): c1 = 2 # 学习因子,一般为2 c2 = 2 w = 0.4 # 自身权重因子 for i in range(pop_size): # 更新速度 self.pop_v[i] = w * self.pop_v[i] + c1 * random.uniform(0, 1) * ( self.p_best[i] - self.pop_x[i]) + c2 * random.uniform(0, 1) * (self.g_best - self.pop_x[i]) # 更新位置 self.pop_x[i] = self.pop_x[i] + self.pop_v[i] # 越界保护 for j in range(self.var_num): if self.pop_x[i][j] < self.bound[0][j]: self.pop_x[i][j] = self.bound[0][j] if self.pop_x[i][j] > self.bound[1][j]: self.pop_x[i][j] = self.bound[1][j] # 更新p_best和g_best if self.fitness(self.pop_x[i]) > self.fitness(self.p_best[i]): self.p_best[i] = self.pop_x[i] if self.fitness(self.pop_x[i]) > self.fitness(self.g_best): self.g_best = self.pop_x[i] |
需要注意的是,在更新位置的过程中有可能超过了我们定义的定义域范围,因此需要进行越界保护,强行将位置拉进定义域内,最后更新个体最好位置和全局最好位置。
main()方法
PSO算法所有的轮子都写好后,我们接下来将它们整合到流程中。代码实现如下:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 |
def main(self): popobj = [] self.ng_best = np.zeros((1, self.var_num))[0] for gen in range(self.NGEN): self.update_operator(self.pop_size) popobj.append(self.fitness(self.g_best)) print('############ Generation {} ############'.format(str(gen + 1))) if self.fitness(self.g_best) > self.fitness(self.ng_best): self.ng_best = self.g_best print('最好的位置:{}'.format(self.ng_best)) print('最大的函数值:{}'.format(self.fitness(self.ng_best))) print("---- End of (successful) Searching ----") plt.figure() plt.title("Figure1") plt.xlabel("iterators", size=14) plt.ylabel("fitness", size=14) t = [t for t in range(self.NGEN)] plt.plot(t, popobj, color='b', linewidth=2) plt.show() |
传入循环迭代指定代数,输出目前找到的最好的位置和函数值,并画图得到整个迭代过程。
完整代码
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 |
import random import numpy as np import matplotlib.pyplot as plt class PSO: def __init__(self, parameters): """ particle swarm optimization parameter: a list type, like [NGEN, pop_size, var_num_min, var_num_max] """ # 初始化 self.NGEN = parameters[0] # 迭代的代数 self.pop_size = parameters[1] # 种群大小 self.var_num = len(parameters[2]) # 变量个数 self.bound = [] # 变量的约束范围 self.bound.append(parameters[2]) self.bound.append(parameters[3]) self.pop_x = np.zeros((self.pop_size, self.var_num)) # 所有粒子的位置 self.pop_v = np.zeros((self.pop_size, self.var_num)) # 所有粒子的速度 self.p_best = np.zeros((self.pop_size, self.var_num)) # 每个粒子最优的位置 self.g_best = np.zeros((1, self.var_num)) # 全局最优的位置 # 初始化第0代初始全局最优解 temp = -1 for i in range(self.pop_size): for j in range(self.var_num): self.pop_x[i][j] = random.uniform(self.bound[0][j], self.bound[1][j]) self.pop_v[i][j] = random.uniform(0, 1) self.p_best[i] = self.pop_x[i] # 储存最优的个体 fit = self.fitness(self.p_best[i]) if fit > temp: self.g_best = self.p_best[i] temp = fit def fitness(self, ind_var): """ 个体适应值计算 """ x1 = ind_var[0] x2 = ind_var[1] x3 = ind_var[2] x4 = ind_var[3] y = x1 ** 2 + x2 ** 2 + x3 ** 3 + x4 ** 4 return y def update_operator(self, pop_size): """ 更新算子:更新下一时刻的位置和速度 """ c1 = 2 # 学习因子,一般为2 c2 = 2 w = 0.4 # 自身权重因子 for i in range(pop_size): # 更新速度 self.pop_v[i] = w * self.pop_v[i] + c1 * random.uniform(0, 1) * ( self.p_best[i] - self.pop_x[i]) + c2 * random.uniform(0, 1) * (self.g_best - self.pop_x[i]) # 更新位置 self.pop_x[i] = self.pop_x[i] + self.pop_v[i] # 越界保护 for j in range(self.var_num): if self.pop_x[i][j] < self.bound[0][j]: self.pop_x[i][j] = self.bound[0][j] if self.pop_x[i][j] > self.bound[1][j]: self.pop_x[i][j] = self.bound[1][j] # 更新p_best和g_best if self.fitness(self.pop_x[i]) > self.fitness(self.p_best[i]): self.p_best[i] = self.pop_x[i] if self.fitness(self.pop_x[i]) > self.fitness(self.g_best): self.g_best = self.pop_x[i] def main(self): popobj = [] self.ng_best = np.zeros((1, self.var_num))[0] for gen in range(self.NGEN): self.update_operator(self.pop_size) popobj.append(self.fitness(self.g_best)) print('############ Generation {} ############'.format(str(gen + 1))) if self.fitness(self.g_best) > self.fitness(self.ng_best): self.ng_best = self.g_best.copy() print('最好的位置:{}'.format(self.ng_best)) print('最大的函数值:{}'.format(self.fitness(self.ng_best))) print("---- End of (successful) Searching ----") plt.figure() plt.title("Figure1") plt.xlabel("iterators", size=14) plt.ylabel("fitness", size=14) t = [t for t in range(self.NGEN)] plt.plot(t, popobj, color='b', linewidth=2) plt.show() if __name__ == '__main__': NGEN = 100 popsize = 100 low = [1, 1, 1, 1] up = [30, 30, 30, 30] parameters = [NGEN, popsize, low, up] pso = PSO(parameters) pso.main() |
在if __name__ == "__main__":
语句后面,我们设定所有的参数。总共跑NGEN=100代,每代的种群大小为100。
最后得到以下结果:
迭代过程如下:
总结
粒子群PSO算法相比遗传算法实现会简单一点,需要设置的参数也少,它的核心就是根据算子更新个体历史最优和全局最优值。寻优效果上比遗传算法更快,但是也不能保证一定能找到全局极值点。本文的代码能够快速运用到其他问题中,可供借鉴和参考。
请问是不是有bug,按道理不需要self.ng_best这个参数,因为g_best根据代码逻辑会获取最优值,但是我实验一下,一步一步对比,如果没有self.ng_best参数,在运行过程中会把pop_x[i]无理由赋值到g_best,导致g_best变为不是最优解。