智能优化算法学习-粒子群算法（Particle Swarm Optimization ...

yukamu

一、简介

     粒子群算法（Particle Swarm Optimization， PSO）是一个经典的智能优化算法，属于进化计算技术（evolutionary computation）的一种（CS领域习惯这么定义）。PSO的思想源于对鸟群捕食行为的研究，模拟鸟群飞行觅食的行为。在捕食的时候，一个鸟群的个体直接按通过集体的协作使群体行为达到最优目的，这是一种基于群体智能的优化方法。
     PSO没有遗传算法的交叉和变异操作，只通过追随当前搜索到的最优值来找寻全局最优（学习率  的设置关键，太大容易陷入局部最优，太小收敛速度不理想）。
情景设想：
     一群鸟正在觅食，在远处有一片玉米地，但是所有的鸟都不知道其具体的位置，只知道自己当前的位置与玉米的距离。所以搜寻当前距离玉米地最近的鸟的周围区域即是找到玉米地的最简单有效的策略。
<hr/>二、PSO算法原理

基本核心：利用群体中的个体对信息的共享，使整个群体的运动在问题求解空间中产生从无序到有序的演化过程，从而获得问题的最优解。
       我们把上面情景中的鸟称之为“粒子”，那么我们所要优化的问题即是找到玉米地。在算法中所有的粒子都具有一个位置向量（在解空间中的位置）与速度向量（决定下次迭代的方向和速度or步长），并根据目标函数来计算当前所在位置的适应度值（可以理解为到玉米地的距离）。在迭代过程中，粒子除了根据自身的“经验”（历史位置）进行学习（个人学习因子  ）,还可以根据种群中最优粒子的“经验”来学习(（社会学习因子  ）)，从而确定下一次迭代时需要如何调整和改变迭代的方向和速度。如此迭代n次后，最终使得种群的粒子逐步趋于最优解。



图-1.PSO示意图

    在一个空间n中，以 表示第i个粒子的位置向量，以 表示第i个粒子的速度向量，对于粒子i在维度d的速度向量迭代公式为：


位置迭代公式：


其中：

表示历史最佳位置

表示群体的历史最佳位置

表示个体学习因子

表示社会学习因子（种群）

[0, 1]范围内的随机数，以增加搜索的随机性
     在上面的迭代公式可以看出，种群中粒子通过不断地向自身和种群的历史信息进行学习，从而找到问题的最优解。但是在后续的研究中，发现PSO的一些迭代更新的问题，即 更新具有随机性，也即这样迭代虽然使PSO算法全局优化能力很强，但是局部搜索能力较差。实际上在算法迭代初期需要PSO有较强的全局优化能力，而在算法的后期整个种群需要更强的局部搜索能力。后来学者通过引用惯性权重，提出了PSO的惯性权重模型（个人理解，类似机器学习中的学习率）：


      一般来说  设置为0.8-1.2之前，算法有较快的收敛速度；但是超过1.2时，容易陷入局部最优。所以在搜索过程中可以对  进行动态调整，一开始赋予较大的正值，随着搜索的进行，可以线性地使其逐渐减小，这样可以保证算法在刚开始粒子有较大的速度步长，并在全局范围内搜索到较好的区域；而在后期，较小的  可以保证粒子能在极值点周围做精细搜索，从而使算法有较大的概率向全局最优收敛。（想到Lagrangian Relaxation的乘子与步长的迭代 ）

的动态迭代公式：


其中 表示最大迭代次数

表示最小惯性权重，大多数情况下设置为0.4

表示最大惯性权重，大多数情况下设置为0.9
t表示当前迭代次数



图-2.PSO算法流程图

终止迭代的条件有两种：①超过最大的迭代次数 ②最优解的误差小于设置的界限
<hr/>三、实例1-粒子群求解无约束优化问题-求解Rastrigin函数

Rastrigin函数是一个典型的非线性多峰函数，在搜索区域内存在许多极大值和极小值，导致寻找全局最优值比较困难，常常用来测试寻优算法的性能。其函数表达式如下：

（这里我们令A=10）
可以看出其全局最小值点为（0，0）



图-3.Rastrigin函数图像

基础设置：
# Work hard everyday!!!
# Student:Sum3 Klaus
# Time: 2022/3/12 14:35
# Email: Sum3rebirth@gmail.com

