文献翻译 (3)：粒子群优化算法

1 缘起

  粒子群优化 (Particle swarm optimization, PSO) 算法是一种基于种群的搜索算法，其模拟了鸟群中鸟类的社会行为。粒子群概念的最初意图是用图形模拟鸟群优美但不可预测的编排，发现控制鸟类同步飞行能力的模式，并通过重新组合成最佳编队来改变航向，后衍生出简单却高效的PSO优化算法。
  在PSO中，个体，亦称粒子通过超维搜索空间“飞行”。搜索空间内粒子位置的变化是基于个体模仿其他个体的成功社会心理倾向。因此，群体中粒子的变化受到其邻居的经验或知识的影响。换言之，一个粒子的搜索行为会受到群体内其他粒子同样行为的影响，所以PSO可以称为是一种共生协作算法。对这种社会行为建模的结果是使得粒子随机返回搜索空间中先前成功区域的搜索过程。

2 示意

  PSO擅长于超维空间中寻找预定义函数的最大或最小值。令$f(X)$表示从向量参数$X$生成实值的函数，而$X$可以是空间中的任意取值。例如$X$是平面坐标$(x,y)$，PSO将返回以下函数的最小值：
$$f(x,y)=(x-3.14{)}^2+(y-2.72{)}^2+\sin (3x+1.41)+\sin (4y-1.73)$$  绘制如下：

  代码如下：

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.animation as animation

raw = np.random.rand(100, 3)
fig = plt.figure()

ax = plt.axes(projection="3d")
x = y = np.arange(start=0, stop=5, step=0.02)
X, Y = np.meshgrid(x, y)
Z = (X - 3.14) ** 2 + (Y - 2.72) ** 2 + np.sin(3 * X + 1.41) + np.sin(4 * Y + 1.73)
ax.plot_surface(X, Y, Z, alpha=0.9, cstride=1, rstride=1, cmap='rainbow')

rotate = lambda angle: ax.view_init(azim=angle)
rot_animation = animation.FuncAnimation(fig, rotate, frames=np.arange(0, 362, 2), interval=100)
rot_animation.save('rotation.gif', dpi=80, writer='imagemagick')

  如图，该函数像一个弯曲的鸡蛋盒，并不属于凸优化，因此可能找到局部最小值而非全局最小值。那如何来找到该函数的最小值呢？
  1）计算每一个点？似乎代价昂贵了些；
  2）随机选择样本点？在找到的最小值点附近，可能找到更小值的点。这一事实也是PSO想要做的。
  与鸟儿捕食相似，首先在空间中找到一组随机点，也就是粒子，然后令它们朝着随机的方向寻找最小值点。在每一次迭代中，每个粒子围绕它找到的最小点以及整个粒子群找到的最小点进行搜索。在确定的多次迭代后，这群粒子曾经探索过的最小值便视作函数的最小值。

3 细节

  假设有$P$个粒子，在第$t$次迭代中的第$i$个粒子记作$X^i(t)$。对应于节2中的问题，有$X^i(t)=(x^i(t),y^i(t))$，对应了其在空间中的坐标。进一步，每个粒子的速度定义为$V^i(t)=(v_x^i(t),v_y^i(t))$。在下一次迭代中，当前粒子的坐标将更新为：
$$X^i(t+1)=X^i(t)+V^i(t+1)$$其等价于
$$\begin{gathered} x^i(t+1)=x^i(t)+v_x^i(t+1) \\ {y}^i(t+1)=y^i(t)+v_y^i(t+1) \end{gathered}$$于此同时，速度的更新为：
$$V^i(t+1)=w{V}^i(t)+c_1{r}_1(pbes{t}^i-X^i(t))+c_2{r}_2(gbest-X^i(t))$$其中$r_1$和$r_2$是0和1之间的随机数；$w,c_1,c_2$是PSO的常量参数；$pbes{t}^i$粒子$X^i(t)$取得最小值的坐标；以及$gbest$是所有粒子取得的最小值的坐标。该更新的好处是通过将该粒子找到的最佳点与当前坐标的差值附加在当前速度$V^i(t)$上，可以推动粒子朝着$pbes{t}^i$移动。$gbest$也是同样的道理。
  我们将$w\in (\mathrm{0..1}]$称为惯性权重参数，其决定了粒子保持原有速度的程度。$c_1$和$c_2$分别称为认知和社交系数，它们用于细化粒子本身与群体的搜索结果之间的权衡。
  PSO区别于其他优化算法的一个很重要的属性是不依赖于目标函数的梯度，这使得其可以适应于$f(X)$求导困难的情况。另一个属性是支持并行化，为了找到最优解，多个粒子可以同时行动，收集它们所找到的最小值，我们需要做的则是更新$gbest$。.这使得map-reduce架构完美匹配PSO。

