粒子群算法python特征选择 粒子群算法python代码

简介

粒子群算法(Particle Swarm optimization，简称PSO)是由Eberhart博士和kennedy博士发明的一种启发式算法，是通过模拟鸟群觅食行为而发展起来的一种基于群体协作的随机搜索算法。通常认为它是群集智能 (Swarm intelligence, SI) 的一种。

算法思想

通过模拟鸟群的捕食过程，把每只鸟看成PSO算法中的一个粒子，也就是我们需要求解问题的可行解。整个种群中的鸟在搜寻食物的过程中，会根据条件不断地改变自己的位置和速度。这个条件就是根据个体(可以看成一只叫“我”的鸟)的历史最优p_best和整个种群中全局最优g_best来进行改变下一时刻的位置。可以想象成，每只鸟会一定程度上追随该时刻找食物最厉害的那只鸟所在的位置，这里的位置就是指目标函数的极值点位置。

核心算子

# 伪代码
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)
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的类，这个类包括算法的所有操作方法：
class PSO:
def __init__(self, parameters):
pass
def fitness(self, ind_var):
pass
def update_operator(self, pop_size):
pass
def main(self):
pass
1
2
3
4
5
6
7
8
9
10
11
12
classPSO:
def__init__(self,parameters):
pass
deffitness(self,ind_var):
pass
defupdate_operator(self,pop_size):
pass
defmain(self):
pass
使用__init__()方法初始化参数，包括自变量可取的最大值，最小值，种群大小，迭代代数；使用fitness()方法作为适应度函数评估该个体的函数值，在这里就是函数y的值；使用update_operator()方法根据位置和速度更新粒子下一时刻的位置，挑选出当前代该粒子和种群中所有粒子的最好位置作为历史记录；使用main()方法实现整个算法的循环。
__init__()方法
__init__()方法的代码如下：
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
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
foriinrange(self.pop_size):
forjinrange(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])
iffit>temp:
self.g_best=self.p_best[i]
temp=fit
初始化方法接受传入的参数，包括最大值，最小值，种群大小和迭代代数。通过这些参数初始化产生一个粒子种群，并记录个体最优和全局最优位置。
fitness()方法
fitness()方法用来评估该位置的函数值，代码如下：
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
1
2
3
4
5
6
7
deffitness(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
returny
update_operator()方法
update_operator()方法里包含核心算子，根据速度和位置更新下一时刻的速度，代码如下：
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]
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
defupdate_operator(self,pop_size):
c1=2# 学习因子，一般为2
c2=2
w=0.4# 自身权重因子
foriinrange(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]
# 越界保护
forjinrange(self.var_num):
ifself.pop_x[i][j]
self.pop_x[i][j]=self.bound[0][j]
ifself.pop_x[i][j]>self.bound[1][j]:
self.pop_x[i][j]=self.bound[1][j]
# 更新p_best和g_best
ifself.fitness(self.pop_x[i])>self.fitness(self.p_best[i]):
self.p_best[i]=self.pop_x[i]
ifself.fitness(self.pop_x[i])>self.fitness(self.g_best):
self.g_best=self.pop_x[i]
需要注意的是，在更新位置的过程中有可能超过了我们定义的定义域范围，因此需要进行越界保护，强行将位置拉进定义域内，最后更新个体最好位置和全局最好位置。
main()方法
PSO算法所有的轮子都写好后，我们接下来将它们整合到流程中。代码实现如下：
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
defmain(self):
popobj=[]
self.ng_best=np.zeros((1,self.var_num))[0]
forgeninrange(self.NGEN):
self.update_operator(self.pop_size)
popobj.append(self.fitness(self.g_best))
print('############ Generation {} ############'.format(str(gen+1)))
ifself.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=[tfortinrange(self.NGEN)]
plt.plot(t,popobj,color='b',linewidth=2)
plt.show()
传入循环迭代指定代数，输出目前找到的最好的位置和函数值，并画图得到整个迭代过程。
完整代码
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()
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
importrandom
importnumpyasnp
importmatplotlib.pyplotasplt
classPSO:
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
foriinrange(self.pop_size):
forjinrange(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])
iffit>temp:
self.g_best=self.p_best[i]
temp=fit
deffitness(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
returny
defupdate_operator(self,pop_size):
"""
更新算子：更新下一时刻的位置和速度
"""
c1=2# 学习因子，一般为2
c2=2
w=0.4# 自身权重因子
foriinrange(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]
# 越界保护
forjinrange(self.var_num):
ifself.pop_x[i][j]
self.pop_x[i][j]=self.bound[0][j]
ifself.pop_x[i][j]>self.bound[1][j]:
self.pop_x[i][j]=self.bound[1][j]
# 更新p_best和g_best
ifself.fitness(self.pop_x[i])>self.fitness(self.p_best[i]):
self.p_best[i]=self.pop_x[i]
ifself.fitness(self.pop_x[i])>self.fitness(self.g_best):
self.g_best=self.pop_x[i]
defmain(self):
popobj=[]
self.ng_best=np.zeros((1,self.var_num))[0]
forgeninrange(self.NGEN):
self.update_operator(self.pop_size)
popobj.append(self.fitness(self.g_best))
print('############ Generation {} ############'.format(str(gen+1)))
ifself.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=[tfortinrange(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算法相比遗传算法实现会简单一点，需要设置的参数也少，它的核心就是根据算子更新个体历史最优和全局最优值。寻优效果上比遗传算法更快，但是也不能保证一定能找到全局极值点。本文的代码能够快速运用到其他问题中，可供借鉴和参考。