本文旨在介绍粒子群优化算法(PSO)及其在地震学上的应用。

PSO原理简介

粒子群优化算法(PSO)是一个用以求解优化问题的随机演化算法。
科学家从鸟群的觅食等获得启发,将信息共享融入到最优化问题求解中,最终便发展出了一套简单,自适应而且灵活多样的优化算法。
粒子群最优化算法基于群体信息共享,从个体和群体搜索解空间的历史中自我学习,以此不断缩小搜索区域。
PSO基本流程可表述为如下步骤:

  • 以空间矢量代表单个粒子,其长度取决于求解问题的自由度。
  • 设定群体中单个粒子数目,各个粒子的初始随机位置\(b^{0}_{i}\)和速度\(v_{i}^{0}\)。对于粒子\(i\)计算其目标函数\(S_i\)
  • 时间向前演化,各个粒子的位置和速度变化与单个粒子的搜索历史和其周边粒子的搜索历史的目标函数相关。时间演变的第\(k+1\)步,粒子\(i\)按照如下公式更新自己的位置和速度:$$\mathbf{V}^{k+1}_i = \mathbf{V}^{k}_i + a_l r_1 (\mathbf{l}^{k}_i - \mathbf{x}^{k}_i) +a_{g} r_2 (\mathbf {g}^{k} - \mathbf{x}^{k}_i)$$
$$\mathbf{x}^{k+1}_i = \mathbf{x}^{k}_i + \mathbf{v}^{k+1}_i $$

其中,$\mathbf{l}^{k}_{i}$表示第$i$个粒子在之前$k$次迭代中目标函数的最佳拟合,$\mathbf{g}^{k}$表示所有粒子在之前$k$次迭代中对目标函数的最佳拟合,$a_l,a_g$表示局部及全局加速常量。$r_1,r_2$是分布在$(0,1)$区间内的随机数,用以对$a_l,a_g$进行加权。

为增强粒子探索解空间的能力以及削弱初始值对粒子探索的影响,科学家在式中加入惯性权重$\omega$。则该式重新表述如下:

$$\mathbf{V}^{k+1}_i = \omega \mathbf{V}^{k}_i + a_l r_1 ( \mathbf{l}^{k}_i - \mathbf{x}^{k}_i) + a_{g} r_2 ( \mathbf{g}^{k} -\mathbf{x}^{k}_i)$$

其中,$\omega$是实数,其在迭代过程中可保持不变亦可逐渐衰减。

为限制粒子速度,前人在式中引进限制因子$\chi$。限制因子的引入客观上控制了搜索半径。
前人研究表明:限制因子的引入能显著提升算法的性能。
$\chi$可通过认知参数$a_l$和社会参数$a_g$进行设定:

$$\mathbf{V}^{k+1}_i = \chi{\mathbf{V}^{k}_i + a_l r_1 ( \mathbf{l}^{k}_i - \mathbf{x}^{k}_i) + a_{g} r_2 ( \mathbf{g}^{k} -\mathbf{x}^{k}_i) }$$ $$\chi = \frac{2}{|2-\varphi-\sqrt{ \varphi^{2} -4 \varphi }|} \quad \text{其中} \quad \varphi = a_l+a_g,\varphi>4,\chi \in [0,1]$$

在PSO参数选择中,前人提供了一些比较理想参数搭配:$(\omega,a_l,a_g)=(0.8,1.8,2.0)$;
$(\omega,a_l,a_g)=(0.729,1.494,1.494)$;$(\omega,a_l,a_g)=(0.6,1.7,1.7)$;$(\chi,a_l,a_g)=(0.729,2.8,1.3)$;
$(\chi,a_l,a_g)=(0.729,2.05,2.05)$;
但是,在个人的实际地球物理反演问题中,参数仍旧需要随研究问题的变化而进行调整。

在由频散曲线起始的面波反演中,目标函数($S$)定义如下:

$$S=\lVert \mathbf{V_R}^{obs} - \mathbf{V_R}^{theo} \rVert / \sqrt{m}$$

其中,$\mathbf{V_R}^{obs}$是$m \times 1$阶的矢量,表示观测到的相速度频散曲线。
$\mathbf{V_R}^{theo}$是$m \times 1$阶的矢量,表示瑞利波理论频散。
$m$是频散曲线上所取的点数(一个频率对应于一个相速度)。$\lVert \rVert$表示矢量的欧几里得范数。

Tip:对于面波反演的主要介绍参照文章 Application of particle swarm optimization to interpret Rayleigh wave dispersion curves.

Pyswarm

Pyswarm 安装

自动安装或升级

  • 若用户按照ObsPy-安装 中所述,安装了Anaconda。则可直接在terminal中输入如下代码以安装:

    1
    pip install --upgrade Pyswarm
  • 若用户拥有setuptools ,则可尝试如下代码:

    1
    easy_install --upgrade pyswarm

网络下载及安装