4 实现

  1）定义优化目标且展示如下：

  代码如下：

import numpy as np
import matplotlib.pyplot as plt


def f(x, y):
    """优化目标"""
    return (x - 3.14) ** 2 + (y - 2.72) ** 2 + np.sin(3 * x + 1.41) + np.sin(4 * y - 1.73)


if __name__ == '__main__':
    x, y = np.array(np.meshgrid(np.linspace(0, 5, 100), np.linspace(0, 5, 100)))
    z = f(x, y)
    x_min = x.ravel()[z.argmin()]
    y_min = y.ravel()[z.argmin()]
    plt.figure(figsize=(8, 6))
    plt.imshow(z, extent=[0, 5, 0, 5], origin='lower', cmap='viridis', alpha=0.5)
    plt.colorbar()
    plt.plot([x_min], [y_min], marker='x', markersize=5, color="white")
    contours = plt.contour(x, y, z, 10, colors='black', alpha=0.4)
    plt.clabel(contours, inline=True, fontsize=8, fmt="%.0f")
    plt.show()

  2）随机初始化粒子及其速度：

# 初始化粒子
    np.random.seed(1)
    n_particles = 20
    X = np.random.rand(2, n_particles) * 5
    V = np.random.randn(2, n_particles) * 0.1

  看看它们在哪里：

  代码如下：

import numpy as np
import matplotlib.pyplot as plt


def f(x, y):
    """优化目标"""
    return (x - 3.14) ** 2 + (y - 2.72) ** 2 + np.sin(3 * x + 1.41) + np.sin(4 * y - 1.73)


if __name__ == '__main__':
    # 初始化粒子
    np.random.seed(1)
    n_particles = 20
    X = np.random.rand(2, n_particles) * 5
    V = np.random.randn(2, n_particles) * 0.1

    x, y = np.array(np.meshgrid(np.linspace(0, 5, 100), np.linspace(0, 5, 100)))
    z = f(x, y)
    x_min = x.ravel()[z.argmin()]
    y_min = y.ravel()[z.argmin()]
    plt.figure(figsize=(8, 6))
    plt.imshow(z, extent=[0, 5, 0, 5], origin='lower', cmap='viridis', alpha=0.5)
    plt.colorbar()
    plt.scatter(X[0], X[1])
    plt.plot([x_min], [y_min], marker='x', markersize=5, color="white")
    contours = plt.contour(x, y, z, 10, colors='black', alpha=0.4)
    plt.clabel(contours, inline=True, fontsize=8, fmt="%.0f")
    plt.show()

  3）计算当前的最小值点：

pbest = X
pbest_obj = f(X[0], X[1])
gbest = pbest[:, pbest_obj.argmin()]
gbest_obj = pbest_obj.min()

  加上方向看看：

  代码如下：

import numpy as np
import matplotlib.pyplot as plt


def f(x, y):
    "Objective function"
    return (x - 3.14) ** 2 + (y - 2.72) ** 2 + np.sin(3 * x + 1.41) + np.sin(4 * y - 1.73)


x, y = np.array(np.meshgrid(np.linspace(0, 5, 100), np.linspace(0, 5, 100)))
z = f(x, y)

x_min = x.ravel()[z.argmin()]
y_min = y.ravel()[z.argmin()]

n_particles = 20
np.random.seed(1)
X = np.random.rand(2, n_particles) * 5
V = np.random.randn(2, n_particles) * 0.1

pbest = X
pbest_obj = f(X[0], X[1])
gbest = pbest[:, pbest_obj.argmin()]
gbest_obj = pbest_obj.min()

fig, ax = plt.subplots(figsize=(8, 6))
fig.set_tight_layout(True)
img = ax.imshow(z, extent=[0, 5, 0, 5], origin='lower', cmap='viridis', alpha=0.5)
fig.colorbar(img, ax=ax)
ax.plot([x_min], [y_min], marker='x', markersize=5, color="white")
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)
p_plot = ax.scatter(X[0], X[1], marker='o', color='blue', alpha=0.5)
p_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([0, 5])
ax.set_ylim([0, 5])
plt.show()

  4）初始化相关参数：

c1 = c2 = 0.1
w = 0.8

  5）设置更新函数：

