2020.06.24:修改错误的Schwefel

---

0. 前言:

本篇主要是纪录Particle swarm optimization (PSO). A tutorial的阅读心得。

2020/06/20 提供较完整的实验数据；修改错误的PSO代码；提供测试代码(含测试函数)

1. 动机:

最近脑子动到想用PSO来训练神经网路，所以蒐集了一些相关资料，预计将分享5篇文章

Particle swarm optimization (PSO). A tutorialAdaptive Particle Swarm OptimizationOrthogonal Learning Particle Swarm OptimizationComparison of particle swarm optimization and backpropagation as training algorithms for neural networksA hybrid particle swarm optimization–back-propagation algorithm for feedforward neural network training

1.1. 关于PSO:

第一篇是敲门砖，主要是描述基本算法框架，并且进一步的提供超参数的设定指南
第二篇是我硕班用的演算法，我只能说这个真的很猛
第三篇的作者群和第二篇相似，趁这次机会强迫自己看一下

1.2. 关于PSO+NN:

第四篇是很早期的PSO+NN，虽然内容简单但胜在可读性绝佳
第五篇我认为是近代PSO+NN的圣经，但我今天测试的结果是不太好，找机会在分享

其实PSO+NN的文献大多是[变种PSO]+[FFNN]这种组合，主要篇幅都在介绍变种PSO(狼群算法、蜜蜂算法)架构，很少提到怎么与NN结合(如何餵值等)

2. 正文

2.1. 什么是PSO

完全没碰过的请直接看PSO 粒子群演算法
简单来说，一群个被称为粒子的潜在解在多维度解空间中找寻最佳解位置，粒子每次移动都会参考自身过往曾找到的的最佳解位置、所有例子的过往最佳解位置，然后再决定移动方向还有距离。

2.2. PSO的计算流程与算式

假设解空间总共有D个维度，则粒子i的结构应设计为

因为这个演算法明子叫作粒子'群'，所以我们会将N个粒子洒在解空间中，所以

编程上，我习惯用row代表粒子，而column代表维度，用Sphere Function作为範例:

因为需要求解的自变数有D=3个，用N=5个粒子下去求解，因此X矩阵的维度会是5x3，每一个粒子都代表一组可能解所以每一次迭代都会有N=5组可能解彼此竞争

2.2.1 对N个粒子作初始化

a. 初始化位置

b. 取得每一粒子的，因为当前t=0，所以
c. 计算所有粒子的适应值，若条件成立

代表粒子j的适应值为群体中最佳的(假设适应值越大越好)，所以

2.2.2 除非满足停止条件(loss曲线连续停滞n代、loss已满足门槛、迭代次数满足)，否则持续执行下列迴圈(a~e=>a~e=>...=>a~e)

a. 更新每一粒子的速率

这边还会限制速率範围，约束条件我放在2.3.5b. 更新每一粒子的位置

并且确认没有超出解空间

c. 计算每一粒子的适应值
d. 如果

也就是指本次适应值刷新个人过往最佳纪录，则
e. 如果

也就是指本次适应值刷新群体过往最佳纪录，则

2.2.3 最终输出等于g

2.3 PSO的参数指南

这边我直接上结果，若想进一步了解请看原文

2.3.1. 速度爆炸(velocity explosion)

从前文2.2.2(b)可知，x(t+1)=x(t)+v(t+1)，若v(t+1)过大会导致:

x(t+1)容易冲出解空间，最终粒子求得的解会是不能使用(超出範围)，这个靠约束条件可以避免一直贴在墙壁上(因为约束条件)，最终粒子求得的解与实际最佳解差距甚远(因为一直在解空间边缘徘徊)移动幅度过大导致无法收敛，最终粒子求得的解与实际最佳解差距甚远(因为一直大幅度移动的在解空间中)

2.3.2. 加速度常数c

2.3.3 乱数R

numpy.random.uniform

2.3.4 初始位置X(0)与初始速度V(0)

U是指uniform distribution(均匀分布)

为防止速度爆炸(velocity explosion)因此初始速度V(0)建议极小值或者直接设0

2.3.5 速度v的约束条件-避免速度爆炸

k是超参数

2.3.6 另一个速率更新公式-避免速度爆炸

这个就是目前PSO相关文献最常看到的速率更新公式，透过w来缩小v(t)，以防止速度爆炸

当w>1代表全域探勘(exploration)当w<1代表区域开发(exploitation)
原文整理出常见的6个w制定策略，我只列出其中三个定值
乱数
线性递减

2.3.7 群大小

原文建议50个粒子，因为当粒子数大于50后对于求解帮助明显下降

2.4 实现

2.4.1 PSO代码

参考于kuhess/pso-ann文献上并没有明确注明R1和R2的维度，但可以从内文推测是1xD或者Nx1(前者感觉比较合理)；然而GITHUB上的写法通常都是NxD