# 速度
# Vi+1 = w * Vi + c1*r1*(Pbest_i - Xi) + c2*r1*(Gbest_i - Xi)
# 位置
# Xi+1 = Xi + Vi+1
# vi 与 xi 分表表示粒子在第i维的速度和位置
# phest_i 与 gbest_i分别表示 某单个粒子在第i维的最佳历史位置，种群的历史最佳位置

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

plt.rcParams['font.sans-serif'] = ['SimHei']
plt.rcParams['axes.unicode_minus'] = False定义适应度值函数fitness_func：
def fitness_func(X):
    """计算粒子的适应度值，X 的维度是 size * 2"""
    A = 10
    pi = np.pi
    x = X[:, 0]
    y = X[:, 1]
    return 2*A + x**2 - A*np.cos(2*pi*x) + y**2 - A*np.cos(2*pi*y)定义速度迭代函数velocity_update：
def velocity_update(V, X, pbest, gbest, c1, c2, w, max_val):
    """
    根据惯性权重速度更新公式，迭代每个粒子的速度
    :param V: 粒子当前的速度矩阵，shape = 20*2
    :param X: 粒子当前的位置矩阵，shape = 20*2
    :param pbest:每个粒子历史最优位置，shape = 20*2
    :param gbest:种群的历史最优位置， shape = 1*2
    :param c1:
    :param c2:
    :param w:
    :param max_val:粒子的最大速度限制
    :return:
    """
    size = X.shape[0]
    r1 = random.random((size, 1))
    r2 = random.random((size, 1))
    V = w*V + c1 * r1 * (pbest * X) + c2 * r2 * (gbest * X)
    # 边界设置，防止越界
    V[V < -max_val] = -max_val
    V[V > max_val] = max_val
    return V定义位置迭代函数position_upodate：
def position_upodate(X, V):
    """
    更新粒子的位置
    :param X: 粒子当前的位置矩阵，shape = 20*2
    :param V: 粒子当前的速度矩阵，shape = 20*2
    :return:
    """
    return X + V主函数PSO：
def pso():
    # pso 的参数
    w = None
    w_min = 0.4
    w_max = 0.9
    c1 = 2
    c2 = 2
    r1 = None                    # (0, 1)的随机数
    r2 = None
    dim = 2                      # 变量的个数（维度）
    size = 20                    # 种群的大小
    iter_num_max = 1000          # 最大迭代次数
    max_val = 0.5                # 粒子的最大速度限制
    best_fitness = float(9e10)   # 初始的适应度值，在迭代过程中不断减小这个值
    fitness_value_lst = []       # 记录每次迭代的种群适应度值

    X = random.uniform(-5, 5, size=(size, dim))
    V = random.uniform(-0.5, 0.5, size=(size, dim))

    p_fitness_values = fitness_func(X)
    g_fitness_best_values = p_fitness_values.min()
    fitness_value_lst.append(g_fitness_best_values)
    pbest = X
    gbest = X[p_fitness_values.argmin()]

    # 开始迭代
    for i in range(1, iter_num_max):
        w = w_max - (w_max - w_min)*i/iter_num_max
        V = velocity_update(V, X, pbest, gbest, c1, c2, w, max_val)
        X = position_upodate(X, V)

        p_fitness_values_new = fitness_func(X)
        g_fitness_best_values_new = p_fitness_values_new.min()

        # 更新每个粒子的历史最优位置
        for j in range(size):
            if p_fitness_values[j] > p_fitness_values_new[j]:
                pbest[j] = X[j]
                p_fitness_values[j] = p_fitness_values_new[j]

        # 更新群体的最优位置
        if g_fitness_best_values > g_fitness_best_values_new:
            gbest = X[p_fitness_values_new.argmin()]
            g_fitness_best_values = g_fitness_best_values_new

        # 记录最优迭代结果
        fitness_value_lst.append(g_fitness_best_values)
        i += 1
    # 输出最优迭代结果
    print(f'最优值是： {fitness_value_lst[-1]}')
    print(f'最优解是： x = {gbest[0]} y = {gbest[1]}')输出结果：



图-4.PSO迭代结果



图-5.PSO算法目标函数值的变化曲线

欢迎各位老师和同学批评指正！
<hr/>

《python最优化算法实战》

xiangtingsl

非常实用！收藏！