pso-svm 算法实现（1）：python DEAP

粒子群优化算法（Particle Swarm Optimization，PSO）属于进化算法的一种，是通过模拟鸟群捕食行为设计的。从随机解出发，通过迭代寻找最优解，通过适应度来评价解的品质。设想这样一个场景：一群鸟在随机搜索食物。在这个区域里只有一块食物。所有的鸟都不知道食物在那里。但是他们知道当前的位置离食物还有多远。那么找到食物的最优策略是什么呢。最简单有效的就是搜寻目前离食物最近的鸟的周围区域。
所有的粒子具有以下两个属性：速度、位置。
PSO 初始化为一群随机粒子(随机解)。然后通过迭代找到最优解。在每一次迭代中，粒子通过跟踪两个”极值”来更新自己。第一个就是粒子本身所找到的最优解pbest。另一个极值是整个种群目前找到的最优解，即全局极值gbest。粒子通过下面的公式来更新自己的速度和位置：

速度：
vi=w∗vi+c1∗rand1∗(pbesti−xi)+c2∗rand2∗(gbesti−xi)vi=w∗vi+c1∗rand1∗(pbesti−xi)+c2∗rand2∗(gbesti−xi)
位置：
xi=xi+vi+1xi=xi+vi+1
以上两式中：
ww为惯性因子，一般取1
c1、c2c1、c2为学习因子，一般取2
rand1、rand2rand1、rand2为两个（0,1）之间的随机数
vi和xivi和xi分别表示粒子第i维的速度和位置
pbesti、gbestipbesti、gbesti分别表示某个粒子最好位置第i维的值、整个种群最好位置第i维的值

注意：以上两式是针对粒子的某一维进行跟新的。对粒子的每一维，都要用上述的式子进行更新

说了这么多，下面把主要部分代码上上来（不上全是因为伸手党太多，哈哈）：

import numpy as np
import random

class PSO:
def __init__(self, dim, size, iter_num, x_max, max_vel, best_fitness_value=float('Inf'), C1 = 2, C2 = 2, W = 1):
self.C1 = C1
self.C2 = C2
self.W = W
self.dim = dim # 粒子的维度
self.size = size # 粒子个数
self.iter_num = iter_num # 迭代次数
self.x_max = x_max
self.max_vel = max_vel # 粒子最大速度
self.best_fitness_value = best_fitness_value
self.best_position = [0.0 for i in range(dim)] # 种群最优位置
self.fitness_val_list = [] # 每次迭代最优适应值

# 对种群进行初始化
self.Particle_list = [Particle(self.x_max, self.max_vel, self.dim) for i in range(self.size)]

# 更新速度
def update_vel(self, part):
for i in range(self.dim):
vel_value = self.W * part.get_vel()[i] + self.C1 * random.random() * (part.get_best_pos()[i] - part.get_pos()[i]) \
+ self.C2 * random.random() * (self.get_bestPosition()[i] - part.get_pos()[i])
if vel_value > self.max_vel:
vel_value = self.max_vel
elif vel_value < -self.max_vel:
vel_value = -self.max_vel
part.set_vel(i, vel_value)

# 更新位置
def update_pos(self, part):
for i in range(self.dim):
pos_value = part.get_pos()[i] + part.get_vel()[i]
part.set_pos(i, pos_value)
value = fit_fun(part.get_pos())
if value < part.get_fitness_value():
part.set_fitness_value(value)
for i in range(self.dim):
part.set_best_pos(i, part.get_pos()[i])
if value < self.get_bestFitnessValue():
self.set_bestFitnessValue(value)
for i in range(self.dim):
self.set_bestPosition(i, part.get_pos()[i])

def update(self):
for i in range(self.iter_num):
for part in self.Particle_list:
self.update_vel(part) # 更新速度
self.update_pos(part) # 更新位置
self.fitness_val_list.append(self.get_bestFitnessValue()) # 每次迭代完把当前的最优适应度存到列表
return self.fitness_val_list, self.get_bestPosition()

