优化算法(三):粒子群算法

引言

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

所以CAS系统中的主体具有4个基本特点（这些特点是粒子群算法发展变化的依据）：

• 首先，主体是主动的、活动的。
• 主体与环境及其他主体是相互影响、相互作用的，这种影响是系统发展变化的主要动力。
• 环境的影响是宏观的，主体之间的影响是微观的，宏观与微观要有机结合。
• 最后，整个系统可能还要受一些随机因素的影响。

粒子群算法就是对一个CAS系统－－－鸟群社会系统的研究得出的。

一、粒子群算法

粒子群算法通过设计一种无质量的粒子来模拟鸟群中的鸟，**粒子仅具有两个属性：速度和位置**，速度代表移动的快慢，位置代表移动的方向。每个粒子在**搜索空间**中单独的搜寻最优解，并将其记为**当前个体极值**，并将个体极值与整个粒子群里的其他粒子共享，**找到最优的那个个体极值作为整个粒子群的当前全局最优解**，粒子群中的所有粒子根据自己找到的当前个体极值和整个粒子群共享的当前全局最优解来**调整自己的速度和位置**。粒子群算法的思想相对比较简单，主要分为：

• 1. 初始化粒子群；
• 2. 评价粒子，即计算适应值；
• 3. 更新个体极值；
• 4. 更新全局最优解；
• 5. 更新粒子的速度和位置并迭代(b)~(e)。

（1）初始化 首先，我们需要设置最大的速度区间，防止超出最大的区间。位置信息即为整个搜索空间，我们在速度区间和搜索空间上随机初始化速度和位置，设置群体规模。

（2）个体极值与全局最优解 个体极值为每个粒子找到的历史上最优的位置信息，并从这些个体历史最优解中找到一个全局最优解，并与历史最优解比较，选出最佳的作为当前的历史最优解。

（3）更新速度和位置的公式

$$V_{id}^{n+1}=\omega V_{id}^n+C_{1}ran{d}_1(P_{id}-X_{id}^n)+C_{2}ran{d}_2(P_{gd}-X_{id}^n))$$

$$X_{id}^{n+1}=X_{id}^n+r{V}_{id}^{n+1}$$

• $\omega$ 是保持原来速度的系数，所以叫做惯性权重；
• $C_1$ 是粒子跟踪自己历史最优值的权重系数，它表示粒子自身的认识，称为“自身认知”，通常设置为2；
• $C_2$ 是粒子跟踪群体最优值的权重系数，它表示粒子对整个群体知识的认识，所以叫做“社会知识”，通常设置为2；
• $rand1$ 、$rand2$ 是 区间 $[0,1]$ 上的随机浮点数；
• $r$ 是对位置更新的时候，在速度前面加的一个系数，这个系数我们叫做”约束因子“。通常设置为1；
• $P_{id}$ 是粒子个体历史上最优值对应的位置信息；$P_{gd}$ 是粒子群中最优值对应的位置信息；

（4）终止条件 有两种终止条件可以选择，一是设置最大迭代次数；二是相邻两代之间的偏差在一个指定的范围内。

（注：位置信息可能是多维度的，比如 $O-xy$ 位置就是一维； $O-xyz$ 位置就是二维…）

二、实例

研究函数 $y=xsin(10\pi x)+1$ 在区间 $[-1,1]$ 上的最值问题。

import numpy as np
import matplotlib.pyplot as plt


class Particle():
    def __init__(self, position, speed, best_fitness, best_position):
        self.position = position               # 个体当前位置
        self.speed = speed                     # 个体当前速度
        self.best_fitness = best_fitness       # 个体历史最优适应值
        self.best_position = best_position     # 个体历史最优适应值对应的位置

class PSO():
    def __init__(self, min_speed, max_speed, Iternumbers, C1, C2, R, W, particleNum):
        self.min_speed = min_speed       # 系统最小速度
        self.max_speed = max_speed       # 系统最大速度
        self.Iternumbers = Iternumbers   # 系统迭代次数
        self.C1 = C1                     # 自身认知
        self.C2 = C2                     # 社会认知
        self.R = R                       # 约束因子
        self.W = W                       # 惯性因子
        self.particleNum = particleNum   # 粒子个数
        self.best_position = None        # 系统历史最优适应值
        self.best_fitness = -100000      # 系统历史最优适应值对应的位置

    def initial(self):    # 初始化粒子群
        particles = []
        for _ in range(self.particleNum): # 每轮创建随机粒子
            # 初始化随机位置
            rand_position = np.random.rand() * (self.max_speed-self.min_speed) + self.min_speed
            # 初始化随即速度
            rand_speed = np.random.rand() * (self.max_speed-self.min_speed) + self.min_speed
            # 计算适应值
            fitness = self.compute_fitness(rand_position)
            # 更新系统历史最优适应值和位置
            if fitness > self.best_fitness:
                self.best_fitness = fitness
                self.best_position = rand_position
            # 创建粒子
            particle = Particle(rand_position, rand_speed, fitness, rand_position)
            particles.append(particle)
        return particles

    def compute_fitness(self, x):
        return x * np.sin(10*np.pi*x) + 1

    def update_speed(self, particle):
        rand1 = np.random.rand()
        rand2 = np.random.rand()
        V_next = self.W * particle.speed + self.C1 * rand1 * (particle.best_position - particle.position) \
                 + self.C2 * rand2 * (self.best_position - particle.position)
        if V_next > self.max_speed:
            V_next = self.max_speed
        elif V_next < self.min_speed:
            V_next = self.min_speed
        else:
            V_next = V_next
        return V_next

    def update_position(self, particle, V_next):
        X_next = particle.position + self.R * V_next
        if X_next > self.max_speed:
            X_next = self.max_speed
        elif X_next < self.min_speed:
            X_next = self.min_speed
        else:
            X_next = X_next
        return X_next

    def update(self, particles):
        temp_fitness = 0     # 记录当前更新的最优适应值及位置
        temp_position = 0
        for particle in particles: # 遍历每个粒子
            V_next = self.update_speed(particle)
            X_next = self.update_position(particle,V_next)
            particle.position = X_next   # 更新粒子的位置及速度
            particle.speed = V_next
            fitness = self.compute_fitness(X_next)
            if fitness > particle.best_fitness:  # 如果适应值大于个体历史适应值，进行更新
                particle.best_fitness = fitness
                particle.best_position = X_next
                if fitness > self.best_fitness:  # 如果大于系统历史适应值，保留当前值
                    temp_fitness = fitness
                    temp_position = X_next
        if temp_fitness > self.best_fitness:  # 更新
            self.best_fitness = temp_fitness
            self.best_position = temp_position
        return particles

    def plot(self, particles):
        for i in range(len(particles)):
            position = particles[i].position
            fitness = self.compute_fitness(position)
            plt.vlines(position,0,fitness,colors='r')
            plt.scatter(position,fitness, c='r')

    def solve(self):
        particles = self.initial()        # 初始化粒子群
        self.plot(particles)
        for _ in range(self.Iternumbers): # 迭代更新
            particles = self.update(particles)
            self.plot(particles)

if __name__ == '__main__':
    pso = PSO(-1, 1, 50, 2, 2, 1, 1, 2)
    X = np.arange(-1,1,0.01)
    Y = pso.compute_fitness(X)
    plt.plot(X, Y, linewidth = 2)
    pso.solve()

    plt.grid()
    plt.show()

粒子群迭代图如下：