PSO粒子群搜索算法步骤及应用实例（一）

案例一：利用PSO求解无约束连续函数优化问题

代码：
PSO_1.py

'''
#案例一：
# 求解Rastrigin函数的极小值，自变量取值范围为[-5,5]
# 函数：fun = 2*a+x**2-a*np.cos(2*np.pi*x) + y**2-a*np.cos(2*np.pi*y)
'''

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
# 绘制三维图像和色图
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D

mpl.rcParams['font.sans-serif'] = ['SimHei']#指定默认字体，可显示中文
mpl.rcParams['axes.unicode_minus'] = False#正常显示图像中的负号


def func(x, y, a=10):
    '''计算目标函数的值'''
    return 2*a+x**2-a*np.cos(2*np.pi*x) + y**2-a*np.cos(2*np.pi*y)



def show_picture():
    '''绘制函数图像'''
    x = np.arange(-5, 5, 0.1)  # [-5,5]
    y = np.arange(-5, 5, 0.1)  # [-5,5]
    x, y = np.meshgrid(x, y)  # 生成网格

    fig = plt.figure()
    ax = Axes3D(fig)
    _ = ax.plot_surface(x, y, func(x, y), cmap=cm.coolwarm)
    # return ax
    plt.show()

# 定义适应度函数，速度更新函数，位置更新函数
def fitness_func(x):
    '''通过矩阵的方式计算粒子的适应度值
    输入：X为 size * 2 维的矩阵，其中size为粒子种群大小
    输出：1*size 维的适应度值向量'''

    a = 10
    return func(x[:, 0], x[:, 1], a)


def velocity_update(v, x, p_best, g_best, c1, c2, w, max_val):
    '''更新每个粒子的速度向量
    输入：
    v-当前粒子群体的速度向量矩阵
    x-当前粒子群体的位置向量矩阵
    p_best-每个粒子历史最优的位置向量矩阵
    g_best-历史最优的粒子位置向量
    c1,c2-“自我认知”和“社会认知”系数
    w-惯性系数
    max_val-限制粒子的最大速度
    输出：
    更新后的粒子群体速度向量矩阵
    '''
    size = x.shape[0]  # 获取粒子个数
    r1 = np.random.random((size, 1))
    r2 = np.random.random((size, 1))
    v_update = w*v + c1*r1*(p_best-x) + c2*r2*(g_best-x)
    #检测并修改超出自变量取值范围的值
    v_update[v_update<-1*max_val] = -1*max_val
    v_update[v_update>max_val] = max_val

    return v_update

def position_update(x,v):
    '''更新粒子群体的位置矩阵'''
    return x + v

#粒子群算法主函数
def pso_solver():
    '''粒子群算法求解主函数'''
    #parameters
    w = 1
    c1 = 2
    c2 = 2
    dim = 2#自变量个数
    size = 20 #种群大小
    iter_num = 1000#算法迭代最大次数
    max_val = 0.5#限制粒子的最大速度为0.5
    fitness_value_list = []#记录搜索过程中全体最优适应度的变化
    fitness_value_position_list = []#记录搜索过程中最优粒子位置的变化
    #PSO算法迭代过程
    x = np.random.uniform(-5,5,size=(size,dim))
    v = np.random.uniform(-1*max_val,max_val,size=(size,dim))
    #每个粒子的局部最优和群体的历史最优
    p_best = x
    p_fitness = fitness_func(p_best)
    g_fitness = min(p_fitness)
    g_best = x[p_fitness.argmin()]#.argmin()--->返回最小值所在的索引下标

    fitness_value_list.append(g_fitness)
    fitness_value_position_list.append(g_best)

    #迭代
    for i in range(1,iter_num):
        v = velocity_update(v,x,p_best,g_best,c1,c2,w,max_val)
        x = position_update(x,v)

        p_fitness_current = fitness_func(x)
        g_fitness_current = p_fitness_current.min() #这一代粒子的最优值
        g_best_current = x[p_fitness_current.argmin()] #这一代粒子的最优位置

        #更新每个粒子的历史最优位置
        for j in range(size):
            if p_fitness[j]>p_fitness_current[j]:
                p_best[j] = x[j] #更新该粒子的历史最优位置
                p_fitness[j] = p_fitness_current[j]
        #更新种群历史最优位置
        if g_fitness > g_fitness_current:
            g_fitness = g_fitness_current
            g_best = g_best_current

            # 记录最优值的改变
            fitness_value_list.append(g_fitness)
            fitness_value_position_list.append(g_best)
        
        i += 1
    #返回历史最优值变化列表和最优值对应粒子位置向量列表
    return fitness_value_list,fitness_value_position_list