import numpy as npimport matplotlib.pyplot as pltnp.random.seed(42)class PSO():    def __init__(self, fit_func, num_dim, num_particle=20, max_iter=500,                 x_max=1, x_min=-1, w_max=0.9, w_min=0.4, c1=2.0, c2=2.0, k=0.2):        self.fit_func = fit_func                self.num_dim = num_dim        self.num_particle = num_particle        self.max_iter = max_iter        self.x_max = x_max        self.x_min = x_min        self.w = w_max                self.w_max = w_max        self.w_min = w_min        self.c1 = c1        self.c2 = c2        self.k = k        self._iter = 1        self.gBest_curve = np.zeros(self.max_iter)        self.X = np.random.uniform(low=self.x_min, high=self.x_max,                                   size=[self.num_particle, self.num_dim])        self.V = np.zeros(shape=[self.num_particle, self.num_dim])        self.v_max = self.k*(self.x_max-self.x_min)        self.pBest_X = self.X.copy()        self.pBest_score = self.fit_func(self.X)        self.gBest_X = self.pBest_X[self.pBest_score.argmin()]        self.gBest_score = self.pBest_score.min()        self.gBest_curve[0] = self.gBest_score.copy()    def opt(self):        while(self._iter<self.max_iter):               R1 = np.random.uniform(size=(self.num_particle, self.num_dim))            R2 = np.random.uniform(size=(self.num_particle, self.num_dim))            w = self.w_max - self._iter*(self.w_max-self.w_min)/self.max_iter            for i in range(self.num_particle):                                self.V[i, :] = w * self.V[i, :] \                        + self.c1 * (self.pBest_X[i, :] - self.X[i, :]) * R1[i, :] \                        + self.c2 * (self.gBest_X - self.X[i, :]) * R2[i, :]                               self.V[i, self.v_max < self.V[i, :]] = self.v_max[self.v_max < self.V[i, :]]                self.V[i, -self.v_max > self.V[i, :]] = -self.v_max[-self.v_max > self.V[i, :]]                                self.X[i, :] = self.X[i, :] + self.V[i, :]                self.X[i, self.x_max < self.X[i, :]] = self.x_max[self.x_max < self.X[i, :]]                self.X[i, self.x_min > self.X[i, :]] = self.x_min[self.x_min > self.X[i, :]]                                        score = self.fit_func(self.X[i, :])                if score<self.pBest_score[i]:                    self.pBest_score[i] = score.copy()                    self.pBest_X[i, :] = self.X[i, :].copy()                    if score<self.gBest_score:                        self.gBest_score = score.copy()                        self.gBest_X = self.X[i, :].copy()            self.gBest_curve[i] = self.gBest_score.copy()                     self._iter = self._iter + 1     def plot_curve(self):        plt.figure()        plt.title('loss curve ['+str(round(self.gBest_curve[-1], 3))+']')        plt.plot(self.gBest_curve, label='loss')        plt.grid()        plt.legend()        plt.show()    

2.4.2 测试代码

测试代码

Benchmark可以在这个网站找或者这个网站，两个网站对于同一个测试函数的命名会不太一样!!测试函数是引用于第二篇(APSO)，我们这系列文章都会用该篇结果作为baseline

