PSO算法

最近实验的涉及到PSO算法(Particle Swarm Optimization)，pso属于进化算法(pso也属于群智算法)。由Kennedy和Eberhart提出，模拟鸟群觅食进行优化。基本核心是利用群体中的个体对信息的共享从而使整个群体的运动在问题求解空间中产生从无序到有序的演化过程，从而获得问题的最优解。

假设一个粒子\(\)\(\vec x_i\)\(\)属于D维空间的向量，群体进行搜索食物，粒子开始行动，以速度\(\)\(\vec v_i\)\(\)飞出，进行搜索，整个过程粒子对比先前搜索到的最优点\(\)\(\vec pbest_i\)\(\)，这个是粒子个体的搜索经验，放大到群体范围，粒子还可以参考全体粒子的最优值(全体的搜索经验，也就是粒子间的搜索协作和共享)\(\)\(\vec gbest\)\(\)

粒子速度和位置的更新:
\(\)\(\vec v_i = w \times \vec v_i + c_1 \times rand() \times (pbest_i – \vec x_i) + c_2 \times rand() \times (gbest – \vec x_i)\)\(\)
\(\)\(\vec x_i = \vec x_i + \vec v_i\)\(\)

算法过程：
1.随机初始化种群中粒子的位置x与速度v
for 搜索迭代循环：
for 粒子循环：
2. 更新粒子的位置x与速度v。
3. 通过约束函数f，进行计算当前位置的函数值fitness，
4. 如果当前粒子x找到的值，比前面已找到的值pbest小，则更新
5. 如果前面找到的值pbest小于全局最小值gbest，则更新
end for
end for

使用PSO寻找函数的最小值：
\(\)\(f(x) = \frac{sin \sqrt{x^2_1 +x^2_2}}{\sqrt{x^2_1 +x^2_1}} + e^{\frac{con 2 \pi x_1 +con 2 \pi x_2}{2}} -2.71289\)\(\)

• 

• 

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

def f(x):
    return np.sin(np.sqrt(np.add.reduce(np.power(x, 2))))/np.sqrt(np.add.reduce(np.power(x, 2))) + \
           np.exp(np.add.reduce(np.cos(2*np.pi*x))/2) - 2.71289

X = np.arange(-2, 2, 0.05)
Y = np.arange(-2, 2, 0.05)
X, Y = np.meshgrid(X, Y)

Z = f(np.array([X, Y]))

fig = plt.figure()

# 立体图
# ax = Axes3D(fig)
# plt.title("PSO")
# ax.plot_surface(X, Y, Z, rstride=1, cstride=1, color='b', )
# ax.set_xlabel('x label', color='r')
# ax.set_ylabel('y label', color='g')
# ax.set_zlabel('z label', color='b')

# 等高线图
cs = plt.contourf(X, Y, Z)
fig.colorbar(cs)
plt.show()

通过PSO算法搜索上面函数的最小值以及搜索过程的动态展示。通过调整参数w会发现，w值越大，全局搜索能力更强，w值越小，全局搜索能力较弱。

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


class PSO(object):

    def __init__(self, func, n_partciles=30, dims=2, iters=200, w=0.9, lr=(0.5, 1.5), rangep=(-2, 2), rangev=(-0.5, 0.5)):
        self.func = func
        self.n_partciles = n_partciles  # 种群大小
        self.dims = dims  # 粒子维度
        self.iters = iters  # 迭代次数
        self.w = w
        self.lr = lr  # 个体粒子的搜索经验，全体的搜索经验
        self.rangep = rangep  # 粒子搜索的限制范围
        self.rangev = rangev  # 粒子速度的限制范围

    def initSwarm(self):
        # 初始化粒子位置与速度，计算初始化下的函数适应度
        particles = np.zeros((self.n_partciles, self.dims))
        v = np.zeros((self.n_partciles, self.dims))
        fitness = np.zeros(self.n_partciles)

        for i in range(self.n_partciles):
            particles[i] = [(np.random.rand() - 0.5) * self.rangep[1] for i in range(self.dims)]
            v[i] = [(np.random.rand() - 0.5) * self.rangep[1] for i in range(self.dims)]
            fitness[i] = self.func(particles[i])

        # 初始化pbest, gbest
        gbest, gbestfitness = particles[fitness.argmin()].copy(), fitness.min()
        pbest, pbestfitness = particles.copy(), fitness.copy()

        return particles, v, fitness, gbest, gbestfitness, pbest, pbestfitness

    def run(self):
        particles, v, fitness, gbest, gbestfitness, pbest, pbestfitness = self.initSwarm()
        # 迭代循环
        result = np.zeros(self.iters)
        fig = plt.figure()
        plt.ion()
        X = np.arange(-2, 2, 0.05)
        Y = np.arange(-2, 2, 0.05)
        X, Y = np.meshgrid(X, Y)
        Z = f(np.array([X, Y]))
        cs = plt.contourf(X, Y, Z)
        fig.colorbar(cs)
        for i in range(self.iters):
            # 速度更新
            for j in range(self.n_partciles):
                v[j] = self.w * v[j] + self.lr[0] * np.random.rand() * (pbest[j] - particles[j]) \
                       + self.lr[1] * np.random.rand() * (gbest - particles[j])
            v[v < self.rangev[0]] = self.rangev[0]
            v[v > self.rangev[1]] = self.rangev[1]
            # 位置更新
            for j in range(self.n_partciles):
                particles[j] += v[j]
            particles[particles < self.rangep[0]] = self.rangep[0]
            particles[particles > self.rangep[1]] = self.rangep[1]

            s = plt.scatter(particles[:, 0], particles[:, 1], c='orange')
            plt.pause(0.4)
            s.remove()

            # 适应度更新
            for j in range(self.n_partciles):
                fitness[j] = self.func(particles[j])

            for j in range(self.n_partciles):
                if fitness[j] < pbestfitness[j]:
                    pbestfitness[j] = fitness[j]
                    pbest[j] = particles[j].copy()

            if pbestfitness.min() < gbestfitness:
                gbestfitness = pbestfitness.min()
                gbest = particles[pbestfitness.argmin()].copy()
            result[i] = gbestfitness

        plt.ioff()
        # 显示图形
        plt.scatter(gbest[0], gbest[1], c='white')
        plt.show()
        return result


def f(x):
    return np.sin(np.sqrt(np.add.reduce(np.power(x, 2))))/np.sqrt(np.add.reduce(np.power(x, 2))) + \
           np.exp(np.add.reduce(np.cos(2*np.pi*x))/2) - 2.71289

pso = PSO(f)
print(pso.run())

PSO算法的思想十分简洁，但是粒子有时候逃不脱经验的束缚，粒子会被群体最优值拉回来。从而得不到全局最优解，就像其他优化算法一样，例如最简单的梯度下降。当局部最优解很多时，很难寻到最优解。

reference
1. Poli R, Kennedy J, Blackwell T. Particle swarm optimization[J]. Swarm intelligence, 2007, 1(1): 33-57.