数学建模之粒子群算法

（以下为培训学习内容）

群体迭代，粒子在解空间追随最优的粒子进行搜索

 

粒子群算法的思想源于对鸟群捕食行为的研究。

鸟宝宝捕食的故事，正是这个粒子群算法存在的原因。因此，如果想更好的了解粒子群算法，我们就要来分析鸟宝宝捕食的故事。首先，我们来分析分析鸟宝宝们的运动状态，即鸟宝宝自身是怎么决定自己的飞翔速度和位置的。

　　(1) 呐首先，我们知道物体是具有惯性的，鸟宝宝在一开始飞翔的时候，无论它下一次想怎么飞，往哪个方向飞，它都有一个惯性，它必须根据当前的速度和方向来进行下一步的调整。对吧，这个可以理解吧，因此，“惯性”——当前的速度vcurrentvcurrent是一个因素。

　　(2) 其次，由于鸟宝宝长期捕食，因此鸟宝宝有经验，它虽然不知道具体哪里是有虫子的存在，但是它能大概知道虫子分布在哪里。比如当鸟宝宝飞到贫瘠的地方，它肯定知道这里是不会有虫子的，因此，在鸟宝宝的心中，它每次飞，都会根据自己的经验来找，比如以往虫子分布的地区，它肯定优先对这部分的地方进行搜索。因此，自己的“认知”——经验，也是一个因素。

　　(3) 最后，每个鸟宝宝发现自己离虫子更接近的时候，便会发出信号给同伴，在遇到这样的信号，其余还在找的鸟宝宝们就会想着，同伴的位置更接近虫子，我要往它那边过去看看。我们称之为“社会共享”，这也是鸟宝宝在移动时考虑的一个因素。

　　综上所述，鸟宝宝每次在决定下一步移动的速度和方向时，脑子里是由这三个因素影响的。我们想，如果能够用一条公式来描述着三个因素的影响的话，那不就能写出每个鸟宝宝的移动方向和速度么。但是，每一个鸟宝宝，对这三个因素的考虑都是不一致的，比如对于捕食经验不高的鸟宝宝，移动的时候会更看重同伴分享的信息，而对于捕食经验高的鸟宝宝，则更看重自己的经验。因此，我们的公式，在“认知”和“社会共享”这部分，是具有随机性（本故事引用https://www.cnblogs.com/Qling/p/9343625.html）

停止的两种方法：1、鸟群飞累了，寻找过程可以了，就会根据判断停到最优解上。2、鸟群寻找的最优解很久没有更新了（发现比之前最优解更优时会更新）。

 

 

D维空间一般最多三维或二维空间进行。

 

                                                               （会往区域最优解和全局最优解中间飞）

 

流程图如下：

构成要素：

1、群体大小m

m是一个整型参数

m很小：陷入局优的可能性很大

比如有一个大农场和一个小农场，就两个鸟可能只发现了小农场。

m很大：POS的优化能力很好，但是当增大到一个数值后，再增加也不会增加效率和优化能力，只会降低效率。

2、权重因子：

 

  

c1=0太依靠群体，不相信自己；c2=0太相信自我，不相信群体。

如果c1，c2都不为0，则成为完全型粒子群算法。

完全型粒子群算法：

3、最大速度：

Vm比较大时----->探索能力增强，但是粒子容易飞过最优解

Vm比较小时----->开发能力增强，但是粒子容易陷入局部最优解

Vm一般设为每维变量变化范围的10%-20%

所以根据领域拓扑结构分为全局粒子群算法和局部粒子群算法

综上所述，相比于其他算法，所需要的参数较少，但是惯性因子和领域定义特别重要

解答过程：

有初始位置和初始速度可以计算每个粒子的适应值：

结算结果如下：

                                 

选取最小值作为历史最优解

个体历史最优解仍然认为自己是最优解

（因为最大速度为60，所以要将超过这个范围的速度限制到60上）

重复上述步骤，将迭代进行下去

粒子群算法的步骤如下：

惯性权重：

                                                                          （最大惯性权重和最小惯性权重都是初始定义的）

收缩因子法用的比较少。

import numpy as np
import matplotlib.pyplot as plt
import random

 
 
# 定义“粒子”类
class parti(object):
    def __init__(self, v, x):
        self.v = v                    # 粒子当前速度
        self.x = x                    # 粒子当前位置
        self.pbest = x                # 粒子历史最优位置
        