if __name__ == "__main__":
    # 绘制函数图像
    show_picture()
    fitness_value_list,fitness_value_position_list = pso_solver()
    print("搜索结果的最优值为:",fitness_value_list[-1])
    print("最优的解为:",fitness_value_position_list[-1])
    
    #绘制迭代过程图
    plt.figure()
    plt.plot(range(len(fitness_value_list)),fitness_value_list)
    plt.title("PSO迭代中最优值的变化")
    plt.show()

    fig = plt.figure()
    ax = Axes3D(fig)
    x = [pos[0] for pos in fitness_value_position_list]
    y = [pos[1] for pos in fitness_value_position_list]
    x_mesh = np.arange(min(x), max(x), 0.01)  
    y_mesh = np.arange(min(y), max(y), 0.01)  
    x_mesh, y_mesh = np.meshgrid(x_mesh, y_mesh)  # 生成网格
    _ = ax.plot_surface(x_mesh, y_mesh, func(x_mesh, y_mesh), cmap=cm.coolwarm)

    # ax.plot3D(x,y,fitness_value_list,c='r')
    ax.scatter(x,y,fitness_value_list,c='r',marker='H',linewidths=1.5)
    plt.show()
    

结果：

案例二：利用PSO求解有约束连续函数优化问题
若约束条件不满足，通常在目标函数中引入惩罚项，迫使模型在迭代计算的过程中向可行域内寻优
代码：
PSO_2.py

'''案例二：PSO求解有约束连续优化问题
问题：min fun = 2*a+x**2-a*np.cos(2*np.pi*x) + y**2-a*np.cos(2*np.pi*y)
      s.t. x + y <= 6
           3*x-2*y <=5

'''

import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
# 绘制三维图像和色图
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
from PSO_1 import *

mpl.rcParams['font.sans-serif'] = ['SimHei']  # 指定默认字体，可显示中文
mpl.rcParams['axes.unicode_minus'] = False  # 正常显示图像中的负号


# 计算约束惩罚项


def calc_e1(x):
    '''x+y<6约束惩罚项
    输入x为单个粒子位置矩阵'''
    e = x[0] + x[1] - 6
    return max(0, e)


def calc_e2(x):
    '''3*x-2*y<=5约束惩罚项
    输入x为单个粒子位置矩阵'''
    e = 3*x[0] - 2*x[1] - 5
    return max(0, e)


def calc_Lj(e1, e2):
    '''计算两个约束的归一化权重'''
    total_sum = e1.sum()+e2.sum()
    if total_sum == 0:
        # 防止分母为0
        return 0, 0
    else:
        return e1.sum()/total_sum, e2.sum()/total_sum

# 定义每个粒子历史最优位置和种群历史最优位置更新函数，需对约束进行处理


