粒子群算法

粒子群算法也称粒子群优化算法PSO（Particle Swarm Optimization），属于进化算法EA（Evolutionary Algorithm）的一种，实现容易、精度高、收敛快，是一种并行算法。它是从随机解出发，通过迭代寻找最优解，也是通过适应度来评价解的品质，它通过追随当前搜索到的最优值来寻找全局最优。

算法原理

设x为粒子起始位置，v为粒子飞行速度，p为搜索到的粒子的最优位置。

PSO初始化为一群随机粒子，然后通过迭代找到最优解。在每一次迭代中，粒子通过跟踪两个极值来更新自己：第一个就是粒子本身所找到的最优解，这个解称为个体极值pbest；另一个极值是整个种群目前找到的最优解，这个极值是全局极值gbest。粒子始终跟随这两个极值并在每一次迭代中更新自己的位置和速度，直到找到最优解。

具体如下：

假设在一个D维的目标搜索空间中，有N个粒子组成一个群落，其中第i个粒子为一个D维向量：

第i个粒子的飞行速度也是一个D维向量：

第i个粒子迄今为止搜索到的最优位置，即个体极值Pbest：

整个粒子群迄今为止搜索到的最优位置，即全局极值Gbest：

在找到这两个最优值时，粒子根据如下公式来更新自己的速度和位置：

其中，c1、c2为学习因子，也称加速常数（acceleration constant）；r1、r2为[0,1]范围内的均匀随机数。 第一个公式右侧包括三部分，

• 第1部分为惯性inertia或动量momentum，代表了粒子有维持自己先前速度的趋势；

• 第2部分为认知cognition，反映了粒子对自己历史经验的记忆或回忆，代表粒子有向自身历史最佳位置逼近的趋势；

• 第3部分为社会social，反映了粒子间协同合作与知识共享的群体历史经验，代表粒子有向群体或邻近历史最佳位置逼近的趋势。

算法流程

1. 初始化粒子群，包括群体规模N、每个粒子的位置xi和速度vi计算每个粒子的适应度值Fit[i]
2. 对每个粒子，用它的适应度值Fit[i]和个体极值pbest(i)比较，如果Fit[i]>pbest(i)，则用Fit[i]替换掉 pbest(i)
3. 对每个粒子，用其适应度Fit[i]和全局极值gbest(i)比较，如果Fit[i]>gbest(i)，则用Fit[i]替换掉gbest(i)
4. 根据公式更新粒子的位置xi和速度vi
5. 如果满足结束条件，误差足够好或达到最大循环次数，退出，否则返回2

代码实现

import numpy as np
import matplotlib.pyplot as plt

'''
******************************目标函数***********************************
'''
#  目标函数 ， 同时也将他设定为 适应度函数   （一维）
def get_fitness(x):
    y = -x**2 - 8*x -1
    return y 


 
'''
******************************初始化***********************************
'''

# ************************初始化简单参数*************************************
w = 0.8    			#惯性因子
c1 = 2     			#自身认知因子
c2 = 2     			#社会认知因子
r1 = 0.6  			#自身认知学习率
r2 = 0.3  			#社会认知学习率
pN = 100             #粒子数量
dim = 1             #搜索维度  
max_iter = 1000    	#最大迭代次数


X = np.zeros((pN, dim))       #初始粒子的位置 
V = np.zeros((pN, dim))       #初始粒子的速度 
p_best = np.zeros((pN, dim), dtype = float)   #单个粒子的历史最佳位置
g_best = np.zeros((1, dim), dtype = float)    #全局最佳位置  
p_bestfit = np.zeros(pN)      #单个粒子的历史最佳适应值  
g_bestfit = -1e15             #全局最佳适应值  





# *********************初始化个体和全局的最优解********************************
for i in range(pN):  
    #初始化每一个粒子的位置和速度
    X[i] = np.random.uniform(0,5,[1, dim])  
    V[i] = np.random.uniform(0,5,[1, dim]) 

    # 得到初始时的每一个粒子的历史最佳位置（暂且将第一次给的值定为最佳位置）
    p_best[i] = X[i] 
    p_bestfit[i] = get_fitness(X[i])  #得到对应的fit值
    

    # 得到初始时的全局最佳位置和相应的适应度值（每遍历一个粒子就更新一次）
    # （这里不直接用max函数求最大值是因为不一定下一个粒子比前一个粒子更优，从而导致冗余计算。下面的那个for循环同理）
  
g_bestfit = max(p_bestfit)
g_best = max(p_best)  	#得到全局最佳位置





'''
***************************迭代***********************************#
'''
fitness = []  #记录每次迭代的最优解，画图时用

for _ in range(max_iter):  

    #  更新每个粒子的位置和速度
    for i in range(pN):  	
 
        # 更新g_best 和 p_best  
        temp = get_fitness(X[i])  #获得当前位置的适应值

        if temp > p_bestfit[i]:      #先更新个体最优  
            p_bestfit[i] = temp  
            p_best[i] = X[i]  
            
            if p_bestfit[i] > g_bestfit:  #再更新全局最优  （注意：这里不直接用max函数求最大值是因为最外面还有一层max_iter迭代循环，从而更容易导致冗余计算；而初始化时的那重循环外面再无循环层，所以即使使用max函数也不会造成什么冗余。）
                g_best = X[i]  
                g_bestfit = p_bestfit[i]  

        V[i] = w*V[i] + c1*r1*(p_best[i] - X[i]) + c2*r2*(g_best - X[i])  
        X[i] = X[i] + V[i] 
    '''
    #  这样写没错，可以代替上面的那层if语句，但不能这样写，因为最外面还有一重max_iter循环，如果下一次迭代不比当前最优解还好的话，就造成了max函数计算冗余
    g_bestfit = max(p_bestfit)
    g_best = max(p_best)  
    '''
	
    fitness.append(g_bestfit)  # 记录每次迭代后找到的全局最优解
    
print(g_best, g_bestfit)






'''
***************************画图***********************************#
'''
plt.figure(1)
plt.title("Figure1")
plt.xlabel("iterators", size=14)
plt.ylabel("fitness", size=14)
t = np.array([t for t in range(0, max_iter)])
fitness = np.array(fitness)
plt.plot(t, fitness, color='b', linewidth=3)
plt.show()

运行结果

很明显对于函数y = -x**2 - 8*x -1，当x=-4时取最大值y=15