from PSO import PSOimport numpy as npimport timeimport pandas as pdnp.random.seed(42)def Sphere(x):    if x.ndim==1:        x = x.reshape(1, -1)    return np.sum(x**2, axis=1)def Schwefel_P222(x):    if x.ndim==1:        x = x.reshape(1, -1)    return np.sum(np.abs(x), axis=1) + np.prod(np.abs(x), axis=1)def Quadric(x):    if x.ndim==1:        x = x.reshape(1, -1)        outer = 0    for i in range(x.shape[1]):        inner = np.sum(x[:, :i+1], axis=1)**2        outer = outer + inner        return outerdef Rosenbrock(x):    if x.ndim==1:        x = x.reshape(1, -1)         left = x[:, :-1].copy()    right = x[:, 1:].copy()        return np.sum(100*(right - left**2)**2 + (left-1)**2, axis=1)def Step(x):    if x.ndim==1:        x = x.reshape(1, -1)    return np.sum(np.round((x+0.5), 0)**2, axis=1)def Quadric_Noise(x):    if x.ndim==1:        x = x.reshape(1, -1)            matrix = np.arange(x.shape[1])+1         return np.sum((x**4)*matrix, axis=1)+np.random.rand(x.shape[0])def Schwefel(x):    if x.ndim==1:        x = x.reshape(1, -1)                 return -1*np.sum(x*np.sin(np.abs(x)**.5), axis=1)def Rastrigin(x):    if x.ndim==1:        x = x.reshape(1, -1)         return np.sum(x**2 - 10*np.cos(2*np.pi*x) + 10, axis=1)def Noncontinuous_Rastrigin(x):    if x.ndim==1:        x = x.reshape(1, -1)        outlier = np.abs(x)>=0.5    x[outlier] = np.round(2*x[outlier])/2        return np.sum(x**2 - 10*np.cos(2*np.pi*x) + 10, axis=1)def Ackley(x):    if x.ndim==1:        x = x.reshape(1, -1)        left = 20*np.exp(-0.2*(np.sum(x**2, axis=1)/x.shape[1])**.5)    right = np.exp(np.sum(np.cos(2*np.pi*x), axis=1)/x.shape[1])        return -left - right + 20 + np.edef Griewank(x):    if x.ndim==1:        x = x.reshape(1, -1)    left = np.sum(x**2, axis=1)/4000    right = np.prod( np.cos(x/((np.arange(x.shape[1])+1)**.5)), axis=1)    return left - right + 1d = 30g = 500p = 20times = 1table = np.zeros((2, 11))for i in range(times):    x_max = 100*np.ones(d)    x_min = -100*np.ones(d)    optimizer = PSO(fit_func=Sphere, num_dim=d, num_particle=p, max_iter=g, x_max=x_max, x_min=x_min)    start = time.time()    optimizer.opt()    end = time.time()    table[0, 0] += optimizer.gBest_score    table[1, 0] += end - start        x_max = 10*np.ones(d)    x_min = -10*np.ones(d)    optimizer = PSO(fit_func=Schwefel_P222, num_dim=d, num_particle=p, max_iter=g, x_max=x_max, x_min=x_min)    start = time.time()    optimizer.opt()    end = time.time()    table[0, 1] += optimizer.gBest_score    table[1, 1] += end - start        x_max = 100*np.ones(d)    x_min = -100*np.ones(d)    optimizer = PSO(fit_func=Quadric, num_dim=d, num_particle=p, max_iter=g, x_max=x_max, x_min=x_min)    start = time.time()    optimizer.opt()    end = time.time()    table[0, 2] += optimizer.gBest_score    table[1, 2] += end - start        x_max = 10*np.ones(d)    x_min = -10*np.ones(d)    optimizer = PSO(fit_func=Rosenbrock, num_dim=d, num_particle=p, max_iter=g, x_max=x_max, x_min=x_min)    start = time.time()    optimizer.opt()    end = time.time()    table[0, 3] += optimizer.gBest_score    table[1, 3] += end - start        x_max = 100*np.ones(d)    x_min = -100*np.ones(d)    optimizer = PSO(fit_func=Step, num_dim=d, num_particle=p, max_iter=g, x_max=x_max, x_min=x_min)    start = time.time()    optimizer.opt()    end = time.time()    table[0, 4] += optimizer.gBest_score    table[1, 4] += end - start        x_max = 1.28*np.ones(d)    x_min = -1.28*np.ones(d)    optimizer = PSO(fit_func=Quadric_Noise, num_dim=d, num_particle=p, max_iter=g, x_max=x_max, x_min=x_min)    start = time.time()    optimizer.opt()    end = time.time()    table[0, 5] += optimizer.gBest_score    table[1, 5] += end - start        x_max = 500*np.ones(d)    x_min = -500*np.ones(d)    optimizer = PSO(fit_func=Schwefel, num_dim=d, num_particle=p, max_iter=g, x_max=x_max, x_min=x_min)    start = time.time()    optimizer.opt()    end = time.time()    table[0, 6] += optimizer.gBest_score    table[1, 6] += end - start        x_max = 5.12*np.ones(d)    x_min = -5.12*np.ones(d)    optimizer = PSO(fit_func=Rastrigin, num_dim=d, num_particle=p, max_iter=g, x_max=x_max, x_min=x_min)    start = time.time()    optimizer.opt()    end = time.time()    table[0, 7] += optimizer.gBest_score    table[1, 7] += end - start        x_max = 5.12*np.ones(d)    x_min = -5.12*np.ones(d)    optimizer = PSO(fit_func=Noncontinuous_Rastrigin, num_dim=d, num_particle=p, max_iter=g, x_max=x_max, x_min=x_min)    start = time.time()    optimizer.opt()    end = time.time()    table[0, 8] += optimizer.gBest_score    table[1, 8] += end - start        x_max = 32*np.ones(d)    x_min = -32*np.ones(d)    optimizer = PSO(fit_func=Ackley, num_dim=d, num_particle=p, max_iter=g, x_max=x_max, x_min=x_min)    start = time.time()    optimizer.opt()    end = time.time()    table[0, 9] += optimizer.gBest_score    table[1, 9] += end - start        x_max = 600*np.ones(d)    x_min = -600*np.ones(d)    optimizer = PSO(fit_func=Griewank, num_dim=d, num_particle=p, max_iter=g, x_max=x_max, x_min=x_min)    start = time.time()    optimizer.opt()    end = time.time()    table[0, 10] += optimizer.gBest_score    table[1, 10] += end - start    table = table / timestable = pd.DataFrame(table)table.columns=['Sphere', 'Schwefel_P222', 'Quadric', 'Rosenbrock', 'Step', 'Quadric_Noise', 'Schwefel',                 'Rastrigin', 'Noncontinuous_Rastrigin', 'Ackley', 'Griewank']table.index = ['mean score', 'mean time']

2.4.2 结果

同APSO那篇作者的测试条件，设D=30；P=20；迭代10,000次；个别运作30次，PSO结果如下

3. 小结

跟APSO提供的解相比，明显看得出有许多地方需要改进