什么是机器学习中的PSO？

2025年2月3日 | 阅读10分钟

粒子群优化 (PSO) 是机器学习和其他领域中用于优化问题的一种计算技术，它模仿鸟群或鱼群的社会行为。在这种技术中，一组候选解，“粒子”，在搜索空间中移动以找到最优解。每个粒子根据其自身和邻居的过往经验来改进其在搜索空间中的位置，并在迭代中收敛到最优解。这种方法对于优化传统优化技术可能不足的复杂多维问题非常有用。PSO 是一种非常简单的算法，易于实现，并应用于神经网络训练、特征选择和机器学习模型的超参数调优。它在机器学习优化相关的各种问题中的应用，得益于其自适应性和鲁棒性。

PSO技术是这样的：每个粒子都有一个内存，其中存储着整个种群的全局最佳位置，以及它自己过去最好的位置。它可能会根据存储在每个个体（粒子）内存中的数据和其速度vi来更新其位置xi(t+1)，如下图所示。

在PSO技术中，X通常表示种群。可以使用随机数生成器按照以下方程初始化此种群：

其中 L 和 U 分别表示问题搜索空间域的下界和上界。X 中的每个个体/粒子 (xi) 的速度 (vi) 可以使用以下方程更新：

此处，

w: 这是一个惯性权重，用于提高收敛速度。

c1 & c2:: 它们代表加速度的系数。

t: 当前时间（迭代次数）。

之后，粒子xi的位置更新为：

此PSO模型的下一步是适应度函数，它将计算每个粒子的质量。然后，找出当前粒子 P 到目前为止已达到的最佳位置。这也被称为最佳个人位置和全局最佳解 P。上述操作将一直持续到满足停止条件。

PSO的最终步骤如下面的伪代码所示：

参数初始化：N - 粒子数量，d - 维度，L - 下边界，U - 上边界
随机解初始化：X，在搜索空间内
随机速度初始化：X 中每个解的 vi
while (不满足终止条件)
for 每个粒子 in X do
计算当前粒子的目标函数值
更新当前粒子的最佳个人位置，xp
更新全局最佳位置（xg）
使用公式 (2) 更新速度
使用公式 (3) 更新当前粒子的位置
end for
end while
返回全局最佳解，xg

现在我们将应用PSO及其基准测试。

导入库

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.animation import FuncAnimation
import math

基准函数

需要使用从测试函数（如 A_Rosenbrock、A_Rastgrin、A_Ackley、A_Sphere、A_Beale 等）开发的一组基准函数来解决此优化问题。我们打算在开发的基准函数上实现PSO。

