PSO粒子群算法的python简单实现备忘录

主要会出现崩溃问题，所以我就把参考的博文写在前面了：

http://blog.csdn.net/niuyongjie/article/details/1569671

http://blog.csdn.net/kunshanyuz/article/details/63683145

http://blog.csdn.net/zuochao_2013/article/details/53431767?ref=myread

粒子群算法源于复杂适应系统（Complex Adaptive System,CAS）。CAS理论于1994年正式提出，CAS中的成员称为主体。比如研究鸟群系统，每个鸟在这个系统中就称为主体。主体有适应性，它能够与环境及其他的主体进行交流，并且根据交流的过程“学习”或“积累经验”改变自身结构与行为。整个系统的演变或进化包括：新层次的产生（小鸟的出生）；分化和多样性的出现（鸟群中的鸟分成许多小的群）；新的主题的出现（鸟寻找食物过程中，不断发现新的食物）。

 我觉得简单来讲，就是有一个函数，你想求它的最大值或最小值，直接求很难求解，你就只好用迂回地方法，设置好多点，然后计算对比他们的值，然后不断迭代，直到找到使得函数得到最大值或最小值的点。

比如我们求取函数y=1-cos(3*x)*exp(-x)的在[0,4]最大值。并在[0,4]之间放置了两个随机的点，这些点的坐标假设为x1=1.5； x2=2.5；这里的点是一个标量，但是我们经常遇到的问题可能是更一般的情况－－x为一个矢量的情况，比如二维的情况 z=2*x1+3*x2²的情况。这个时候我们的每个粒子为二维，记粒子P1＝(x11,x12),P2=(x21,x22),P3=(x31,x32)，......Pn=(xn1,xn2)。这里n为粒子群群体的规模，也就是这个群中粒子的个数，每个粒子的维数为2。更一般的是粒子的维数为q，这样在这个种群中有n个粒子，每个粒子为q 维。

由n个粒子组成的群体对Q维（就是每个粒子的维数）空间进行搜索。每个粒子表示为：x_i＝（x_i1,x_i2,x_i3,...,x_iQ），每个粒子对应的速度可以表示为v_i=(v_i1,v_i2,v_i3,....,v_iQ)，每个粒子在搜索时要考虑两个因素：

1. 自己搜索到的历史最优值 p_{i ，}p_i=(p_i1,p_i2,....,p_iQ)，i=1,2,3,....,n。

2. 全部粒子搜索到的最优值p_g，p_g=(p_g1,p_g2,....,p_gQ)，注意这里的p_g只有一个。

#coding= utf-8import numpy as npimport randomimport matplotlib.pyplot as plt"""PSO粒子类（本程序只设置了最大迭代轮数，而未设置其他限制条件）#p:self；粒子数量；搜索维度（也就是一个粒子里有几个变量）；最大迭代次数"""class PSO():    def __init__(self, pN, dim, max_iter):        self.w = 0.8    #保持原来速度的系数，所以叫做惯性权重        self.c1 = 2     #粒子跟踪自己历史最优值的权重系数，它表示粒子自身的认识，所以叫“认知”。通常设置为2        self.c2 = 2     #粒子跟踪群体最优值的权重系数，它表示粒子对整个群体知识的认识，所以叫做“社会知识”，经常叫做“社会”。通常设置为2        self.r1 = 0.6   #[0,1]区间内均匀分布的伪随机数        self.r2 = 0.3   #[0,1]区间内均匀分布的伪随机数        self.r3 = 1     #对位置更新的时候，在速度前面加的一个系数，这个系数我们叫做约束因子。通常设置为1        self.pN = pN    #粒子数量        self.dim = dim  #搜索维度        self.max_iter = max_iter  #最大迭代次数        self.X = np.zeros((self.pN, self.dim))      #所有粒子的位置        self.V = np.zeros((self.pN, self.dim))      #所有粒子的速度        self.pbest = np.zeros((self.pN, self.dim))  #个体经历的最佳位置        self.gbest = np.zeros((1, self.dim))        #全局最佳位置        self.p_fit = np.zeros(self.pN)  #每个个体的历史最佳适应值        self.fit = 1e10                 #全局最佳适应值    """    #适应度的计算（对于函数而言，其实就是计算y的值），这个例子用的就是简单地将其平方再加和    #p:self；粒子数量；粒子的位置信息    #r:适应度的值    """    def function(self, x):        sum = 0        length = len(x)        x = x ** 2  #将每一个维度的x值求平方        #对各个维度的平方求加和        for i in range(length):            sum += x[i]        return sum    """    #初始化粒子信息    #p:self    #r:无；在过程中更新粒子所需的信息    """    def init_Population(self):        #遍历所有粒子        for i in range(self.pN):            #遍历所有维度            for j in range(self.dim):                self.X[i][j] = random.uniform(0, 1)  #随机初始化位置信息                self.V[i][j] = random.uniform(0, 1)  #随机初始化速度信息            #初始化个体经历的最佳位置            self.pbest[i] = self.X[i]            #每个个体的历史最佳适度            tmp = self.function(self.X[i])            self.p_fit[i] = tmp            #对比各个粒子的适应度，选择最高的作为全局最优            if (tmp < self.fit):                self.fit = tmp                self.gbest = self.X[i]    """    #算法主体，计算更新位置信息、速度信息和适应度，并给出最终结果    #p:self    #r:最佳适应度（也就是所求的y值）    """    def iterator(self):        fitness = []  #制作一个数组存放每轮的最优值        #在最大迭代数内迭代        for t in range(self.max_iter):            #遍历粒子数，计算适应度            for i in range(self.pN):                temp = self.function(self.X[i])                #如果所得的适应度比之前得到的个体最优适应度要小，就更新个体最优                if (temp < self.p_fit[i]):                    self.p_fit[i] = temp                    self.pbest[i] = self.X[i]                    #如果更新的个体最优小于全局最优，那么就用个体最优更新全局最优                    if (self.p_fit[i] < self.fit):  # 更新全局最优                        self.gbest = self.X[i]                        self.fit = self.p_fit[i]            #更新完适应度的个体最优和全局最优后，根据公式更新位置学习和速度学习            for i in range(self.pN):                self.V[i] = self.w * self.V[i] + self.c1 * self.r1 * (self.pbest[i] - self.X[i]) \                            + self.c2 * self.r2 * (self.gbest - self.X[i])                self.X[i] = self.X[i] + self.r3*self.V[i]            #将每轮得到的最佳适应度存放在fitness数组中            fitness.append(self.fit)            print(self.fit)  #输出最优值        return fitness# ----------------------程序执行-----------------------my_pso = PSO(pN=30, dim=5, max_iter=50)my_pso.init_Population()fitness = my_pso.iterator() # -------------------画图--------------------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, 50)])fitness = np.array(fitness)plt.plot(t, fitness, color='b', linewidth=3)plt.show()