class PSO(object):
    def __init__(self, interval, tab, partisNum=10, iterMax=1000, w=1, c1=2, c2=2):
        self.interval = interval                                            # 给定状态空间 - 即待求解空间
        self.tab = tab                                              # 求解最大值还是最小值的标签: 'min' - 最小值；'max' - 最大值
        self.iterMax = iterMax                                              # 最大迭代求解次数
        self.w = w                                                          # 惯性因子
        self.c1, self.c2 = c1, c2                                           # 学习因子
        self.v_max = (interval[1] - interval[0]) * 0.1                      # 设置最大迁移速度
        self.partis_list, self.gbest = self.initPartis(partisNum)                 # 完成粒子群的初始化，并提取群体历史最优位置
        self.x_seeds = np.array(list(parti_.x for parti_ in self.partis_list))    # 提取粒子群的种子状态 
        self.solve()                                                              # 完成主体的求解过程
        self.display()                                                            # 数据可视化展示
        
    def initPartis(self, partisNum):
        partis_list = list()
        for i in range(partisNum):
            v_seed = random.uniform(-self.v_max, self.v_max)
            x_seed = random.uniform(*self.interval)
            partis_list.append(parti(v_seed, x_seed))                       #确定目标找最大值还是最小值
        temp = 'find_' + self.tab
        if hasattr(self, temp):                                            
            gbest = getattr(self, temp)(partis_list)
        else:
            exit('>>>tab标签传参有误："min"|"max"<<<')
        return partis_list, gbest
        
    def solve(self):
        for i in range(self.iterMax):
            for parti_c in self.partis_list:
                f1 = self.func(parti_c.x)
                # 更新粒子速度，并限制在最大迁移速度之内
                parti_c.v = self.w * parti_c.v + self.c1 * random.random() * (parti_c.pbest - parti_c.x) + self.c2 * random.random() * (self.gbest - parti_c.x)
                if parti_c.v > self.v_max: 
                    parti_c.v = self.v_max
                elif parti_c.v < -self.v_max: 
                    parti_c.v = -self.v_max
                # 更新粒子位置，并限制在待解空间之内
                if self.interval[0] <= parti_c.x + parti_c.v <=self.interval[1]:
                    parti_c.x = parti_c.x + parti_c.v 
                else:
                    parti_c.x = parti_c.x - parti_c.v
                f2 = self.func(parti_c.x)
                getattr(self, 'deal_'+self.tab)(f1, f2, parti_c)             # 更新粒子历史最优位置与群体历史最优位置      
        
    def func(self, x):                                                       # 状态产生函数 - 即待求解函数
        value = np.sin(x**2) * (x**2 - 5*x)
        return value
        
    def find_min(self, partis_list):                                         # 按状态函数最小值找到粒子群初始化的历史最优位置
        parti = min(partis_list, key=lambda parti: self.func(parti.pbest))
        return parti.pbest
        
    def find_max(self, partis_list):
        parti = max(partis_list, key=lambda parti: self.func(parti.pbest))   # 按状态函数最大值找到粒子群初始化的历史最优位置
        return parti.pbest
        
    def deal_min(self, f1, f2, parti_):
        if f2 < f1:                          # 更新粒子历史最优位置
            parti_.pbest = parti_.x
        if f2 < self.func(self.gbest):
            self.gbest = parti_.x            # 更新群体历史最优位置
            
    def deal_max(self, f1, f2, parti_):
        if f2 > f1:                          # 更新粒子历史最优位置
            parti_.pbest = parti_.x
        if f2 > self.func(self.gbest):
            self.gbest = parti_.x            # 更新群体历史最优位置
            
    def display(self):
        print('solution: {}'.format(self.gbest))
        plt.figure(figsize=(8, 4))
        x = np.linspace(self.interval[0], self.interval[1], 300)
        y = self.func(x)
        plt.plot(x, y, 'g-', label='function')
        plt.plot(self.x_seeds, self.func(self.x_seeds), 'b.', label='seeds')
        plt.plot(self.gbest, self.func(self.gbest), 'r*', label='solution')
        plt.xlabel('x')
        plt.ylabel('f(x)')
        plt.title('solution = {}'.format(self.gbest))
        plt.legend()
        plt.savefig('PSO.png', dpi=500)
        plt.show()
        plt.close()
 
        
if __name__ == '__main__':
    PSO([-9, 5],'max')#

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.max(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()

（参考别人的代码，暂存整理用的）