def create_benchmark(function):
  
  if function == "A_Rosenbrock":
    a=(-10,10)
    x = np.linspace(-10,10,100)
    y = np.linspace(-10,10,100)
    x,y = np.meshgrid(x,y)
    def funcd(x,y): return 100*((x-y**2))**2+(x-1)**2
    return 100*((x-y**2))**2+(x-1)**2,x,y,a,funcd

  if function == "A_Rastgrin":
    a=(-5.12,5.12)
    x = np.linspace(-5.12, 5.12, 100)     
    y = np.linspace(-5.12, 5.12, 100) 
    x,y = np.meshgrid(x,y)
    def funcd(x,y) : return (x**2 - 10 * np.cos(2 * np.pi * x))+(y**2 - 10 * np.cos(2 * np.pi * y))
    return (x**2 - 10 * np.cos(2 * np.pi * x))+(y**2 - 10 * np.cos(2 * np.pi * y)) + 20, x,y,a,funcd  
  
  if function == "A_Ackley":
    a=(-32.768, 32.768)
    minimum_r, maximum_r = -32.768, 32.768
    xaxis = np.arange(minimum_r, maximum_r, 2.0)
    yaxis = np.arange(minimum_r, maximum_r, 2.0)
    x, y = np.meshgrid(xaxis, yaxis)
    def funcd(x,y) : return -20.0 * np.exp(-0.2 * np.sqrt(0.5 * (x**2 + y**2))) - np.exp(0.5 * (np.cos(2 * np.pi * x) + np.cos(2 * np.pi * y))) + np.e + 20
    return -20.0 * np.exp(-0.2 * np.sqrt(0.5 * (x**2 + y**2))) - np.exp(0.5 * (np.cos(2 * np.pi * x) + np.cos(2 * np.pi * y))) + np.e + 20,x,y,a,funcd  
  
  if function == "A_Sphere":
    a=(-100,100)
    x, y = np.array(np.meshgrid(np.linspace(-100,100,100), np.linspace(-100,100,100)))
    def funcd(x,y) : return (4/3)*math.pi*x**3
    return (4/3)*math.pi*x**3, x,y,a,funcd  
  
  if function == "A_Beale":
    a=(-4.5,4.5)
    x, y = np.array(np.meshgrid(np.linspace(-4.5,4.5,100), np.linspace(-4.5,4.5,100)))
    def funcd(x,y) : return (1.5 - x + x*y)**2 + (2.25 - x + x*y**2)**2 + (2.625 - x + x*y**3)**2
    return (1.5 - x + x*y)**2 + (2.25 - x + x*y**2)**2 + (2.625 - x + x*y**3)**2,x,y,a,funcd
  
  if function == "A_Goldstein-Price":
    a=(-2,2)
    x, y = np.array(np.meshgrid(np.linspace(-2,2,100), np.linspace(-2,2,100)))
    def funcd(x,y) : return (1+(x+y+1)**2*(19-14*x+3*x**2-14*y+6*x*y+3*y**2))*(30+(2*x-3*y)**2*(18-32*x+12*x**2+48*y-36*x*y+27*y**2))
    return (1+(x+y+1)**2*(19-14*x+3*x**2-14*y+6*x*y+3*y**2))*(30+(2*x-3*y)**2*(18-32*x+12*x**2+48*y-36*x*y+27*y**2)),x,y,a,funcd
  
  if function == "A_Booth":
    a=(-10,10)
    x, y = np.array(np.meshgrid(np.linspace(-10,10,100), np.linspace(-10,10,100)))
    def funcd(x,y) : return (x+(2*y)-7)**2+(((2*x)+y+5)**2)
    return (x+(2*y)-7)**2+(((2*x)+y+5)**2),x,y,a,funcd

  if function == "A_Bukin":
    a=(-15,-5)
    x, y = np.array(np.meshgrid(np.linspace(-15,-5,100), np.linspace(-3,3,100)))
    def funcd(x,y) : return 100*np.sqrt(np.abs(y-(0.01*(x**2))))+np.abs(x*10)*0.01
    return 100*np.sqrt(np.abs(y-(0.01*(x**2))))+np.abs(x*10)*0.01,x,y,a,funcd

  if function == "A_Matyas":
    a=(-10,10)
    x, y = np.array(np.meshgrid(np.linspace(-10,10,100), np.linspace(-10,10,100)))
    def funcd(x,y) : return 0.26*(x**2+y**2)-(0.48*x*y)
    return 0.26*(x**2+y**2)-(0.48*x*y),x,y,a,funcd

  if function == "A_Levi":
    a=(-10,10)
    x, y = np.array(np.meshgrid(np.linspace(-10,10,100), np.linspace(-10,10,100)))
    def funcd(x,y) : return np.sin(3*np.pi*x)**2 + (x-1)**2*(1+np.sin(3*np.pi*y)*np.sin(3*np.pi*y))+ (y-1)*(y-1)*(1+np.sin(2*np.pi*y)*np.sin(2*np.pi*y))
    return np.sin(3*np.pi*x)**2 + (x-1)**2*(1+np.sin(3*np.pi*y)*np.sin(3*np.pi*y))+ (y-1)*(y-1)*(1+np.sin(2*np.pi*y)*np.sin(2*np.pi*y)),x,y,a,funcd

  if function == "A_Himmelblau":
    a=(-5,5)
    x, y = np.array(np.meshgrid(np.linspace(-5,5,100), np.linspace(-5,5,100)))
    def funcd(x,y) : return (((x**2+y-11)**2) + (((x+y**2-7)**2)))
    return (((x**2+y-11)**2) + (((x+y**2-7)**2))),x,y,a,funcd

  if function == "A_Three hump camel":
    a=(-5,5)
    x, y = np.array(np.meshgrid(np.linspace(-5,5,100), np.linspace(-5,5,100)))
    def funcd(x,y) : return ((2*(x**2)-1.05*(x**4)+((x**6)/6)+(x*y)+(y**2)))
    return ((2*(x**2)-1.05*(x**4)+((x**6)/6)+(x*y)+(y**2))),x,y,a,funcd

  if function == "A_Easom":
    a=(-100,100)
    x, y = np.array(np.meshgrid(np.linspace(-100,100,100), np.linspace(-100,100,100)))
    def funcd(x,y) : return -(np.cos(x)*np.cos(y))*np.exp(-(x-math.pi)**2-(y-math.pi)**2)
    return -(np.cos(x)*np.cos(y))*np.exp(-(x-math.pi)**2-(y-math.pi)**2),x,y,a,funcd

  if function == "A_Cross in tray":
    c=(-10,10)
    x, y = np.array(np.meshgrid(np.linspace(-10,10,100), np.linspace(-10,10,100)))
    a=np.fabs(100-np.sqrt(x*x+y*y)/np.pi)
    b=np.fabs(np.sin(x)*np.sin(y)*np.exp(a))+1
    def funcd(x,y) : return -0.0001*(np.fabs(np.sin(x)*np.sin(y)*np.exp(np.fabs(100-np.sqrt(x*x+y*y)/np.pi)))+1)**0.1
    return funcd(x,y),x,y,c,funcd

  if function == "A_Egg holder":
    c=(-512,512)
    x, y = np.array(np.meshgrid(np.linspace(-512,512,100), np.linspace(-512,512,100)))
    a=np.sqrt(np.fabs(y+x/2+47))
    b=np.sqrt(np.fabs(x-(y+47)))
    def funcd(x,y) : return -(y+47)*np.sin(np.sqrt(np.fabs(y+x/2+47)))-x*np.sin(np.sqrt(np.fabs(x-(y+47))))
    return -(y+47)*np.sin(a)-x*np.sin(b),x,y,c,funcd

  if function == "A_McCormick":
    a=(-1.5,4)
    x, y = np.array(np.meshgrid(np.linspace(-1.5,4,100), np.linspace(-3,4,100)))
    def funcd(x,y) : return np.sin(x+y)+(x-y)**2-(1.5*x)+(2.5*y)+1
    return np.sin(x+y)+(x-y)**2-(1.5*x)+(2.5*y)+1,x,y,a,funcd

  if function == "A_Schaffer":
    a=(-50,50)
    x = np.linspace(-50,50)
    y = np.linspace(-50,50)
    x,y = np.meshgrid(x,y)
    def funcd(x,y) : return 0.5+(np.sin(x**2-y**2)-0.5)**2/((1+0.001*(x**2+y**2)))**2
    return 0.5+(np.sin(x**2-y**2)-0.5)**2/((1+0.001*(x**2+y**2)))**2,x,y,a,funcd

  if function == "A_Stiblinski-tag":
    a=(-5,5)
    x = np.arange(-5, 5, 0.25)
    y = np.arange(-5, 5, 0.25)
    x, y = np.meshgrid(x, y)
    def funcd(x,y) : return 0.5 * ((x**4 + y**4) - 16 * (x**2 + y**2)+ 5 * (x + y))
    return 0.5 * ((x**4 + y**4) - 16 * (x**2 + y**2)+ 5 * (x + y)),x,y,a,funcd

  if function == "A_Salomon":
    a=(-10,10)
    x = np.arange(-100, 100)
    y = np.arange(-10, 10)
    x, y = np.meshgrid(x, y)
    r = x**2+y**2
    def funcd(x,y) : return 1-np.cos(2 * np.pi * np.sqrt(x**2+y**2) + 0.1 * np.sqrt(x**2+y**2))
    return 1-np.cos(2 * np.pi * np.sqrt(r) + 0.1 * np.sqrt(r)),x,y,a,funcd

  if function == "A_Objective":
    a = (0,5)
    x, y = np.array(np.meshgrid(np.linspace(0,5,100), np.linspace(0,5,100)))

    def funcd(x,y) : return (x-3.14)**2 + (y-2.72)**2 + np.sin(3*x+1.41) + np.sin(4*y-1.73)
    return (x-3.14)**2 + (y-2.72)**2 + np.sin(3*x+1.41) + np.sin(4*y-1.73), x,y,a,funcd

  if function == "A_Parabolic":
    a = (0,5)
    x, y = np.array(np.meshgrid(np.linspace(0,5,100), np.linspace(0,5,100)))
    def funcd(x,y) : return x**2+y**2
    return x**2+y**2, x,y,a,funcd