def update():
    global V, X, pbest, pbest_obj, gbest, gbest_obj

    r1, r2 = np.random.rand(2)
    V = w * V + c1 * r1 * (pbest - X) + c2 * r2 * (gbest.reshape(-1, 1) - X)
    X = X + V
    obj = f(X[0], X[1])
    pbest[:, (pbest_obj >= obj)] = X[:, (pbest_obj >= obj)]
    pbest_obj = np.array([pbest_obj, obj]).min(axis=0)
    gbest = pbest[:, pbest_obj.argmin()]
    gbest_obj = pbest_obj.min()

  6）更新一次看看：

  更新两次：

  更新5次：

  更新10次：

  更新20次：

  更新30次：

  代码如下：

import numpy as np
import matplotlib.pyplot as plt


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


x, y = np.array(np.meshgrid(np.linspace(0, 5, 100), np.linspace(0, 5, 100)))
z = f(x, y)

x_min = x.ravel()[z.argmin()]
y_min = y.ravel()[z.argmin()]

c1 = c2 = 0.1
w = 0.8

n_particles = 20
np.random.seed(1)
X = np.random.rand(2, n_particles) * 5
V = np.random.randn(2, n_particles) * 0.1

pbest = X
pbest_obj = f(X[0], X[1])
gbest = pbest[:, pbest_obj.argmin()]
gbest_obj = pbest_obj.min()


def update():
    global V, X, pbest, pbest_obj, gbest, gbest_obj

    r1, r2 = np.random.rand(2)
    V = w * V + c1 * r1 * (pbest - X) + c2 * r2 * (gbest.reshape(-1, 1) - X)
    X = X + V
    obj = f(X[0], X[1])
    pbest[:, (pbest_obj >= obj)] = X[:, (pbest_obj >= obj)]
    pbest_obj = np.array([pbest_obj, obj]).min(axis=0)
    gbest = pbest[:, pbest_obj.argmin()]
    gbest_obj = pbest_obj.min()


# 这里设置更新次数
for i in range(1):
    update()

fig, ax = plt.subplots(figsize=(8, 6))
fig.set_tight_layout(True)
img = ax.imshow(z, extent=[0, 5, 0, 5], origin='lower', cmap='viridis', alpha=0.5)
fig.colorbar(img, ax=ax)
ax.plot([x_min], [y_min], marker='x', markersize=5, color="white")
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)
p_plot = ax.scatter(X[0], X[1], marker='o', color='blue', alpha=0.5)
p_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([0, 5])
ax.set_ylim([0, 5])
plt.show()

  7）画一个动图看看：

  完整代码：

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


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


x, y = np.array(np.meshgrid(np.linspace(0, 5, 100), np.linspace(0, 5, 100)))
z = f(x, y)

x_min = x.ravel()[z.argmin()]
y_min = y.ravel()[z.argmin()]

c1 = c2 = 0.1
w = 0.8

n_particles = 20
np.random.seed(1)
X = np.random.rand(2, n_particles) * 5
V = np.random.randn(2, n_particles) * 0.1

pbest = X
pbest_obj = f(X[0], X[1])
gbest = pbest[:, pbest_obj.argmin()]
gbest_obj = pbest_obj.min()


def update():
    global V, X, pbest, pbest_obj, gbest, gbest_obj

    r1, r2 = np.random.rand(2)
    V = w * V + c1 * r1 * (pbest - X) + c2 * r2 * (gbest.reshape(-1, 1) - X)
    X = X + V
    obj = f(X[0], X[1])
    pbest[:, (pbest_obj >= obj)] = X[:, (pbest_obj >= obj)]
    pbest_obj = np.array([pbest_obj, obj]).min(axis=0)
    gbest = pbest[:, pbest_obj.argmin()]
    gbest_obj = pbest_obj.min()


fig, ax = plt.subplots(figsize=(8, 6))
fig.set_tight_layout(True)
img = ax.imshow(z, extent=[0, 5, 0, 5], origin='lower', cmap='viridis', alpha=0.5)
fig.colorbar(img, ax=ax)
ax.plot([x_min], [y_min], marker='x', markersize=5, color="white")
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)
p_plot = ax.scatter(X[0], X[1], marker='o', color='blue', alpha=0.5)
p_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([0, 5])
ax.set_ylim([0, 5])


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


anim = FuncAnimation(fig, animate, frames=list(range(1, 30)), interval=500, blit=False, repeat=True)
anim.save("PSO.gif", dpi=120, writer="imagemagick")

print("PSO found best solution at f({})={}".format(gbest, gbest_obj))
print("Global optimal at f({})={}".format([x_min, y_min], f(x_min, y_min)))

参考文献

【1】Complete Step-by-step Particle Swarm Optimization Algorithm from Scratch
【2】https://machinelearningmastery.com/a-gentle-introduction-to-particle-swarm-optimization/