def update_pbest(pbest, pbest_fitness, pbest_e, xi, xi_fitness, xi_e):
    '''判断是否需要跟新粒子的历史最优位置
    PS:此处针对的都是单个粒子,而不是粒子群体
    pbest-该粒子当前历史最优位置 1*2维向量
    pbest_fitness-历史最优位置对应的目标函数值 1*1
    pbest_e--历史最优位置对应的约束惩罚项 1*1
    xi-当前位置 1*2维向量
    xi_fitness-当前位置的适应度函数值1*1
    xi_e-当前位置的约束惩罚项1*1'''

    if pbest_e <= 1e-7 and xi_e <= 1e-7:
        # 该粒子历史最优位置和当前位置都没有违反约束
        if xi_fitness < pbest_fitness:
            return xi, xi_fitness, xi_e
        else:
            return pbest, pbest_fitness, pbest_e
    elif pbest_e > 1e-7 and xi_e < 1e-7:
        # 该粒子历史最优位置违反约束，返回当前位置
        return xi, xi_fitness, xi_e
    elif pbest_e < 1e-7 and xi_e > 1e-7:
        # 该粒子历史最优位置没有违反约束，返回历史最优位置（不变）
        return pbest, pbest_fitness, pbest_e
    elif pbest_e > 1e-7 and xi_e > 1e-7:
        # 该粒子历史最优位置和当前位置都违反约束，选择适应度跟优的返回
        if pbest_fitness < xi_fitness:
            return pbest, pbest_fitness, pbest_e
        else:
            return xi, xi_fitness, xi_e


def update_gbest(gbest, gbest_fitness, gbest_e, pbest, pbest_fitness, pbest_e):
    '''更新全局最优位置
    gbest-种群历史最优位置-gbest_fitness-种群历史最优适应度值-gbest_e种群历史最优位置的约束惩罚项(单个)
    pbest-当前粒子群体的位置矩阵,pbest_fitness-粒子群体的适应度值向量,pbest_e-粒子群体的约束惩罚项向量(粒子种群)
    '''
    # 合并矩阵info:(各列的含义)
    #     0.种群每个粒子的历史最优位置对应的适应度
    #     1.历史最优位置对应的惩罚项
    #     2.当前适应度
    #     3.当前目标函数
    #     4.约束一的惩罚项
    #     5.约束二的惩罚项
    #     6.惩罚项加权和

    pbest2 = np.concatenate(
        [pbest, pbest_fitness.reshape(-1, 1), pbest_e.reshape(-1, 1)], axis=1)
    # 找出所有没有违反约束的粒子个体,1e-7是考虑到计算机的精度设置,值等同于0
    pbest2_1 = pbest2[pbest2[:, -1] < 1e-7]
    if len(pbest2_1) > 0:  # 存在没有违反约束的粒子
        index = pbest2_1[:, 2].argmin()
        pbesti, pbesti_fitness, pbesti_e = pbest2_1[index,
                                                    :2], pbest2_1[index, 2], pbest2_1[index, -1]
    else:  # 所有粒子都违反了约束
        index = pbest2[:, 2].argmin()
        pbesti, pbesti_fitness, pbesti_e = pbest2[index,
                                                  :2], pbest2[index, 2], pbest2[index, -1]
    # 种群当前最优解和种群历史全局最优的比较
    best, best_fitness, best_e = update_pbest(
        gbest, gbest_fitness, gbest_e, pbesti, pbesti_fitness, pbesti_e)
    return best, best_fitness, best_e