def update():

    "Function to do one iteration of particle swarm optimization"
    global V, X, pbest, pbest_obj, gbest, gbest_obj
    # Update params
    r1, r2 = np.random.rand(2)
    V = w * V + c1*r1*(ptbest - X) + c2*r2*(gbest.reshape(-1,1)-X)
    X = X + V
    obj = funcd(X[0], X[1])
    pbest[:, (ptbest_obj>= obj)] = X[:, (ptbest_obj>= obj)]
    ptbest_obj= np.array([pbest_obj, obj]).min(axis=0)
    gtbest = pbest[:, pbest_obj.argmin()]
    gtbest_objective= pbest_obj.min()    

def animate(i):
    "Steps of PSO: algorithm update and show in a plot"
    title = 'Iteration {:02d}'.format(i)
    # Update params
    update()
    # Set picture
    ax.set_title(title)
    pbest_plot.set_offsets(pbest.T)
    pt_plot.set_offsets(X.T)
    pt_arrow.set_offsets(X.T)
    pt_arrow.set_UVC(V[0], V[1])
    gbest_plot.set_offsets(gbest.reshape(1,-1))
    return ax, pbest_plot, pt_plot, pt_arrow, gbest_plot

以下是一些标准基准函数，通常用于评估优化方法的性能。这些函数通常用于基准测试优化算法的性能评估。

