基于粒子群算法(PCO)的函数寻优代码详解

最新推荐文章于 2022-11-28 14:11:34 发布

师范大学生

最新推荐文章于 2022-11-28 14:11:34 发布

阅读量1.2k

点赞数

分类专栏：算法程序文章标签：算法 python

本文链接：https://blog.csdn.net/m0_37872216/article/details/112531503

版权

算法程序专栏收录该内容

17 篇文章 5 订阅

订阅专栏

前言

粒子群算法与遗传算法，蚁群算法都属于启发式算法，在原理上有相当的相似性。并且粒子群算法的原理依据更为简单。
粒子群算法的原理依据和鸟群捕食的原理比较相似：假设有一群鸟一起在天上飞行捕食，每只鸟都在一边飞一边干饭。这时候鸟群中有些位置的鸟，由于所在的地方比较好，那一片天空的飞虫之类的食物多，比其他鸟吃得更饱一些。这样其他的鸟看见就坐不住了，它们会调整自己的速度和位置，让自己去靠近那些吃得饱的鸟。
最后，大部分的鸟会集中到食物最丰富的区域（理想情况下）。
由粒子群算法的理论依据可以得知，该算法适合用来解决函数寻优问题。

问题

该次解决的问题和基于遗传算法（deap）的非线性函数寻优与编码方式浅谈中的问题一样，都是求最小值。
源代码如下，根据问题做了一定改动，并优化了部分代码。
由于代码的较大部分与蚁群算法相同，这里只对粒子群算法的特色做些讲解。

import random
import numpy as np
import matplotlib.pyplot as plt
 
 
class PSO:
    def __init__(self, parameters):
        """
        particle swarm optimization
        parameter: a list type, like [NGEN, pop_size, var_num_min, var_num_max]
        """
        # 初始化
        self.NGEN = parameters[0]    # 迭代的代数
        self.pop_size = parameters[1]    # 种群大小
        self.var_num = len(parameters[2])     # 变量个数
        self.bound = []                 # 变量的约束范围
        self.bound.append(parameters[2]) # bound是个二维数组
        self.bound.append(parameters[3])
 
        self.pop_x = np.zeros((self.pop_size, self.var_num))    # 所有粒子的位置列表
        self.pop_v = np.zeros((self.pop_size, self.var_num))    # 所有粒子的速度列表
        self.p_best = np.zeros((self.pop_size, self.var_num))   # 每个粒子最优的位置列表
        self.g_best = np.zeros((1, self.var_num))   # 全局最优的位置
        # 初始化第0代初始全局最优解
        temp = 0
        for i in range(self.pop_size): # 对粒子的位置列表遍历
            for j in range(self.var_num): # 对每个粒子的变量遍历
                self.pop_x[i][j] = random.uniform(self.bound[0][j], self.bound[1][j])# 粒子的位置是边界内随机数
                self.pop_v[i][j] = random.uniform(0, 1)# 粒子的速度是一个较小随机数
            self.p_best[i] = self.pop_x[i]      # 储存每个个体的最优值，个体是粒子当前的位置
            fit = self.fitness(self.p_best[i]) # 对粒子进行适应度计算
            if i == 0:                         # 为temp赋予一个初始值（第一个粒子的适应度）
                temp = fit
            if fit > temp:                      # 用self.g_best记录当前种群里最优的粒子
                self.g_best = self.p_best[i]  
                temp = fit

    def fitness(self, ind_var):
        """
        个体适应值计算
        """
        x1 = ind_var[0]
        x2 = ind_var[1]
        y = (1.5-x1+x1*x2)*(1.5-x1+x1*x2)+ \
            (2.25-x1+x1*x2*x2)*(2.25-x1+x1*x2*x2)+ \
            (2.625-x1+x1*x2*x2*x2)*(2.625-x1+x1*x2*x2*x2)
        return -y #因为要求最小值，所以最简单的解决方式是令评价函数取负求最大值
 
    def update_operator(self, pop_size):
        """
        更新算子：更新下一时刻的位置和速度
        """
        c1 = 2     # 学习因子，一般为2
        c2 = 2
        w = 0.4    # 自身权重因子
        for i in range(pop_size):
            # 更新速度，得到位置所需要的偏移量，令粒子向种群内最优值和本身的最优值逼近
            self.pop_v[i] = w * self.pop_v[i] + c1 * random.uniform(0, 1) * (
                        self.p_best[i] - self.pop_x[i]) + c2 * random.uniform(0, 1) * (self.g_best - self.pop_x[i])
            # 使用更新的速度来更新位置
            self.pop_x[i] = self.pop_x[i] + self.pop_v[i]
            # 越界保护
            for j in range(self.var_num):
                if self.pop_x[i][j] < self.bound[0][j]:#最小值
                    self.pop_x[i][j] = self.bound[0][j]
                if self.pop_x[i][j] > self.bound[1][j]:#最大值
                    self.pop_x[i][j] = self.bound[1][j]
            # 更新p_best和g_best
            # 如果更新后的某个粒子优度比之前的高，则将该更新后的粒子记录
            # 如果更新后的某个粒子优度比之前记录的所有的粒子优度要高，则记录为最优粒子
            if self.fitness(self.pop_x[i]) > self.fitness(self.p_best[i]):
                self.p_best[i] = self.pop_x[i]
            if self.fitness(self.pop_x[i]) > self.fitness(self.g_best):
                self.g_best = self.pop_x[i]

    def main(self):#主运行函数
        popobj = []#列表记录每代的最优粒子
        self.ng_best = np.zeros((1, self.var_num))[0]
        for gen in range(self.NGEN):#循环迭代
            self.update_operator(self.pop_size)#运行粒子群算法算子
            popobj.append(self.fitness(self.g_best))#将该代最好粒子加入到popobj
            print('############ Generation {} ############'.format(str(gen + 1)))
            if self.fitness(self.g_best) > self.fitness(self.ng_best):
                self.ng_best = self.g_best.copy()
            print('最好的位置：{}'.format(self.ng_best))
            print('最小的函数值：{}'.format(self.fitness(self.ng_best)))
        print("---- End of (successful) Searching ----")
 
        plt.figure()# 使用plt模块绘制可视化图形
        plt.title("Figure1")# 题目
        plt.xlabel("iterators", size=14)# 横坐标名称
        plt.ylabel("fitness", size=14)# 纵坐标名称
        t = [t for t in range(self.NGEN)]# 横坐标的个数
        plt.plot(t, popobj, color='b', linewidth=2)# 横坐标的量，纵坐标的量（来自于列表）
        plt.show()#显示
 
 
