粒子群算法求解无约束优化问题Rastrigin函数

一、标准测试函数（Rastrigin函数）

Rastrigin函数是一个典型的非线性多峰函数，在搜索区域内存在许多极大值和极小值，导致寻找全局最小值比较困难，常用来测试寻优算法的性能。Rastrigin函数表达式和函数图像如下：
$$f(x_1,x_2)=20+x_1^2+x_2^2-10\cos (2\pi x_1)-10\cos (2\pi x_2)$$

其中全局最小值点为（0，0），最小值为0

二、使用粒子群算法求解Rastrigin函数全局最小值

用PSO求$\min f(x_1,x_2)$在 x 1 , x 2 ∈ { − 5 , 5 } x_1,x_2 \in \{-5,5\} x1​,x2​∈{−5,5}的最小值

运行算法100次，pso得到的最优解，最优值如下：

最优解是：x=1e-05, y=1e-05 最优值是：0.00000
最优解是：x=2e-05, y=0.0 最优值是：0.00000
最优解是：x=1e-05, y=1e-05 最优值是：0.00000
最优解是：x=0.0, y=0.0 最优值是：0.00000
最优解是：x=-0.0, y=0.0 最优值是：0.00000
最优解是：x=0.0, y=2e-05 最优值是：0.00000
最优解是：x=3e-05, y=1e-05 最优值是：0.00000
最优解是：x=0.0, y=-4e-05 最优值是：0.00000
最优解是：x=3e-05, y=-0.0 最优值是：0.00000
最优解是：x=2e-05, y=0.0 最优值是：0.00000
最优解是：x=1e-05, y=2e-05 最优值是：0.00000
最优解是：x=2e-05, y=-0.0 最优值是：0.00000
最优解是：x=2e-05, y=1e-05 最优值是：0.00000
最优解是：x=-4e-05, y=2e-05 最优值是：0.00000
最优解是：x=-0.0, y=1e-05 最优值是：0.00000
最优解是：x=-2e-05, y=-0.0 最优值是：0.00000
最优解是：x=-0.0, y=-1e-05 最优值是：0.00000
最优解是：x=1e-05, y=0.0 最优值是：0.00000
最优解是：x=-1e-05, y=-2e-05 最优值是：0.00000
最优解是：x=-3e-05, y=1e-05 最优值是：0.00000
最优解是：x=2e-05, y=-2e-05 最优值是：0.00000
最优解是：x=-1e-05, y=-0.0 最优值是：0.00000
最优解是：x=-0.0, y=9e-05 最优值是：0.00000
最优解是：x=-0.0, y=0.0 最优值是：0.00000
最优解是：x=-1e-05, y=2e-05 最优值是：0.00000
最优解是：x=3e-05, y=-2e-05 最优值是：0.00000
最优解是：x=0.0, y=0.0 最优值是：0.00000
最优解是：x=-1e-05, y=1e-05 最优值是：0.00000
最优解是：x=-1e-05, y=1e-05 最优值是：0.00000
最优解是：x=1e-05, y=-1e-05 最优值是：0.00000
最优解是：x=3e-05, y=1e-05 最优值是：0.00000
最优解是：x=-2e-05, y=-2e-05 最优值是：0.00000
最优解是：x=2e-05, y=-4e-05 最优值是：0.00000
最优解是：x=-2e-05, y=-7e-05 最优值是：0.00000
最优解是：x=3e-05, y=4e-05 最优值是：0.00000
最优解是：x=-2e-05, y=2e-05 最优值是：0.00000
最优解是：x=-2e-05, y=-1e-05 最优值是：0.00000
最优解是：x=3e-05, y=-0.0 最优值是：0.00000
最优解是：x=3e-05, y=1e-05 最优值是：0.00000
最优解是：x=3e-05, y=3e-05 最优值是：0.00000
最优解是：x=-1e-05, y=5e-05 最优值是：0.00000
最优解是：x=1e-05, y=2e-05 最优值是：0.00000
最优解是：x=1e-05, y=-0.0 最优值是：0.00000
最优解是：x=-1e-05, y=-2e-05 最优值是：0.00000
最优解是：x=-1e-05, y=0.0 最优值是：0.00000
最优解是：x=-4e-05, y=1e-05 最优值是：0.00000
最优解是：x=2e-05, y=0.0 最优值是：0.00000
最优解是：x=-2e-05, y=-1e-05 最优值是：0.00000
最优解是：x=0.0, y=2e-05 最优值是：0.00000
最优解是：x=-4e-05, y=-0.0 最优值是：0.00000
最优解是：x=1e-05, y=-0.0 最优值是：0.00000
最优解是：x=-1e-05, y=-0.0 最优值是：0.00000
最优解是：x=-3e-05, y=-0.0 最优值是：0.00000
最优解是：x=1e-05, y=-0.0 最优值是：0.00000
最优解是：x=-4e-05, y=1e-05 最优值是：0.00000
最优解是：x=2e-05, y=-1e-05 最优值是：0.00000
最优解是：x=0.0, y=-0.0 最优值是：0.00000
最优解是：x=-1e-05, y=-2e-05 最优值是：0.00000
最优解是：x=0.0, y=2e-05 最优值是：0.00000
最优解是：x=-1e-05, y=-1e-05 最优值是：0.00000
最优解是：x=-1e-05, y=-5e-05 最优值是：0.00000
最优解是：x=-0.0, y=-1e-05 最优值是：0.00000
最优解是：x=-1e-05, y=-3e-05 最优值是：0.00000
最优解是：x=-5e-05, y=1e-05 最优值是：0.00000
最优解是：x=-0.0, y=3e-05 最优值是：0.00000
最优解是：x=1e-05, y=4e-05 最优值是：0.00000
最优解是：x=-1e-05, y=1e-05 最优值是：0.00000
最优解是：x=-2e-05, y=-0.0 最优值是：0.00000
最优解是：x=-0.0, y=0.0 最优值是：0.00000
最优解是：x=-1e-05, y=0.0 最优值是：0.00000
最优解是：x=2e-05, y=-1e-05 最优值是：0.00000
最优解是：x=4e-05, y=-3e-05 最优值是：0.00000
最优解是：x=-0.0, y=-1e-05 最优值是：0.00000
最优解是：x=2e-05, y=-0.0 最优值是：0.00000
最优解是：x=2e-05, y=4e-05 最优值是：0.00000
最优解是：x=0.0, y=1e-05 最优值是：0.00000
最优解是：x=-0.0, y=2e-05 最优值是：0.00000
最优解是：x=1e-05, y=2e-05 最优值是：0.00000
最优解是：x=1e-05, y=-1e-05 最优值是：0.00000
最优解是：x=-2e-05, y=-3e-05 最优值是：0.00000
最优解是：x=-2e-05, y=0.0 最优值是：0.00000
最优解是：x=-1e-05, y=-2e-05 最优值是：0.00000
最优解是：x=0.0, y=1e-05 最优值是：0.00000
最优解是：x=2e-05, y=-3e-05 最优值是：0.00000
最优解是：x=1e-05, y=7e-05 最优值是：0.00000
最优解是：x=2e-05, y=-1e-05 最优值是：0.00000
最优解是：x=1e-05, y=-1e-05 最优值是：0.00000
最优解是：x=-1e-05, y=2e-05 最优值是：0.00000
最优解是：x=2e-05, y=0.0 最优值是：0.00000
最优解是：x=-2e-05, y=-3e-05 最优值是：0.00000
最优解是：x=-0.0, y=1e-05 最优值是：0.00000
最优解是：x=-0.0, y=0.0 最优值是：0.00000
最优解是：x=-0.0, y=-1e-05 最优值是：0.00000
最优解是：x=1e-05, y=-1e-05 最优值是：0.00000
最优解是：x=1e-05, y=2e-05 最优值是：0.00000
最优解是：x=-1e-05, y=5e-05 最优值是：0.00000
最优解是：x=0.0, y=-1e-05 最优值是：0.00000
最优解是：x=7e-05, y=-4e-05 最优值是：0.00000
最优解是：x=-1e-05, y=-4e-05 最优值是：0.00000
最优解是：x=2e-05, y=4e-05 最优值是：0.00000
最优解是：x=-1e-05, y=0.0 最优值是：0.00000