List_of_functions= ["A_Rosenbrock", "A_Rastgrin",
                    "A_Ackley","A_Sphere","A_Beale",
                    "A_Goldstein-Price","A_Booth",
                    "A_Bukin","A_Matyas","A_Levi",
                    "A_Himmelblau","A_Three hump camel",
                    "A_Easom","A_Cross in tray","A_Egg holder",
                    "A_McCormick","A_Schaffer","A_Stiblinski-tag",
                    "A_Salomon","A_Objective","A_Parabolic"] 

# It's possible that some features, including A_Salomon and A_Bukin, aren't functioning properly in the current version. To carry out appropriate optimization tests, it is essential to make sure that each function's implementation faithfully follows its mathematical formulation.                             

虽然“A_Rastgrin”是首选，但我们选择了“A_Objective”函数作为我们优化过程的基准。这是应用于优化问题的两个不同的基准函数。

选择测试

# Now, we will select the test from the number of tests that are to be used for benchmarking.
# We have applied here the Object function for our purpose.

z,x,y,a,funcd = create_benchmark("A_Objective") 
# If you wish to use ?A_Rastgrin? you can use the underlying code.
# z, x, y, a, funcd = create_benchmark("A_Rastgrin")

我们将尝试找到由基准函数 'x'、'y' 和 'z' 表示的全局最小值和最大值点。现在，由于我们可以找到最高和最低点并与其他点进行比较，因此基准测试变得非常容易，因为在视觉上我们只需要查看图中更接近最小值的点。

# Finding the minimum and maximum value for global

#minimum value of x
minimum_x = x.ravel()[z.argmin()] 

# minimum value of y
minimum_y = y.ravel()[z.argmin()]

# maximum value of x
maximum_x= x.ravel()[z.argmax()]

# maximum value of y
maximum_y = y.ravel()[z.argmax()]

绘制等高线图

我们想生成一个二维图，其中我们将根据颜色图表示函数值。基准函数 'z' 将表示为图像。

它还将叠加等高线，显示函数的等高线，并在先前确定的全局最小值点设置一个红色的“x”。

plt.figure(figsize=(8,6))
plt.imshow(z, extent=[a[0],a[1],a[0],a[1]], origin='lower', cmap='viridis', alpha=0.5)
plt.colorbar()
plt.plot([minimum_x], [minimum_y], marker='x', markersize=10, color="red")
contours = plt.contour(x, y, z, 10, colors='black', alpha=0.4)
plt.clabel(contours, inline=True, fontsize=8, fmt="%.0f")
plt.show()

输出

上面是二维图，其中红色的“X”表示基准的最小值，图中的黑线表示函数的特定等高线值。这种图的颜色表示将帮助您理解值的密度。

参数和初始化

现在我们将设置控制PSO算法行为和收敛的参数。

# Hyper-parameter of the algorithm
c1 = 0.8
c2 = 0.8
w = 0.6

我们需要初始化PSO的粒子。

# Create the  particles
number_of_particles = 40
np.random.seed(100)
# Values of particles
X = np.random.rand(2, number_of_particles) * 5
# Set random velocity for each particle
V = np.random.randn(2, number_of_particles) * 0.1

# Initialize data
ptbest = X
ptbest_obj= funcd(X[0], X[1])
gtbest = pbest[:, pbest_obj.argmin()]
gtbest_objective= pbest_obj.min()

绘制点

我们将实现的粒子群优化算法发展可视化的基础将是上面绘制的基准函数的等高线图。稍后，我们将在该等高线图本身上绘制点。

# Set up base figure: The contour map
fig, ax = plt.subplots(figsize=(8,6))
fig.set_tight_layout(True)

