粒子群算法（PSO）基本思想
粒子群（PSO）算法最早是由美国电气工程师Eberhart和社会心理学家Kennedy在1995年基于群鸟觅食提出来的。
群鸟觅食其实是一个最佳决策的过程， 与人类决策的过程相似。Boyd和Re chars on探索了人类的决策过程，并提出了个体学习和文化传递的概念。根据他们的研究成果，人们在决策过程中常常会综合两种重要的信息：第一种是他们自己的经验，即他们根据以前自己的尝试和经历，已经积累了一定的经验，知道怎样的状态会比较好；第二种是其他人的经验，即从周围人的行为获取知识，从中知道哪些选择是正面的，哪些选择是消极的。
同样的道理，群鸟在觅食的过程中，每只鸟的初始状态都处于随机位置，且飞翔的方向也是随机的。每只鸟都不知道食物在哪里，但是随着时间的推移，这些初始处于随机位置的鸟类通过群内相互学习、信息共享和个体不断积累字觅食物的经验，自发组织积聚成一个群落，并逐渐朝唯一的目标-—食物前进。每只鸟能够通过一定经验和信息估计目前所处的位置对于能寻找到食物有多大的价值，即多大的适应值；每只鸟能够记住自己所找到的最好位置，称之为局部最优(p best) 。此外， 还能记住群鸟中所有个体所能找到的最好位置， 称为全局最优(g best) ， 整个鸟群的觅食中心都趋向全局最优移动， 这在生物学上称之为“同步效应”。通过鸟群觅食的位置不断移动，即不断迭代，可以使鸟群朝食物步步进逼。
PSO算法的优化
近年来， 一些学者将PSO算法推广到约束优化问题，其关键在于如何处理好约束， 即解的可行性。如果约束处理得不好，其优化的结果往往会出现不能收敛和结果是空集的状况。基于PSO算法的约束优化工作主要分为两类：
- ①罚函数法。罚函数的目的是将约束优化问题转化成无约束优化问题。
- ②将粒子群的搜索范围都限制在条件约束簇内，即在可行解范围内寻优。
Parsopoulos等采用罚函数法，利用非固定多段映射函数对约束优化问题进行转化，再利用PSO算法求解转化后的问题，仿真结果显示PSO算法相对遗传算法更具有优越性，但其罚函数的设计过于复杂，不利于求解；Hu等采用可行解保留策略处理约束，即一方面更新存储中所有粒子时仅保留可行解，另一方面在初始化阶段所有粒子均从可行解空间取值， 然而初始可行解空间对于许多问题是很难确定的，Ray等提出了具有多层信息共享策略的粒子群原理来处理约束， 根据约束矩阵采用多层Pareto排序机制来产生优良粒子、进而用一些优良的粒子来决定其余个体的搜索方向。
PSO算法的优缺点
- 1)不依赖于问题信息，采用实数求解，算法通用性强。
- 需要调整的参数少，原理简单，容易实现，这是PSO算法的最大优点。
- 3)协同搜索，同时利用个体局部信息和群体全局信息指导搜索。
- 收敛速度快， 算法对计算机内存和CPU要求不高。
- 5)更容易飞越局部最优信息。对于目标函数仅能提供极少搜索最优值的信息，在其他算法无法辨别搜索方向的情况下，PSO算法的粒子具有飞越性的特点使其能够跨过搜索平面上信息严重不足的障碍，飞抵全局最优目标值。比如Generalized Rosenbrock函数全局最小值在原占附近.但是此函数全局最优值与可到达的局部最优值之间右一条独长的山公，曲面山谷中点的最速下降方向几乎与到函数最小值的最佳方向垂直，找到全局最小值的可能性微乎其微， 但是PSO算法完全有可能找到全局最优值。
同时， PSO算法的缺点也是显而易见的：
1)算法局部搜索能力较差，搜索精度不够高。
2)算法不能绝对保证搜索到全局最优解，主要有两方面的原因：
①有时粒子群在俯冲过程中会错失全局最优解。粒子飞翔过程中的俯冲动作使搜索行为不够精细，不容易发现全局最优目标值，所以对粒子的最大飞翔速度进行限制既是为了使粒子不要冲出搜索区域的边界，同时也是为了使搜索行为不至于太粗糙。②应用PSO算法处理高维复杂问题时，算法可能会早熟收敛，也就是粒子群在没有找到全局最优信息之前就陷入停顿状态，飞翔的动力不够，粒子群丧失了多样性，各粒子之间的抱合力增强，紧紧地聚集在一起，并且它们的飞翔速度几乎为零，虽然此时粒子距离全局最优解并不远，但是几乎为零的飞翔速度使其跳出停滞不前的状态，各个粒子力不从心。这些停滞不前的早熟点未必都是局部最优点，也可能是位于局部最优点邻域内的其他点，这一点与梯度搜索法不同，梯度搜索法如果出现早熟，通常只会陷人局部最优点，而不可能陷人局部最优点邻域内的其他点，这一点与梯度搜索算法不同，梯度搜索算法如果出现早熟，通常只会陷入局部最优点，而·不可能陷入局部最优点领域内的其它点。3)算法搜索性能对参数具有一定的依赖性。对于特定的优化问题，如果用户经验不足，参数调整的确是个棘手的问题。参数值的大小直接影响到算法是否收敛以及求解结果的精度。4)PSO算法是一种概率算法，算法理论不完善，缺乏独特性，理论成果偏少。从数学角度严格证明算法结果的正确性和可靠性还比较困难；缺少算法结构设计和参数选取的实用性指导原则，特别是全局收敛研究和大型多约束非线性规划的研究成果非常少。
PSO算法程序设计
PSO算法实现的流程图如下图所示：
PSO算法设计的具体步骤如下：
- 步骤1：初始化粒子群（速度和位置）、惯性因子、加速常数、最大迭代次数、算法终止的最小允许误差。
- 步骤2：评价每个粒子的初始适应值。
- 步骤3：将初始适应值作为当前每个粒子的局部最优值，并将各适应值对应的位置作为每个粒子的局部最优值所在的位置。
- 步骤4：将最佳初始适应值作为当前全局最优值，并将最佳适应值对应的位置作为全局最优值所在的位置。
- 步骤5：依据式(8-1)更新每个粒子当前的飞翔速度。
- 步骤6对每个粒子的飞翔速度进行限幅处理，使之不能超过设定的最大飞翔速度。
- 步骤7依据式(8-2)更新每个粒子当前所在的位置。
- 步骤8比较当前每个粒子的适应值是否比历史局部最优值好，如果好，则将当前粒子适应值作为粒子的局部最优值，其对应的位置作为每个粒子的局部最优值所在的位置。
- 步骤9在当前群中找出全局最优值，并将当前全局最优值对应的位置作为粒子群的全局最优值所在的位置。
- 步骤10：重复步骤（5）~（9），直到满足设定的最小误差或最大迭代次数
- 步骤11：输出粒子群的全局最优值和其对应的位置以及每个粒子的局部最优值和其对应的位置。
python简单实现
import numpy as np import matplotlib.pyplot as plt class PSO(object): def __init__(self, population_size, max_steps): self.w = 0.6 # 惯性权重 self.c1 = self.c2 = 2 self.population_size = population_size # 粒子群数量 self.dim = 2 # 搜索空间的维度 self.max_steps = max_steps # 迭代次数 self.x_bound = [-10, 10] # 解空间范围 self.x = np.random.uniform(self.x_bound[0], self.x_bound[1], (self.population_size, self.dim)) # 初始化粒子群位置 self.v = np.random.rand(self.population_size, self.dim) # 初始化粒子群速度 fitness = self.calculate_fitness(self.x) self.p = self.x # 个体的最佳位置 self.pg = self.x[np.argmin(fitness)] # 全局最佳位置 self.individual_best_fitness = fitness # 个体的最优适应度 self.global_best_fitness = np.min(fitness) # 全局最佳适应度 def calculate_fitness(self, x): return np.sum(np.square(x), axis=1) def evolve(self): fig = plt.figure() for step in range(self.max_steps): r1 = np.random.rand(self.population_size, self.dim) r2 = np.random.rand(self.population_size, self.dim) # 更新速度和权重 self.v = self.w*self.v+self.c1*r1*(self.p-self.x)+self.c2*r2*(self.pg-self.x) self.x = self.v + self.x plt.clf() plt.scatter(self.x[:, 0], self.x[:, 1], s=30, color='k') plt.xlim(self.x_bound[0], self.x_bound[1]) plt.ylim(self.x_bound[0], self.x_bound[1]) plt.pause(0.01) fitness = self.calculate_fitness(self.x) # 需要更新的个体 update_id = np.greater(self.individual_best_fitness, fitness) self.p[update_id] = self.x[update_id] self.individual_best_fitness[update_id] = fitness[update_id] # 新一代出现了更小的fitness，所以更新全局最优fitness和位置 if np.min(fitness) < self.global_best_fitness: self.pg = self.x[np.argmin(fitness)] self.global_best_fitness = np.min(fitness) print('best fitness: %.5f, mean fitness: %.5f' % (self.global_best_fitness, np.mean(fitness))) pso = PSO(100, 100) pso.evolve() plt.show()
输出
参考文献
《matlab在数学建模中的应用》
