粒子群算法(PSO)的Python实现（求解多元函数的极值）

电玩女神 2022-12-07 01:50 199阅读 0赞

PSO算法算是寻优算法中比较简单的一种，其大概思想是：  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L0N5cmlsX0tJ_size_16_color_FFFFFF_t_70_pic_center]  
现在我们计算：  
![在这里插入图片描述][20200914215615556.gif_pic_center]  
的最大值，每一个变量的取值范围都是(1，25)。  
Python代码为：

# -*- coding: utf-8 -*-
    """ @Time ： 2020/9/13 10:08 @Author ：KI @File ：pso.py @Motto：Hungry And Humble """
    import math
    import random
    import numpy as np
    import matplotlib.pyplot as plt
    import pylab as mpl
    mpl.rcParams['font.sans-serif'] = ['SimHei']
    
    
    class PSO:
        def __init__(self, dimension, time, size, low, up, v_low, v_high):
            # 初始化
            self.dimension = dimension  # 变量个数
            self.time = time  # 迭代的代数
            self.size = size  # 种群大小
            self.bound = []  # 变量的约束范围
            self.bound.append(low)
            self.bound.append(up)
            self.v_low = v_low
            self.v_high = v_high
            self.x = np.zeros((self.size, self.dimension))  # 所有粒子的位置
            self.v = np.zeros((self.size, self.dimension))  # 所有粒子的速度
            self.p_best = np.zeros((self.size, self.dimension))  # 每个粒子最优的位置
            self.g_best = np.zeros((1, self.dimension))[0]  # 全局最优的位置
    
            # 初始化第0代初始全局最优解
            temp = -1000000
            for i in range(self.size):
                for j in range(self.dimension):
                    self.x[i][j] = random.uniform(self.bound[0][j], self.bound[1][j])
                    self.v[i][j] = random.uniform(self.v_low, self.v_high)
                self.p_best[i] = self.x[i]  # 储存最优的个体
                fit = self.fitness(self.p_best[i])
                # 做出修改
                if fit > temp:
                    self.g_best = self.p_best[i]
                    temp = fit
    
        def fitness(self, x):
            """ 个体适应值计算 """
            x1 = x[0]
            x2 = x[1]
            x3 = x[2]
            x4 = x[3]
            x5 = x[4]
            y = math.floor((x2 * np.exp(x1) + x3 * np.sin(x2) + x4 + x5) * 100) / 100
            # print(y)
            return y
    
        def update(self, size):
            c1 = 2.0  # 学习因子
            c2 = 2.0
            w = 0.8  # 自身权重因子
            for i in range(size):
                # 更新速度(核心公式)
                self.v[i] = w * self.v[i] + c1 * random.uniform(0, 1) * (
                        self.p_best[i] - self.x[i]) + c2 * random.uniform(0, 1) * (self.g_best - self.x[i])
                # 速度限制
                for j in range(self.dimension):
                    if self.v[i][j] < self.v_low:
                        self.v[i][j] = self.v_low
                    if self.v[i][j] > self.v_high:
                        self.v[i][j] = self.v_high
    
                # 更新位置
                self.x[i] = self.x[i] + self.v[i]
                # 位置限制
                for j in range(self.dimension):
                    if self.x[i][j] < self.bound[0][j]:
                        self.x[i][j] = self.bound[0][j]
                    if self.x[i][j] > self.bound[1][j]:
                        self.x[i][j] = self.bound[1][j]
                # 更新p_best和g_best
                if self.fitness(self.x[i]) > self.fitness(self.p_best[i]):
                    self.p_best[i] = self.x[i]
                if self.fitness(self.x[i]) > self.fitness(self.g_best):
                    self.g_best = self.x[i]
    
        def pso(self):
            best = []
            self.final_best = np.array([1, 2, 3, 4, 5])
            for gen in range(self.time):
                self.update(self.size)
                if self.fitness(self.g_best) > self.fitness(self.final_best):
                    self.final_best = self.g_best.copy()
                print('当前最佳位置：{}'.format(self.final_best))
                temp = self.fitness(self.final_best)
                print('当前的最佳适应度：{}'.format(temp))
                best.append(temp)
            t = [i for i in range(self.time)]
            plt.figure()
            plt.plot(t, best, color='red', marker='.', ms=15)
            plt.rcParams['axes.unicode_minus'] = False
            plt.margins(0)
            plt.xlabel(u"迭代次数")  # X轴标签
            plt.ylabel(u"适应度")  # Y轴标签
            plt.title(u"迭代过程")  # 标题
            plt.show()
    
    
    if __name__ == '__main__':
        time = 50
        size = 100
        dimension = 5
        v_low = -1
        v_high = 1
        low = [1, 1, 1, 1, 1]
        up = [25, 25, 25, 25, 25]
        pso = PSO(dimension, time, size, low, up, v_low, v_high)
        pso.pso()

运行结果：  
![在这里插入图片描述][20200914215831386.png_pic_center]  
收敛过程：  
![在这里插入图片描述][watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L0N5cmlsX0tJ_size_16_color_FFFFFF_t_70_pic_center 1]  
可以看出，不到10次就收敛了。

matlab代码：

z=@(x)-(x(2)*exp(x(1))+x(3)*sin(x(2))+x(4)*x(5));
    x0=[1;1;1;1;1];
    [x,feval]=fmincon(z,x0,[],[],[],[],[1;1;1;1;1],[25;25;25;25;25])

运行结果：

x =
    
       25.0000
       25.0000
       13.1400
        1.0002
        1.0002
    
    
    feval =
    
      -1.8001e+12

[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L0N5cmlsX0tJ_size_16_color_FFFFFF_t_70_pic_center]: /images/20221123/9036496cfc004a5c9934418093b45f46.png
[20200914215615556.gif_pic_center]: /images/20221123/1592d76ff4904ff1bcedaf6f6843e91a.png
[20200914215831386.png_pic_center]: /images/20221123/501ee56fc823471ba69052a9f0b7bf19.png
[watermark_type_ZmFuZ3poZW5naGVpdGk_shadow_10_text_aHR0cHM6Ly9ibG9nLmNzZG4ubmV0L0N5cmlsX0tJ_size_16_color_FFFFFF_t_70_pic_center 1]: /images/20221123/2c0171dd46a34ae5880cf1db565e4b0b.png