img = ax.imshow(z, extent=[a[0],a[1],a[0],a[1]], origin='lower', cmap='viridis', alpha=0.5)
fig.colorbar(img, ax=ax)
ax.plot([minimum_x], [minimum_y], marker='x', markersize=50, color="red")
contours = ax.contour(x, y, z, 10, colors='black', alpha=0.4)
ax.clabel(contours, inline=True, fontsize=8, fmt="%.0f")
pbest_plot = ax.scatter(pbest[0], pbest[1], marker='o', color='black', alpha=0.5)
pt_plot = ax.scatter(X[0], X[1], marker='o', color='blue', alpha=0.5)
pt_arrow = ax.quiver(X[0], X[1], V[0], V[1], color='blue', width=0.005, angles='xy', scale_units='xy', scale=1)
gbest_plot = plt.scatter([gbest[0]], [gbest[1]], marker='*', s=100, color='black', alpha=0.4)
ax.set_xlim([a[0],a[1]])
ax.set_ylim([a[0],a[1]])

输出

大多数粒子已经从等高图上的最小点散开。这是第一次迭代，请在迭代后再次检查分数。

print("PSO found best solution at funcd({})={}".format(gbest, gbest_obj))
print("Global optimal at funcd({})={}".format([minimum_x,minimum_y], 
funcd(minimum_x,minimum_y)))
# Set iteration of update function and frame interval

输出

PSO解决方案

在 x = 2.99421688 和 y = 1.38382531 时，PSO 发现了解决方案。
此时，目标函数的值为 0.3668918708631955。

全局最优解

对于基准函数，x = 3.1818181818181817 和 y = 3.131313131313131 是已知的全局最优解。
A_Objective 函数现在的值为 -1.8082706615747688。

生成的解似乎不完全是全局最优解。尽管如此，PSO方法找到了一个非常接近全局最优解的解。PSO的目标是找到或足够接近全局最优解，并且在某些情况下它可能会收敛到局部最优解。

注意：重要的是要记住，许多因素会影响PSO成功找到全局最优解的程度，包括超参数的选择、优化景观的复杂性以及粒子初始化。通过改变PSO技术、调整超参数或多次运行过程，可以提高达到全局最优解的可能性。

PSO 迭代

现在，让我们直观地看一下 PSO 迭代。

iteration = list(range(1,50))
anim = FuncAnimation(fig, animate, frames=iteration, interval=100, blit=False, repeat=True)
anim
from IPython.display import HTML
HTML(anim.to_html5_video())
# If It does not run properly, It is recommended to re-run all the cells.

输出

因此可以推断，粒子在迭代过程中向红色的“X”或最小值移动，这表明该方法得到了优化，并且更接近最小值的粒子会产生有利的优化信号。

# the benchmark scores
print("PSO found best solution at funcd({})={}".format(gbest, gbest_obj))
print("Global optimal at funcd({})={}".format([minimum_x,minimum_y], funcd(minimum_x,minimum_y)))

输出

尽管不完全精确，但PSO算法会收敛到一个非常接近全局最优值的值。尽管粒子群优化算法的目标是在搜索空间内找到全局最优解，但在某些情况下它可能会收敛到局部最优解。PSO的优化在每次迭代中都在改进。

下一主题机器学习特征工程

← 上一个下一个 →

什么是机器学习中的PSO？

导入库

基准函数

选择测试

绘制等高线图

参数和初始化

绘制点

注意：重要的是要记住，许多因素会影响PSO成功找到全局最优解的程度，包括超参数的选择、优化景观的复杂性以及粒子初始化。通过改变PSO技术、调整超参数或多次运行过程，可以提高达到全局最优解的可能性。

PSO 迭代

联系信息

关注我们

教程

面试题

在线编译器

Python

Java

.Net Framework

AI, ML and Data Science

Cloud Technology

B.Tech and MCA

Web Technology

PHP

Software Testing

Technical Interview

Java Interview

Python

Web Interview

Database Interview

B.Tech / MCA

Important Interview

Software Testing Interview

Company Interviews

Online Compilers

Multiple Choice Questions

机器学习

监督式学习

分类

杂项

相关教程

面试题

什么是机器学习中的PSO？

导入库

基准函数

选择测试

绘制等高线图

参数和初始化

绘制点

注意：重要的是要记住，许多因素会影响PSO成功找到全局最优解的程度，包括超参数的选择、优化景观的复杂性以及粒子初始化。通过改变PSO技术、调整超参数或多次运行过程，可以提高达到全局最优解的可能性。

PSO 迭代

相关帖子

心理模型对数据科学家和 ML/AI 从业者的重要性

机器学习的 A/B 测试

CatBoost 与 XGBoost

覆盖网络

蛙跳算法

机器学习在信号处理中的应用

预测建模与机器学习

机器学习中的数据可视化

CycleGAN

将 Scrapy 爬虫部署到 ScrapingHub

订阅 Tpoint Tech

联系信息

关注我们

教程

面试题

在线编译器