if __name__ == '__main__':
    NGEN = 100#迭代代数
    popsize = 100#种群内粒子个数
    low = [-4.5,-4.5]#变量下界，2个变量
    up = [4.5,4.5]#变量上界
    parameters = [NGEN, popsize, low, up]#参数列表
    pso = PSO(parameters)#生成对象
    pso.main()#运行主程序

基本概念

和遗传算法中的交叉变异迭代，蚁群算法中的信息素迭代不同，粒子群算法中迭代的是位置和速度。
粒子群中，粒子等效于遗传算法中的染色体，等效于蚁群算法中的蚂蚁。显然，粒子群等效于遗传算法中的种群，等效于蚁群算法中的蚁群。

初始化函数

初始化函数中前面几个属性的概念与基于蚁群算法(ACO)的函数寻优代码详解中的一致，这里不多讲解。需要注意的是粒子群算法中需要两个“最优位置”。
self.p_best是一个和种群列表一样大小的列表，该列表存储的是每个粒子目前为止的最优值。显然粒子的位置改变具有一定随机性，并不能够保证每次迭代都比上一代的位置更好。
self.g_best是一个和一个粒子一样大小的列表，该列表存储的是当前粒子群内位置最好的粒子。

		self.p_best = np.zeros((self.pop_size, self.var_num))   # 每个粒子最优的位置列表
        self.g_best = np.zeros((1, self.var_num))   # 全局最优的位置

接下来需要构建粒子与粒子群，该步骤与蚁群算法的相同。本人在代码中做了一点改进。原先的代码中令temp =-1，作为一个临时变量的初始值，然后令后续的粒子优度与其对比，筛选出最好的粒子。
但是这样会带来一个问题，如果粒子优度都是正数，自然运行不会出问题，甚至不需要设定为-1，设定为0就可以了；但是如果粒子优度为负数，并且初始种群的优度都达不到-1，就尴尬了（这样就得不到最优粒子了，因为大家都比-1菜，都选不上）。
因此我在代码中先初始化temp为0，但是现在不能用。等到第一个粒子有了优度之后将第一个粒子的优度直接赋值给temp，这样优度的比较就只在粒子群内部进行了，和外部设置无关。
通过这里我也认识到，在程序中设置一个参考定量是件需要无比谨慎的事情。

		temp = 0
        for i in range(self.pop_size): # 对粒子的位置列表遍历
            for j in range(self.var_num): # 对每个粒子的变量遍历
                self.pop_x[i][j] = random.uniform(self.bound[0][j], self.bound[1][j])# 粒子的位置是边界内随机数
                self.pop_v[i][j] = random.uniform(0, 1)# 粒子的速度是一个较小随机数
            self.p_best[i] = self.pop_x[i]      # 储存每个个体的最优值，个体是粒子当前的位置
            fit = self.fitness(self.p_best[i]) # 对粒子进行适应度计算
            if i == 0:                         # 为temp赋予一个初始值（第一个粒子的适应度）
                temp = fit
            if fit > temp:                      # 用self.g_best记录当前种群里最优的粒子
                self.g_best = self.p_best[i]  
                temp = fit