用户可直接从pipy.python.org ,下载 pyswarm相关的压缩包,解压缩后进入其路径,运行setup.py即可。

  • 需安装权限的运行方式

    1
    python setup.py install
  • 不需额外权限的运行方式(需在用户的python库中进行)

    1
    python setup.py install --user
  • 或在用户自定义路径my_directory下安装

    1
    python setup.py install --install-lib my_directory
  • 当系统提示需要安装权限时(Unix)

    1
    sudo python setup.py install

源码编译

为获得最新的工作版本,用户可从GitHub 上直接clone代码并安装即可。

pycharm 范例

完整代码如下:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
from pyswarm import pso
def banana(x):
x1 = x[0]
x2 = x[1]
return x1**4 - 2*x2*x1**2 + x2**2 + x1**2 - 2*x1 + 5
def con(x):
x1 = x[0]
x2 = x[1]
return [-(x1 + 0.25)**2 + 0.75*x2]
lb = [-3, -1]
ub = [2, 6]
xopt, fopt = pso(banana, lb, ub, f_ieqcons=con)
# Optimum should be around x=[0.5, 0.76] with banana(x)=4.5 and con(x)=0

代码分段说明:

  • 载入粒子群优化模块

    1
    from pyswarm import pso
  • 定义将要获得最优最小值的待优化函数myfunction(x, *args, **kwargs).本函数需提供x,其类型应近似于1-D数组。本函数应返回一个标量值,该标量值便是应达到最小。本范例中,待优化函数即为banana

    1
    2
    3
    4
    def banana(x):
    x1 = x[0]
    x2 = x[1]
    return x1**4 - 2*x2*x1**2 + x2**2 + x1**2 - 2*x1 + 5

    数学形式为

    $$f_{banana}(\mathbf{x}) = x_{1}^4 - 2x_{2}x_{1}^2 + x_{2}^2+x_{1}^2-2x_{1}+5$$
  • 对于解是否需要满足某限制条件是可选参数,该范例中选择了con函数作为限制。

    1
    2
    3
    4
    def con(x):
    x1 = x[0]
    x2 = x[1]
    return [-(x1 + 0.25)**2 + 0.75*x2] #** constrain function

    该式表明,所有可接受的最优化解均应满足条件:

    $$ -(x_1 + 0.25 )^2 + 0.75x_2 \geq 0$$
  • 限定函数解的定义域的上下限
    1
    2
    lb = [-3, -1]
    ub = [2, 6]
    该式标定出了解的$x$的定义域:$$-3 \leq x[0] \leq 2 \\ -1 \leq x[1] \leq 6$$
  • 参数准备工作完成后即可调用pso函数:
    1
    xopt, fopt = pso(banana, lb, ub, f_ieqcons=con)

Pyswarm参数解析

所有调用参数如下:

1
2
3
pso(func, lb, ub, ieqcons=[], f_ieqcons=None, args=(), kwargs={},
swarmsize=100, omega=0.5, phip=0.5, phig=0.5, maxiter=100, minstep=1e-8,
minfunc=1e-8, debug=False)

  1. 必选参数

    1
    2
    3
    4
    5
    6
    func: function
    需要寻求其最小值的函数
    lb : array
    自变量定义域的下边界
    up : array
    自变量定义与的上边界
  2. 可选参数

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    ieqcons: list
    自变量需要满足的限制条件的列表,其长度为n。
    且需满足条件ieqcons[j](x,*args)>= 0.0 (j=1,...n)。
    参数默认为空列表([])。
    f_ieqcons: function
    自变量需要满足的限制条件,该函数将返回一个一维数组,其中
    每个参数均不小于零。参数默认为空(None)。
    args : tuple
    额外的需传递给对象或限制函数的参数,默认为空元组(())。
    kargs : dict
    额外的需传递给对象和限制函数的参数,默认为空字典({})。
    swarmsize : int
    粒子群中粒子个数,默认100个。
    omega : scalar
    粒子速度变化的权重参数即惯性权重,默认值为0.5
    phip : scalar
    粒子从该粒子自己认知历史的最佳位置迭代搜索的权重参数即认知参数,默认为0.5。
    phig : scalar
    粒子从粒子群的最佳位置迭代搜索的权重参数即社会参数,默认为0.5。
    maxiter: int
    粒子群迭代操作的最多次数,默认为100次。
    minstep: scalar
    迭代终止前,粒子群最佳位置前进的最小距离,默认1e-8。
    minfunc: scalar
    迭代终止前,粒子群最佳位置对应的函数减小的最小数值,默认1e-8。
    debug : boolean
    若为Ture,则将输出每一次迭代的结果,默认为False。

Note:在本例中,如下代码是等价的:

1
2
3
4
5
xopt, fopt = pso(banana, lb, ub, ieqcons=[con])
and
xopt, fopt = pso(banana, lb, ub, f_ieqcons=con)

参考文献

修改历史

  • 2017-01-21:初版
  • 2017-03-11:添加pyswarm工具