粒子群优化算法（Particle Swarm Optimization）
粒子群算法的思想源于对鸟/鱼群捕食行为的研究，模拟鸟集群飞行觅食的行为，鸟之间通过集体的协作使群体达到最优目的，是一种基于Swarm Intelligence的优化方法。它没有遗传算法的“交叉”(Crossover) 和“变异”(Mutation) 操作，它通过追随当前搜索到的最优值来寻找全局最优。粒子群算法与其他现代优化方法相比的一个明显特色就是所需要调整的参数很少、简单易行，收敛速度快，已成为现代优化方法领域研究的热点。
- 粒子群算法的基本思想
设想这样一个场景：一群鸟在随机搜索食物。已知在这块区域里只有一块食物；所有的鸟都不知道食物在哪里；但它们能感受到当前的位置离食物还有多远。那么找到食物的最优策略是什么呢？
1. 搜寻目前离食物最近的鸟的周围区域
2. 根据自己飞行的经验判断食物的所在。
PSO正是从这种模型中得到了启发，PSO的基础是信息的社会共享
- 算法介绍
每个寻优的问题解都被想像成一只鸟，称为“粒子”。所有粒子都在一个D维空间进行搜索。
所有的粒子都由一个fitness function 确定适应值以判断目前的位置好坏。
每一个粒子必须赋予记忆功能，能记住所搜寻到的最佳位置。
每一个粒子还有一个速度以决定飞行的距离和方向。这个速度根据它本身的飞行经验以及同伴的飞行经验进行动态调整。
粒子速度更新公式包含三部分： 第一部分为“惯性部分”，即对粒子先前速度的记忆；第二部分为“自我认知”部分，可理解为粒子i当前位置与自己最好位置之间的距离；第三部分为“社会经验”部分，表示粒子间的信息共享与合作，可理解为粒子i当前位置与群体最好位置之间的距离。
- 粒子群算法流程
第1步 在初始化范围内，对粒子群进行随机初始化，包括随机位置和速度
第2步 计算每个粒子的适应值
第3步 更新粒子个体的历史最优位置
第4步 更新粒子群体的历史最优位置
第5步 更新粒子的速度和位置
第6步 若未达到终止条件，则转第2步
粒子群算法流程图如下：
- 计算示例
以Ras函数（Rastrigin's Function）为目标函数，求其在x1,x2∈[-5,5]上的最小值。这个函数对模拟退火、进化计算等算法具有很强的欺骗性，因为它有非常多的局部最小值点和局部最大值点，很容易使算法陷入局部最优，而不能得到全局最优解。如下图所示，该函数只在(0,0)处存在全局最小值0。
# -*- coding: utf-8 -*- import sys reload(sys) sys.setdefaultencoding('utf-8') import numpy as np import matplotlib.pyplot as plt # 目标函数定义 def ras(x): y = 20 + x[0]**2 + x[1]**2 - 10*(np.cos(2*np.pi*x[0])+np.cos(2*np.pi*x[1])) return y # 参数初始化 w = 1.0 c1 = 1.49445 c2 = 1.49445 maxgen = 200 # 进化次数 sizepop = 20 # 种群规模 # 粒子速度和位置的范围 Vmax = 1 Vmin = -1 popmax = 5 popmin = -5 # 产生初始粒子和速度 pop = 5 * np.random.uniform(-1,1,(2,sizepop)) v = np.random.uniform(-1,1,(2,sizepop)) fitness = ras(pop) # 计算适应度 i = np.argmin(fitness) # 找最好的个体 gbest = pop # 记录个体最优位置 zbest = pop[:,i] # 记录群体最优位置 fitnessgbest = fitness # 个体最佳适应度值 fitnesszbest = fitness[i] # 全局最佳适应度值 # 迭代寻优 t = 0 record = np.zeros(maxgen) while t < maxgen: # 速度更新 v = w * v + c1 * np.random.random() * (gbest - pop) + c2 * np.random.random() * (zbest.reshape(2,1) - pop) v[v > Vmax] = Vmax # 限制速度 v[v < Vmin] = Vmin # 位置更新 pop = pop + 0.5 * v; pop[pop > popmax] = popmax # 限制位置 pop[pop < popmin] = popmin ''' # 自适应变异 p = np.random.random() # 随机生成一个0~1内的数 if p > 0.8: # 如果这个数落在变异概率区间内，则进行变异处理 k = np.random.randint(0,2) # 在[0,2)之间随机选一个整数 pop[:,k] = np.random.random() # 在选定的位置进行变异 ''' # 计算适应度值 fitness = ras(pop) # 个体最优位置更新 index = fitness < fitnessgbest fitnessgbest[index] = fitness[index] gbest[:,index] = pop[:,index] # 群体最优更新 j = np.argmin(fitness) if fitness[j] < fitnesszbest: zbest = pop[:,j] fitnesszbest = fitness[j] record[t] = fitnesszbest # 记录群体最优位置的变化 t = t + 1 # 结果分析 print zbest plt.plot(record,'b-') plt.xlabel('generation') plt.ylabel('fitness') plt.title('fitness curve') plt.show()
输出结果为如下，可以看到它陷入了局部最优解
[-9.95932740e-01 2.82240596e-04]
删除自适应变异部分的注释，运行后结果如下，可以看出收敛到全局最优解。
[0.00081551 -0.00011209]
- 参数选择与优化
参数w,c1,c2的选择分别关系粒子速度的3个部分：惯性部分、社会部分和自身部分在搜索中的作用。如何选择、优化和调整参数，使得算法既能避免早熟又能比较快的收敛，对工程实践有着重要意义。
1. 惯性权重w描述粒子上一代速度对当前代速度的影响。w值较大，全局寻优能力强，局部寻优能力弱；反之，则局部寻优能力强。当问题空间较大时，为了在搜索速度和搜索精度之间达到平衡，通常做法是使算法在前期有较高的全局搜索能力以得到合适的种子，而在后期有较高的局部搜索能力以提高收敛精度。所以w不宜为一个固定的常数。
wmax最大惯性权重，wmin最小惯性权重，run当前迭代次数，runmax为算法迭代总次数。较大的w有较好的全局收敛能力，较小的w则有较强的局部收敛能力。因此，随着迭代次数的增加，惯性权重w应不断减少，从而使得粒子群算法在初期具有较强的全局收敛能力，而晚期具有较强的局部收敛能力。
2. 学习因子c2=0称为自我认识型粒子群算法，即“只有自我，没有社会”，完全没有信息的社会共享，导致算法收敛速度缓慢；学习因子c1=0称为无私型粒子群算法，即“只有社会，没有自我”，会迅速丧失群体多样性，容易陷入局部最优解而无法跳出；c1,c2都不为0，称为完全型粒子群算法，完全型粒子群算法更容易保持收敛速度和搜索效果的均衡，是较好的选择。
3. 群体大小m是一个整数，m很小时陷入局部最优解的可能性很大；m很大时PSO的优化能力很好，但是当群体数目增长至一定水平时，再增长将不再有显著作用，而且数目越大计算量也越大。群体规模m 一般取20～40，对较难或特定类别的问题 可以取到100～200。
4. 粒子群的最大速度Vmax对维护算法的探索能力与开发能力的平衡很重要，Vmax较大时，探索能力强，但粒子容易飞过最优解；Vmax较小时，开发能力强，但是容易陷入局部最优解。Vmax一般设为每维变量变化范围的10%-20%
参考：