在上述代码中，跟粒子的速度设置了最大值，确保粒子速度的变化在某一范围内。接下来，我们可以用以上代码来进行测试了，这里测试用的是Hölder table function测试函数，在我的另一篇博客用python绘制评估优化算法性能的测试函数里有提到，其图像如下：

它的最优解有如下4个：

接下了我们通过以下代码进行测试：

from OptAlgorithm.PSO import PSO
import matplotlib.pyplot as plt
import numpy as np

dim = 2
size = 20
iter_num = 1000
x_max = 10
max_vel = 0.5

pso = PSO(dim, size, iter_num, x_max, max_vel)
fit_var_list, best_pos = pso.update()
print("最优位置:" + str(best_pos))
print("最优解:" + str(fit_var_list[-1]))
plt.plot(np.linspace(0, iter_num, iter_num), fit_var_list, c="R", alpha=0.5)
plt.show()

以上代码中，设置粒子维度为2维，种群大小即粒子个数为20，迭代次数为1000，粒子初始化的位置在（-10,10）之间，粒子的最大速度为0.5

以下是一次算法的收敛图像结果，：

打印的最优解和适应度如下：

可以看出其最优位置是函数四个最优解里面的一个，运行多次会出现不同的结果。甚至可能会得到如下图所示的函数收敛图像：

注意看此时函数的适应度在-9.5左右，继续迭代也没有发生变化，此时即陷入了局部最优。

接下来解决pso-svm算法的实现。

在pso-svm 算法中，参数寻优步骤大概如下：

（1）PSO算法寻找适用于SVM模型的核函数类型
Step 1：初始化粒子群规模m，设定算法的权重因子，终止条件和初始粒子编码；
Step 2：将每个粒子的个体极值设置为当前位置，利用**适应度函数**计算每个粒子的适应度值，取适应度好的粒子做，
对应的个体极值作为最初的全局极值；
Step 3：按照粒子的位置和速度更新公式进行迭代计算，更新粒子的位置和速度；
Step 4：按照粒子的**适应度函数**计算每次迭代后每个粒子的适应度值；
Step 5：将每个粒子的适应度值与其个体极值的适应度值作比较，如果更优的话，则更新个体极值，否则保留原值；
Step 6：将更新后的每个粒子的个体极值与全局极值比较，如果更优的话，则更新全局极值，否则保留原值；
Step 7：判断是否满足终止条件，若达到最大迭代次数或者所得解收敛或者所得解已经达到了预期的效果，就终止迭代，否则返回Step 3
Step 8：得到使得模型最佳的参数组合，用于构建子最优模型。

根据上述步骤，“适应度函数”就是判断svm模型的准确率和泛化能力的函数，在SVR模型中，一般采用MSE均方误差函数作为适应度函数，模型越好，那么适应度越高。更详细的说，就是对每次迭代过程中的每个粒子，将其表示的参数值代入到模型中，训练得到SVM模型，然后利用测试数据来评估模型的MSE，这个MSE就作为最终得到的适应度值，用这个值来比较粒子的好坏，即比较这一组参数的好坏。
在DEAP框架的PSO代码中，有如下代码：

creator.create("FitnessMax", **base.Fitness**, weights=(1.0,))
creator.create("Particle", list, fitness=creator.FitnessMax, speed=list, smin=None, smax=None, best=None)

第一行创建了适应度函数，第二行创建了粒子，为了修改适应度函数为SVM相关的MSE，考察相关的base.Fitness 类。
链接如下：http://deap.gel.ulaval.ca/doc/dev/api/base.html#fitness
翻译原文：

Fitness

class deap.base.Fitness([values])
这个fitness(适应度值)可以测量一个解的好坏。如果values是一个元祖类型，那么fitness使用这些值进行初始化，否则为空或者无效。
参数：values fitness的初始值元祖，可选项。
fitness可以使用这些符号比较：>, <, >=, <=, ==, !=。这些算子的比较是按照字典序进行的。最大化和最小化通过权重和适应度值之间的乘法来处理。对不同大小的适应度值进行比较，如果每个元素的适应度都相等，那么适应度越高越好。
不同类型的适应度在 Creating Types（http://deap.gel.ulaval.ca/doc/dev/tutorials/basic/part1.html#creating-types）手册中创建。
注意：如果比较最小化的适应度值a>b,如果a小于b那么返回True。