def pso_solver2():
    '''pso主函数'''
    # PSO parameters
    w = 1
    c1 = 2
    c2 = 2
    dim = 2  # 自变量个数
    size = 20  # 种群大小
    iter_num = 1000  # 算法迭代最大次数
    max_val = 0.5  # 限制粒子的最大速度为0.5
    fitness_value_list = []  # 记录搜索过程中全体最优适应度的变化
    fitness_value_position_list = []  # 记录搜索过程中最优粒子位置的变化
    info = np.zeros((size, 7))
    # 随机初始化位置矩阵和速度矩阵
    x = np.random.uniform(-5, 5, size=(size, dim))
    v = np.random.uniform(-1*max_val, max_val, size=(size, dim))
    # 每个粒子的局部最优和群体的历史最优
    p_best = x
    # 计算每个粒子的适应度
    for i in range(size):
      #   info[i, 3] = fitness_func(x[i, :])  # 目标函数值
        info[i, 3] = func(x[i][0], x[i][1])  # 目标函数值
        info[i, 4] = calc_e1(x[i, :])  # 约束1惩罚项
        info[i, 5] = calc_e2(x[i, :])  # 约束2惩罚项
    # 计算惩罚项权重及适应度
    L1, L2 = calc_Lj(info[:, 4], info[:, 5])  # 计算惩罚项权重
    info[:, 6] = L1*info[:, 4] + L2*info[:, 5]  # 惩罚项的加权求和
    # info[:,2] = info[:,3] + L1*info[:,4] + L2*info[:,5]#适应度值
    info[:, 2] = info[:, 3]+info[:, 6]

    # 初始化历史最优
    info[:, 0] = info[:, 2]  # 最优位置适应度值
    info[:, 1] = info[:, 6]  # 最优位置对应的惩罚项

    gbest = x[info[:, 0].argmin()]  # 全局最优粒子的位置
    gbest_fitness = info[info[:, 0].argmin(), 0]  # 全局最优粒子对应的适应度值
    gbest_e = info[info[:, 0].argmin(), 1]  # 全局最优粒子对应的惩罚项

    fitness_value_list.append(gbest_fitness)
    fitness_value_position_list.append(gbest)

    # 迭代过程
    for _ in range(iter_num):
        v = velocity_update(v, x, p_best, gbest, c1, c2, w, max_val)
        x = position_update(x, v)
        # 计算每个粒子的目标函数和约束惩罚项
        # 计算每个粒子的适应度
        for i in range(size):
            # info[i, 3] = fitness_func(x[i, :])  # 目标函数值
            info[i, 3] = func(x[i][0], x[i][1])  # 目标函数值
            info[i, 4] = calc_e1(x[i, :])  # 约束1惩罚项
            info[i, 5] = calc_e2(x[i, :])  # 约束2惩罚项
        # 计算惩罚项权重及适应度
        L1, L2 = calc_Lj(info[:, 4], info[:, 5])  # 计算惩罚项权重
        info[:, 6] = L1*info[:, 4] + L2*info[:, 5]  # 惩罚项的加权求和
        # info[:,2] = info[:,3] + L1*info[:,4] + L2*info[:,5]#适应度值
        info[:, 2] = info[:, 3]+info[:, 6]

        # 更新每个粒子的历史最优位置
        for i in range(size):
            pbesti, pbesti_fitness, pbesti_e = update_pbest(
                p_best[i], info[i, 0], info[i, 1], x[i], info[i, 2], info[i, -1])
            p_best[i] = pbesti
            info[i, 0] = pbesti_fitness
            info[i, 1] = pbesti_e

        # 跟新全局最优
        gbest, gbest_fitness, gbest_e = update_gbest(
            gbest, gbest_fitness, gbest_e, x, info[:, 2], info[:, -1])
        # 记录种群最优值的变化
        fitness_value_list.append(gbest_fitness)
        fitness_value_position_list.append(gbest)

    return fitness_value_list, fitness_value_position_list


if __name__ == "__main__":
    fitness_value_list, fitness_value_position_list = pso_solver2()
    print("PSO搜索的最优值为:", fitness_value_list[-1])
    print("最优值为:", fitness_value_position_list[-1])

    # 绘制迭代过程图
    plt.figure()
    plt.plot(range(len(fitness_value_list[:30])), fitness_value_list[:30])
    plt.title("PSO迭代中最优值的变化")
    plt.show()

    fig = plt.figure()
    ax = Axes3D(fig)
    x = [pos[0] for pos in fitness_value_position_list]
    y = [pos[1] for pos in fitness_value_position_list]
    x_mesh = np.arange(min(x), max(x), 0.01)
    y_mesh = np.arange(min(y), max(y), 0.01)
    x_mesh, y_mesh = np.meshgrid(x_mesh, y_mesh)  # 生成网格
    _ = ax.plot_surface(x_mesh, y_mesh, func(x_mesh, y_mesh), cmap=cm.coolwarm)

    # ax.plot3D(x,y,fitness_value_list,c='r')
    ax.scatter(x, y, fitness_value_list, c='r', marker='H', linewidths=1.5)
    plt.show()

案例三：利用PSO求解TSP问题（离散优化）
未完待续~~~~~~~

参考资料

1.《智能优化方法》_汪定伟
2.《Python最优化算法实战》_苏振裕