适应度函数

这次要求函数的最小值，而程序中比较的方式都是比较谁的优度更大，也就是比谁的函数值更大。
这里显然有两种修改思路：1.将程序中所有比优度大大的语句都改成比优度小的语句。该做法显然是可行的，但是在理论上有些违背“寻求优度最大”的基础概念，并且这样更改修改的幅度大一些。2.第二种方法参考自之前的遗传算法，如果对函数整体取负，那么求最小值就成了求最大值，与原先的代码思路相对契合。
于是在适应度计算函数中，将返回值取负。

 		def fitness(self, ind_var):
        """
        个体适应值计算
        """
        x1 = ind_var[0]
        x2 = ind_var[1]
        y = (1.5-x1+x1*x2)*(1.5-x1+x1*x2)+ \
            (2.25-x1+x1*x2*x2)*(2.25-x1+x1*x2*x2)+ \
            (2.625-x1+x1*x2*x2*x2)*(2.625-x1+x1*x2*x2*x2)
        return -y #因为要求最小值，所以最简单的解决方式是令评价函数取负求最大值

核心算子

粒子群核心算子主要分两步。
第一步，更新该粒子的速度。从以下代码中可以看到，速度的求解方式为：在部分保留自己上一个速度的同时，令自己的值向着该粒子历史上的最优值和粒子群中的最优值逼近。
第二步，直接将速度作为偏移量叠加到位置上得到新的位置。
可以看出，如果当前粒子位置与当前种群内的最优粒子位置和自己的历史最好位置差距较大时，该粒子就会增加移动的幅度，令自己快速逼近到一个较优值上；反之，当前粒子位置与当前种群内的最优粒子位置和自己的历史最好位置差距较小时，则该粒子就会减少移动的幅度。

# 更新速度，得到位置所需要的偏移量，令粒子向种群内最优值和本身的最优值逼近
            self.pop_v[i] = w * self.pop_v[i] + c1 * random.uniform(0, 1) * (
                        self.p_best[i] - self.pop_x[i]) + c2 * random.uniform(0, 1) * (self.g_best - self.pop_x[i])

在更新位置后再用位置计算优度，更新优度存储表。
一个粒子的位置，就是该粒子的解。

		def update_operator(self, pop_size):
        """
        更新算子：更新下一时刻的位置和速度
        """
        c1 = 2     # 学习因子，一般为2
        c2 = 2
        w = 0.4    # 自身权重因子
        for i in range(pop_size):
            # 更新速度，得到位置所需要的偏移量，令粒子向种群内最优值和本身的最优值逼近
            self.pop_v[i] = w * self.pop_v[i] + c1 * random.uniform(0, 1) * (
                        self.p_best[i] - self.pop_x[i]) + c2 * random.uniform(0, 1) * (self.g_best - self.pop_x[i])
            # 使用更新的速度来更新位置
            self.pop_x[i] = self.pop_x[i] + self.pop_v[i]
            # 越界保护
            for j in range(self.var_num):
                if self.pop_x[i][j] < self.bound[0][j]:#最小值
                    self.pop_x[i][j] = self.bound[0][j]
                if self.pop_x[i][j] > self.bound[1][j]:#最大值
                    self.pop_x[i][j] = self.bound[1][j]
            # 更新p_best和g_best
            # 如果更新后的某个粒子优度比之前的高，则将该更新后的粒子记录
            # 如果更新后的某个粒子优度比之前记录的所有的粒子优度要高，则记录为最优粒子
            if self.fitness(self.pop_x[i]) > self.fitness(self.p_best[i]):
                self.p_best[i] = self.pop_x[i]
            if self.fitness(self.pop_x[i]) > self.fitness(self.g_best):
                self.g_best = self.pop_x[i]

测试结果

在这里插入图片描述
可以看到由于问题比较简单，所以在不到20代的时候就收敛到最优解了。
最后想一下这个学习因子的事，如果学习因子较大，粒子的动作幅度就相对较大，收敛的速度快，但是容易陷入局部最优；学习因子较小，粒子的动作幅度就相对较小，收敛的速度慢，但是能做到尽可能增加搜索空间，相对容易找到全局最优解。
以上是个人猜想，后来一查资料发现，这两个学习因子都取[0,4]区间里的整数就行，啊这。
实际上惯性因子的影响力更大，惯性因子大时全局搜索能力强，惯性因子小时局部搜索能力强。
本质上也还是，自身变动小，收敛慢，搜索空间大，全局搜索能力强；自身变动大，收敛快，搜索空间小，局部搜索能力强。