【PSO】基本PSO算法

基础粒子群算法原理

PSO是一种基于群体的随机化技术，通过初始化一组随机解，通过迭代搜索最优解，PSO算法通过模拟社会，将每个可能产生的解表述为群中的一个微粒，每个微粒具有独自的位置向量和速度向量，以及和目标函数有关的适应度，所有粒子在搜索空间中以一定速度飞行，通过追随当前搜索到的最优值来找到全局最优值.
PSO模拟社会根据如下三条规则进行：

• 飞离最近的个体
• 飞向目标
• 飞向群体的中心

PSO采用“群体”和“进化”的概念，根据个体的适应值进行操作，但是PSO中没有进化算子，粒子的飞行速度根据个体飞行经验和群体的飞行经验进行动态调整.

设在$S$维的目标搜索空间中，有$m$个粒子组成一个群体，其中第$i$个粒子表示为一个$S$维的向量$x_i=(x_{i1},x_{i2},\dots ,x_{is}),i=1,2,\dots ,m$，每个粒子的位置是一个潜在的解分量，目标函数值$f(x_i)$给出了解质量信息.
第$i$个粒子的速度矢量是$V=(V_{i1},V_{i2},\dots ,V_{is})$，记第$i$粒子目前搜索到的最优位置为$P_{is}=(P_{i1},P_{i2},\dots ,P_{is})$，整个粒子群目前找到的最优位置为$P_{gs}=(P_{1s},P_{2s},\dots ,P_{ms})$.
粒子$i$的当前最好位置由如下方程确定
$$p_i(t+1)=\{\begin{array}{l}{p}_i(t)\to f(x_i(t+1))\ge f(p_i(t))\\ x_i(t+1)\to f(x_i(t+1))<f(p_i(t))\end{array}\text{\hspace{1em}(1)}$$
Kennedy和Eberhart设置如下进化方程
$$\begin{array}{rl}& v_{is}(t+1)=v_{is}(t)+c_1{r}_{1s}(t)(p_{is}(t)-x_{is}(t))+c_2{r}_{2s}(t)(p_{gs}(t)-x_{is}(t))\\ & x_{is}(t+1)=x_{is}(t)+v_{is}(t+1)\end{array}$$
其中，$i\in [1,m],s\in [1,S]$；$c_1\ge 0$和$c_2\ge 0$是学习因子，$c_1$调节找到自身最好位置的步长，$c_2$调节找到全局最好位置的步长；$r_1$和$r_2$互为独立的服从$[0,1]$均匀分布的随机数；$v_{is}\in [-v_{\max },v_{\max }]$，如果限制搜索空间为$[-x_{\max },x_{\max }]$，可以设置$v_{\max }=k{x}_{\max },(0<k\le 1.0)$.
Y.shi和Eerhart改进了迭代方程
$$v_{is}(t+1)=w\cdot v_{is}(t)+c_1{r}_{1s}(p_{is}(t)-x_{is}(t))+c_2{r}_{2s}(t)(p_{gs}(t)-x_{gs}(t))$$
其中动量$w$为非负值，控制前一时刻的速度和当前速度的关系，如果$w$较大，前速度影响较大，模型的全局搜索能力较强，如果$w$较小，前速度影响较小，局部搜索能力较强，在处理局部最优问题时，需要调整$w$的值.
算法终止条件为达到最大迭代次数或者满足阈值要求.

算法流程

(1). 初始化规模为$m$的粒子群，设置位置和速度矢量. 对于任意的$i,s$，在$[-x_{\max },x_{\max }]$内服从均匀分布产生$x_{is}$；对于任意的$i,s$，在$[-v_{\max },v_{\max }]$内服从均匀分布产生$v_{is}$，$y_i=x_i,\forall i$.
(2). 计算每个粒子的适应值.
(3). 对每个粒子将其适应值和其经历过的最好位置$p_{is}$的适应值比较，如果当前适应值更好则将其作为当前最好位置.
(4). 对每个粒子将其适应值和全局经历过的最好位置$p_{gs}$的适应值比较，如果当前适应值更好则将其作为当前最好位置.
(5). 根据方程$(1)$更新粒子速度和位置
(6). 如果满足终止条件，输出解，否则转(2).

算法收敛性

算法全局收敛需要满足如下假设
(1). $f(D(z,\xi )\le f(z))$，并且如果$\xi \in S$, 则
$$f(D(z,\xi ))\le f(\xi )$$
随机算法全局收敛表示序列 { f ( z k ) } k = 1 ∞ \{f(z_k)\}_{k=1}^\infty {f(zk​)}k=1∞​需要收敛到$\psi$，其中
ψ = inf ⁡ { x ∣ v ( z ∈ S ∣ f ( S ) < x ) > 0 } \psi=\inf\{x\mid v(z\in S\mid f(S)<x)>0\} ψ=inf{x∣v(z∈S∣f(S)<x)>0}
由于Lebesgue测度$v(x\mid f(z)<x)>0$，可以避免病态情况.
(2). 对于$S$任意Borel子集$A$，如果测度$v(A)>0$，则有
$$\prod _{k=0}^{\infty }(1-{\mu }_k(A))=0$$
其中，${\mu }_k(A)$是由测度${\mu }_k$得到$A$的概率.
(3). 略（具体见参考资料page 158）
可以证明粒子群算法满足(1), 但是不满足(2)，所以算法不是全局收敛的.

案例：极值求解

PSO算法流程如图所示

使用fmincon和fminbnd以及Global Search进行求解，找到一维函数$f(x)$的全局最小值
$$f(x)=x\sin (x)+x\cos (x),x\in [0,10]$$
在Octave Online中实现代码如下
fmincon设置不同起始点寻找最优解，其中十字符号表示起点，圆圈符号表示在该起点下找到的最优值(最小值)的坐标.

pkg load optim
clear all
clc

%% objective function
f=@(x) x.*sin(x)+x./cos(x)
lb=0;
ub=10;
x0=[0 1 3 6 8 10];
fig=figure;
for i=1:6
    x(i)=fmincon(f, x0(i), [], [], [], [], lb, ub, [], optimset('Algorithm', 'sqp', 'Disp', 'none'))
    subplot(2, 3, i)
    ezplot(f, [lb, ub]);
    hold on
    plot(x0(i), f(x0(i)), 'k+')
    plot(x(i), f(x(i)), 'ro')
    hold off
    title(['Initial Point at ', num2str(x0(i))])
    if i==1 || i==4
        ylabel('x*sin(x)+x/cos(x)')
    end
end

使用fminbnd函数求解，不需要设置起始点的位置

pkg load optim
clear all
clc
close all

%% objective function
f=@(x) x.*sin(x)+x./cos(x)
lb=0;
ub=10;

fig=figure;
x2=fminbnd(f, lb, ub)
ezplot(f, [lb, ub])
hold on
plot(x2, f(x2), 'ro')
hold off
ylabel('x*sin(x)+x/cos(x)')
title({'Solution using fminbnd.', 'Required no initial point'})

Matlab中提供了全局最优搜索函数createOptimProblem，求解代码如下

problem=createOptimProblem('fmincon', 'objective', f, 'x0', x0(1), 'lb', lb, 'ub', ub, 'options', optimset('Algorithm', 'SQP', 'Disp', 'none'));
gs=GlobalSearch;
xgs=run(gs, problem)
fig=figure;
ezplot(f, [lb, ub]);
hold on
plot(xgs, f(xgs), 'ro')
hold off
ylabel('x sin(x)+x/cos(x)')
title('Solution using global Search.')

参考资料

Matlab优化算法案例分析与应用 清华大学出版社 余胜威