接下来查看 Creating Types手册，学习其他类型适应度函数。链接如上。
翻译如下：

Creating Types

这个手册展示如何用creator创建类型和用toolbox初始化。（这里只介绍和PSO算法相关的）

Fitness 适应度

提供的Fitness 类是一个抽象的类，需要weights属性来保证功能。使用负的weights来构建一个最小化的fitness，正的weights构建最大化的fitness。例如，在接下来的代码中，在creator中，使用单目标最小化fitness，命名为FitnessMin。

creator.create("FitnessMin", base.Fitness, weights=(-1.0,))

create函数需要至少两个参数，一个是新创建的类，另一个是基类。随后的参数是这个类的属性。根据之前介绍的Fitness文档，属性weights必须是元祖，所以多目标和单目标的fitness可以使用相同的方式使用。一个FitnessMulti 可以使用如下方式创建：

creator.create("FitnessMulti", base.Fitness, weights=(-1.0, 1.0))

这行代码产生一个适应度值，最小化第一个目标，最大化第二个目标。weights属性也可以用来该表每个目标的重要性。这意味着weights可以使任何的实数，用来决定最大化还是最小化。有一个例子，在NSGA-II 选择算法中，在拥挤距离排序中使用了这种weights。

Individual 个体

简单思考不同进化算法（GA，GP，PSO，ES，。。。）后，我们注意到有很多种类的个体是可能的，假设开发人员无法获得所有的类型。这里有一个指南关于如何使用creator创建这些个体，和如何使用Toolbox初始化。
warning： 在继承numpy.ndarray之前，你应该阅读Inheriting from Numpy手册，并且看一下One Max Problem: Using Numpy 示例。

Particle粒子

一个粒子是另一种类型的个体，它通常有速度，还要记录最好的位置。这种个体的创建方式同样继承于list。这时，属性速度speed，最优位置best，速度极值都被添加到对象中。而且，一个初始化函数initParticle()被注册，接受参数particle类，大小，阈值，速度极值从而产生个体（粒子）。

import randomfrom deap import base
from deap import creator
from deap import toolscreator.create("FitnessMax", base.Fitness, weights=(1.0, 1.0))
creator.create("Particle", list, fitness=creator.FitnessMax, speed=None,smin=None, smax=None, best=None)def initParticle(pcls, size, pmin, pmax, smin, smax):part = pcls(random.uniform(pmin, pmax) for _ in xrange(size))part.speed = [random.uniform(smin, smax) for _ in xrange(size)]part.smin = sminpart.smax = smaxreturn parttoolbox = base.Toolbox()
toolbox.register("particle", initParticle, creator.Particle, size=2,pmin=-6, pmax=6, smin=-3, smax=3)

调用toolbox.individual（）将很容易地返回一个完整的粒子，速度矢量和最大化两个目标适应度属性。

Population种群：Swarm

swarm在PSO算法（particle swarm optimization）中使用。它不同的地方在于它包含了一个通信网络。最简单的网络是全连接，在全连接网络中每一个粒子都知道任何粒子都曾经访问过的最好的位置best。它的实现通常是复制全局最优位置到gbest属性，复制全局最优适应度值到gbestfit属性。

creator.create("Swarm", list, gbest=None, gbestfit=creator.FitnessMax)
toolbox.register("swarm", tools.initRepeat, creator.Swarm, toolbox.particle)

调用 toolbox.population()将会返回一个完整的swarm。在每一次评估之后，gbest 和 gbestfit按照算法，应该是要反映最好的位置和适应度值。

看完了这些文档，还存在两点问题：
（1）找不到调用其它适应度函数的接口；
（2）找不到粒子编译码的设置；
综上，放弃这种做法，寻找其他方法。