运行算法1次函数收敛图（左图），运行算法100次，目标函数值稳定在0（右图）：

三、完整代码（Python）

from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl

mpl.rcParams['font.sans-serif'] = ['Times New Roman']  # 指定默认字体
mpl.rcParams['axes.unicode_minus'] = False  # 解决保存图像是负号'-'显示为方块的问题


class PSO:
    def __init__(self):
        self.w = 1  # 惯性权重
        self.c1 = 2  # 加速系数
        self.c2 = 2  # 加速系数
        self.dimension = 2
        self.pop_size = 100
        self.max_iteration = 5000
        self.max_velocity = 0.5  # 限制粒子的最大速度为0.5

    def evaluate(self, X):
        x = X[:, 0]
        y = X[:, 1]
        return 2 * 10 + x ** 2 - 10 * np.cos(2 * np.pi * x) + y ** 2 - 10 * np.cos(2 * np.pi * y)

    def update_velocity(self, X, V, pbest, gbest):
        """
           根据速度更新公式更新每个粒子的速度
           :param V: 粒子当前的速度矩阵，20*2 的矩阵
           :param X: 粒子当前的位置矩阵，20*2 的矩阵
           :param pbest: 每个粒子历史最优位置，20*2 的矩阵
           :param gbest: 种群历史最优位置，1*2 的矩阵
        """
        r1 = np.random.random((self.pop_size, 1))
        r2 = np.random.random((self.pop_size, 1))
        V = self.w * V + self.c1 * r1 * (pbest - X) + self.c2 * r2 * (gbest - X)  # 直接对照公式写就好了
        # 防止越界处理, 若越过边界，则退回边界/ 若越过边界，则重新初始化
        # V[V < -self.max_velocity] = -self.max_velocity
        # V[V > self.max_velocity] = self.max_velocity
        # 若越过边界，则重新初始化
        V[V < -self.max_velocity] = np.random.uniform(low=-.2, high=.2)
        V[V > self.max_velocity] = np.random.uniform(low=-.2, high=.2)

        return V

    def update_position(self, X, V):
        X = X + V
        return X

    def run(self):
        trace = []

        X = np.random.uniform(-5, 5, size=(self.pop_size, self.dimension))
        V = np.random.uniform(-0.2, 0.2, size=(self.pop_size, self.dimension))

        pbest = X
        gbest = X[0]

        X_fitness = self.evaluate(X)
        pbest_fitness = X_fitness
        gbest_fitness = X_fitness[0]

        for i in range(0, self.max_iteration):
            V = self.update_velocity(X, V, pbest, gbest)  # 更新速度
            X = self.update_position(X, V)  # 更新位置
            X_fitness = self.evaluate(X)
            best_fitness = np.min(X_fitness)
            # 更新每个粒子的历史最优位置
            for j in range(self.pop_size):
                if X_fitness[j] < pbest_fitness[j]:
                    pbest[j] = X[j]
                    pbest_fitness[j] = X_fitness[j]

            # 更新群体的最优位置
            if best_fitness < gbest_fitness:
                gbest = X[np.argmin(X_fitness)]
                gbest_fitness = best_fitness
            trace.append((gbest, gbest_fitness))

        x, y = gbest
        print(f'x={x:.2f}, y={y:.2f} 全局最小值：{gbest_fitness:.2f}')
        self.plot_objective_value(trace)

    def plot_objective_value(self, trace):
        obj_list = [fitness for gbest, fitness in trace]
        plt.plot(np.arange(self.max_iteration), obj_list, color="#191970", linewidth=2., alpha=1.)
        plt.grid(True)
        plt.show()


if __name__ == '__main__':
    PSO().run()