粒子群优化算法(Partical swarm optimization, PSO)是由Kennedy等^[1]提出的一种基于群体智能的随机优化方法, 为改善标准PSO算法在处理复杂多峰搜索问题时的全局寻优能力, 学者们对此进行深入地研究和分析, 吕强等设计粒子行为相应地改善了算法的寻优能力^[2]; Wang等提出了用拉丁矩阵设计粒子初始分布^[3]; Lu等提出了使用控制理论设计粒子群结构的新方法^[4].从PSO算法^[2-5]不难看出: PSO的目标函数作为评价函数参与调节目标函数变量的过程.从控制理论角度看, PSO算法相当于控制器, 而目标函数相当于控制对象, PSO算法解优化问题可以看成控制系统解优化问题.这为PSO算法开辟新的发展空间.一些学者已开始控制系统用于优化问题的研究, Ustundag用简单的模糊逻辑控制器解决单反馈控制系统的一维优化问题^[6], Jaewook等用两个阶段方法解决有少数局部极小值的全局优化问题^[7], Nader等提出了一类空间分布系统用于解决带线性不等式约束的大型多参数二次规划问题^[8]. Angelia和Ion用下降方法构造不同的神经网络解决简单的约束优化问题^[9-10].但这些方法仍无法应用于复杂的全局最优化问题^[11-13].我们知道:对用于优化的控制系统, 影响其输出的最重要的因素是控制器和控制对象(目标函数).本文结合PSO算法设计了新模糊控制器, 它可以使每个粒子具有不同的规则分布, 其模糊规则提高算法的全局搜索能力与局部搜索能力, 改善PSO的早熟收敛缺陷; 在不影响控制对象全局最优解的条件下, 提出了一种变换函数对其进行简化, 减少寻优过程陷入局部极小值的可能; 构造的新填充函数能使寻优过程很好地跳出局部极小点至比它更小的点处^[14-15].本文提出的基于变换函数与填充函数的模糊粒子群优化算法在保持粒子群算法简单, 收敛速度快等特点的同时, 有效地解决了陷入局部t最优的缺点, 测试函数验证了该算法的有效性.

1.   多回路模糊控制系统

本文考虑的优化问题为

这里$f(x)$为可导连续函数.

1.1   单回路控制系统

问题(1)的单回路控制系统如图 1所示, 其中$f(x)$为被控对象, $R$为常数输入.这里采用$T(a, b, x)$改变目标函数的表现形式而不改变目标函数的极值点位置, 使目标函数的最优值与局部极小值拉开差距, 减少计算过程陷入局部最优的可能, 使寻优问题变得快速简单.变换函数的引入必须保证目标函数变换前后的极值点位置不变.

图 1  单回路控制系统

Fig. 1  The single loop control system

下载: 全尺寸图片 幻灯片

本文中反馈变换为$T(a, b, x)=a(x-b)/{\sqrt{1+(a(x-b))^2}}~(a>0)$, 其为可导连续函数.目标函数的极值点为$\frac{\partial f(x)}{\partial x}=0$的解, 而变换后的函数为$T[f(x)]={a(f(x)-b)}/{\sqrt{1+(a(f(x)-b))^2}}, (a>0)$, 它的极值点是$\frac{\partial T[f(x)]}{\partial x}=\left(\frac{\partial T}{\partial f}\right)\times\frac{\partial f(x)}{\partial x}=0$的解.

由式(2)知$\frac{\partial T}{\partial f} >0$, 所以$\frac{\partial T[f(x)]}{\partial x} =\frac{\partial T}{\partial f}\times\frac{\partial f(x)}{\partial x}=0\Longleftrightarrow\frac{\partial f(x)}{\partial x}=0$, 也就是说$T[f(x)]$与$f(x)$有相同的极值点, $T$变换函数不改变原目标函数的极值点位置.

从图 2可以看出, $b=0$时, 变换函数$T(a, b, x)$在零点附近基本成正比, 而远离零点的区域变换函数的值趋于1.当变换函数作用于目标函数时, 其作用前和作用后的目标函数曲线在零附近形状变化不大, 而目标函数值远离零的目标函数曲线趋于1, 其形状变化大, 呈水平状.这样只要关注由变换得到的新目标函数在零附近的全局最小点即可, 从而减少计算过程陷入局部极小点的可能; 且$a$越大, 有效搜索范围越小; 反之, $a$越小, 有效搜索范围越大.以一维目标函数为例讲述变换函数参数$b$的作用:由图 2可知, 当$a=1$时, 如果变换函数的自变量大于10, 则变换函数值趋于1, 因此其自变量有效变化区间为[-10, 10].如图 3所示, 当$b=0$时, $ T[f(x)]=f(x)/\sqrt{1+f(x)^2}$, $f(x)=10$的直线$l_i$与$f(x)$有两个交点, 此时的有效搜索范围应为$\Delta_i$, 而当$b=1$时, $ T[f(x)]=(f(x)-1)/\sqrt{1+(f(x)-1)^2}$, $f(x)-1=10$即$f(x)=11$的直线$l_j$与$f(x)$有两个交点, 此时的有效搜索范围应为$\Delta_j$, 由图 3可知: $\Delta_i < \Delta_j$, 即对于变换函数$ T(a, b, x)=a(x-b)/\sqrt{1+(a(x-b))^2}~ (a>0)$, 当$a$固定时, $b$越大, 有效搜索范围$\Delta_j$越大; 反之, $b$越小, 有效搜索范围$\Delta_i$越小.

图 2  变换函数图

Fig. 2  Transformation function curve

下载: 全尺寸图片 幻灯片

图 3  目标函数平面示意图

Fig. 3  Objective function diagram

下载: 全尺寸图片 幻灯片

当取$a=1$, $b=0$, 则$ T(a, b, x)=T(1, 0, x)=x/\sqrt{1+x^2}$. $T$变换前后的对比如图 4和5所示, 其中图 4(a)和5(a)为表 1中的$F4$函数和$F11$函数, 可以看出它们有很多局部极值点, 计算过程很容易陷入局部极值点; 而图 4(b)和5(b)为对应$T$变换后$T[f(x)]$图, 可以看出原函数中很多局部极值被削减, 而变换后最优值更好地突显出来, 使计算过程更快速高效的找到最优值.所以变换函数$T(a, b, x)=a(x-b)/\sqrt{1+(a(x-b))^2}$, 用于反馈环节, 减少算法陷入局部极值点的可能性.

图 4  $F4$变换前后的曲线图即$F4$和$T(1,0,F4)$

Fig. 4  The curves of $F4$ before and after transformation: $F4$ and $T(1,0,F4)$

下载: 全尺寸图片 幻灯片

图 5  $F11$变换前后的曲线图即$F11$和$T(1,0,F11)$

Fig. 5  The stereogram of $F11$ before and after transformation: $F11$ and $T(1,0,F11)$

下载: 全尺寸图片 幻灯片

表 1  测试函数

Table 1  Test functions

测试函数	维数	可行域	最优值/最优点
$f_1=\sum\limits_{i=1}^n(x_i^2-10\cos(2\pi x_i)+10)$	30	$[-5.12,5.12]^D$	0.0/0, 0, $\cdots$, 0
$f_2=-20\exp(-0.2\sqrt{\sum\limits_{i=1}^n\frac{x_i^2}{n}})-\exp(\sum\limits_{i=1}^n\frac{\cos(2\pi x_i)}{n})+20+\exp(1)$	30	$[-32,32]^D$	0.0/0, 0, $\cdots$, 0
$f_3=\sum\limits_{i=1}^n\frac{x_i^2}{4\,000}-\prod\limits_{i=1}^n\cos(\frac{x_i}{\sqrt{i}})+1$	30	$[-600,600]^D$	0.0/0, 0, $\cdots$, 0
$f_4=\frac{\pi}{n}(10\sin^2(\pi y_1)+(y_n-1)^2+\sum\limits_{i=1}^{n-1}(y_i-1)^2(1+10\sin^2(\pi y_{i+1}))) y_i=1+\frac{(1+x_i)}{4}$	30	$[-50,50]^D$	$ 0.0/-1,\ -1, \cdots,\ -1 $
$ f_5=\frac{1}{10}(\sin^2(3\pi x_1)+(x_n-1)^2(1+\sin^2(2\pi x_n))+ \quad \sum\limits_{i=1}^{n-1}(x_n-1)^2(1+\sin^2(3\pi x_{i+1})))$	30	$[-50,50]^D$	0.0/1, 1, $\cdots$, 1
$f_6=\sum\limits_{i=1}^n{\frac{(x_i^4-16x_i^2+5x_i)}{n}}$	30	$[-5,5]^D$	$-78.3323/-2.90353,-2.90353, \cdots,-2.90353$
$f_7=\sum\limits_{i=1}^{n-1}\left({100(x_i^2-x_{i+1})^2+(x_i-1)^2}\right)$	30	$[-5,10]^D$	0.0/1, 1, $\cdots$, 1
$f_8=\sum\limits_{i=1}^nx_i^2$	30	$[-100,100]^D$	0.0/0, 0, $\cdots$, 0
$f_9=\sum\limits_{i=1}^nx_i^4$	30	$[-1.28,1.28]^D$	0.0/0, 0, $\cdots$, 0
$f_{10}=\sum\limits_{i=1}^n\left[{\sum\limits_{j=1}^i x_j}\right]^2$	30	$[-100,100]^D$	0.0/0, 0, $\cdots$, 0
$f_{11}=\sum\limits_{i=1}^n\left[\sum\limits_{k=0}^{k\max}\left(a^k\cos(2\pi b^k(x_i+0.5))\right)\right]- \quad n\sum\limits_{k=0}^{k\max}\left(a^k\cos(2\pi b^k0.5)\right) a=0.5,b=3,k\max=20$	30	$[-0.5,0.5]^D$	0.0/0, 0, $\cdots$, 0
$f_{12}=418.9829\times n-\sum\limits_{i=1}^nx_i\sin(\sqrt{|x_i|})$	30	$[-500,500]^D$	0.000381827/420.97, 420.97, $\cdots$, 420.97
$f_{13}=\sum\limits_{i=1}^ni\times x_i^4+{\rm Random}(0,1)$	30	$[-1.28,1.28]^D$	0.0/0, 0, $\cdots$, 0
$f_{14}=\frac{\pi}{n}\Big(10\sin^2(\pi y_1)+(y_n-1)^2+\sum\limits_{i=1}^{n-1}(y_i-1)^2(1+ \quad 10\sin^2(\pi y_{i+1}))\Big) +\sum\limits_{i=1}^nu_i(x_i,10,100,4), \\ y_i=1+\frac{1+x_i}{4} u_i(x_i,a,k,m)=\left\{\!\!\!{array}{ll} k(x_i-a)^m, & x_i>a \\ 0, &-a \leq x_i \leq a \\ k(-x_i-a)^m, & x_i <a {array}\right. $	30	$[-50,50]^D$	$0.0/-1,-1,\cdots,-1 $
$f_{15}=\frac{1}{10}(\sin^2(3\pi x_1)+(x_n-1)^2(1+\sin^2(2\pi x_n))+\quad \sum\limits_{i=1}^{n-1}(x_n-1)^2(1+\sin^2(3\pi x_{i+1})))+\sum\limits_{i=1}^nu_i(x_i,5,100,4)$	30	$[-50,50]^D$	0.0/1, 1, $\cdots$, 1

下载: 导出CSV 

| 显示表格

1.2   多回路模糊控制系统

基于上述单回路控制系统, 我们可得多回路模糊控制系统如图 6所示, 一个子系统相当于一个粒子, FNN1, $\cdots$, FNN$m$为模糊控制器区别于$m$个不同的隶属度函数, 同时体现系统的全局搜索能力与局部搜索能力, $f(x)$为目标函数, ${\rm Best}(f(x))$为子系统的最优值(局部最优值), 各个子系统共享所有子系统的最优值, 根据共享信息迅速调整该子系统的寻优方向, ${\rm min}(f(x))$为整个系统的最优值(全局最优值).以$T$变换后的函数值作为反馈.

图 6  多回路模糊控制系统框图

Fig. 6  A multi-loop distributed fuzzy control system

下载: 全尺寸图片 幻灯片

2.   模糊粒子群算法

2.1   粒子群算法

在标准PSO算法中, 假设在一个$n$维的问题空间中, 包含$m$个粒子, 每个粒子作为搜索空间中待优化问题的一个可行解, 通过粒子之间的协作与竞争来寻找问题的最优解.在第$t$次迭代中, 第$\sigma$个粒子的当前位置表示为$x_\sigma(t)=(x_{\sigma1}(t), x_{\sigma2}(t), \cdots, x_{\sigma n}(t))$, 当前速度表示为$v_\sigma(t)=(v_{\sigma1}(t), v_{\sigma2}(t), \cdots, v_{\sigma n}(t))$, $\sigma=1, 2, \cdots, m$.在每次迭代中, 粒子个体搜索到的最好位置为$p_\sigma(t)=(p_{\sigma1}(t), p_{\sigma2}(t), \cdots, p_{\sigma n}(t))$, 称为个体最优; 群体中所有粒子搜索到的最好位置为$p_g(t)=(p_{g1}(t), p_{g2}(t), \cdots, p_{gn}(t))$, 称为全局最优.优化问题的过程可看作粒子不断更新的过程, 每个粒子根据式(3)以当前速度、个体最优和全局最优来调整自己的飞行速度和方向, 进化方程如下:

其中, $r1$、$r2$为[0, 1]上的随机数; $c1$、$c2$为学习因子是一常数, 表示粒子受自身及全局的影响程度, 调节向个体最优和全局最优方向飞行的最大移动步长.式(3)中的第一部分称为动量部分, 表示粒子对当前自身运动状态的信任程度, $\omega$称为惯性权重, 一般取[0.4, 0.9]之间的数^[1-5], 使其依据自身速度进行惯性运动; 第二部分称为个体认知部分, 引导粒子飞向自身的最好位置; 第三部分称为社会认知部分, 表示粒子间的信息共享与相互协作, 它引导粒子飞向群体中的最好位置.这三个部分之间的相互平衡和制约决定了算法的主要搜索性能.

2.2   模糊粒子群算法

因为传统PSO在全局极值$p_g$或个体极值$p_\sigma$得到局部最优解时, 粒子群不会在解空间中再次进行搜索, 而且其他粒子将迅速向局部最优解靠拢, 所以使算法容易出现过早收敛导致不能得到最优解.又因为传统PSO的全局搜索模式相对单一, 仅仅使用全局极值点的信息, 而没有加入其他的参考信息, 这样不利于种群多样性且阻碍算法扩大搜索的范围.

我们知道如果粒子数足够多, 并且每个粒子的初始位置分布完全分散, 则可以通过式(3)获得目标函数的全局最优解, 但它需要大量的CPU时间.为了克服这个问题, 我们必须减少粒子的个数, 使有限粒子的初始位置尽可能分离.因为模糊决策通常具有分布特征, 所以它可以使每个粒子具有规则的分布.基于此本文利用模糊控制器中所使用的模糊规则改进粒子群算法, 提高算法的全局搜索能力与局部搜索能力, 改善早熟收敛缺陷.根据Takagi--Sugeno模糊模型, 本文模糊控制系统的模糊规则采用如下:

$R_\sigma^l$表示第$\sigma$个模糊子系统的第$l$个模糊规则, $N_\sigma$表示第$\sigma$个模糊子系统的模糊规则数, $p$表示第$\sigma$个模糊子系统的模糊变量数, $x_{\sigma {\rm best}}$是第$\sigma$个模糊子系统的最好位置向量$f(x_{\sigma {\rm best}})=\mathop {\min }_{t\geq t_1}f(x_\sigma(t))$, 即是PSO中的$p_\sigma$, $x_{\rm allbest}$是所有模糊子系统的最好位置向量$f(x_{\rm allbest})=\mathop {\min}_{\sigma}\mathop{\min }_{t\geq t_1}f(x_\sigma(t))$, 即是PSO中的$p_g$, 根据PSO算法^[3-4], $x_\sigma(t)$是第$\sigma$个模糊子系统粒子的当前位置, $\omega_0$, $\omega_1$, $\omega_2$, $\omega_3$与$\omega_4$是第$\sigma$个模糊子系统的调整参数, 其中$\omega_0$, $\omega_1$, $\omega_2$分别代表第$\sigma$个模糊子系统的变量$x_\sigma$的零阶项系数, 一阶项系数, 变化差的系数, 这相当于PI控制器系数; 而$\omega_3$与$\omega_4$是相当于PSO系数. $r_1$和$r_2$是[0, 1]范围内的随机数, 也就是说, 第$\sigma$个模糊子系统相当于PI控制和PSO的结合.采用加权平均法去模糊化, 第$\sigma$个模糊子系统的全局模型可描述为

其中, $h_{\sigma l}(\xi)=\frac{\prod_{i=1}^p\mu_{M_{\sigma p}^l}}{\sum_{l=1}^{N_\sigma}\prod_{i=1}^p\mu_{M_{\sigma p}^l}}$, $\sum_{l=1}^{N_\sigma}h_{\sigma l}(\xi)=1$, $\mu_{M_{\sigma p}^l}$为第$P$个模糊变量的隶属度函数, $\sigma=1, 2, \cdots, m$. $\mu_{M_{\sigma p}^l}=R_\sigma-T(a_{\sigma l}, b_{\sigma l}, f)$, $R_\sigma$为第$\sigma$个子系统的输入.

在本文中$h_{\sigma l}(\xi)=\frac{\prod_{i=1}^p\mu_{M_{\sigma p}^l}}{\sum_{l=1}^{N_\sigma}\prod_{i=1}^p\mu_{M_{\sigma p}^l}}=\frac{\mu_{M_{\sigma1}^l}}{\sum_{l=1}^{N_\sigma}\mu_{M_{\sigma1}^l}} (p=1)$.结合图 7阐述隶属度函数的设计思想:以常用的三角形隶属度函数为例, 图 7(a)中的隶属度函数宽度区间记为$\sigma_i$, 图 7(b)中的隶属度函数宽度区间记为$\sigma_j$, 若$\sigma_i < \sigma_j$, 相同数量的隶属度函数的搜索范围$\Delta_i < \Delta_j$.也就是说, 模糊控制系统的隶属度函数宽度区间$\sigma_i$越小, 该系统的局部搜索能力越强, 精度越高; 相反, 隶属度函数宽度区间$\sigma_j$越大, 该系统的全局搜索能力越强, 搜索范围越大.根据这种想法, 作者为$m$个模糊控制子系统设置$m$个隶属度函数, 如图 8所示.记第$i$个模糊控制子系统的隶属度函数宽度区间为$\sigma_i$, 第$j$个模糊控制子系统的隶属度函数宽度区间为$\sigma_j$, 令$\sigma_i < \sigma_j~(i < j, i, j=1, 2, \cdots, m)$.这样每个模糊控制子系统的隶属度函数宽度区间都不一样, 隶属度函数宽度区间窄的控制系统局部搜索能力强, 精度高; 隶属度函数宽度区间宽的控制系统全局搜索能力强, 搜索范围大, 很好地弥补了传统粒子群算法的缺陷.另外, 与粒子群算法思想类似, 所有模糊子系统之间信息共享, 第$\sigma$个模糊子系统的当前变量$x_\sigma$以上次迭代所得的局部最优值$x_{\sigma {\rm best}}$和全局最优值$x_{\rm allbest}$作为导向.

图 7  不同宽度的隶属度函数

Fig. 7  The membership functions with different width

下载: 全尺寸图片 幻灯片

图 8  多回路控制系统关系图

Fig. 8  The relationship with all subsystems

下载: 全尺寸图片 幻灯片

由此可知对于$m$个模糊控制子系统, 只需要设置变换函数的参数$a$的值就可以改变隶属度函数的宽度区间, 即$a_{\sigma_{il}}>a_{\sigma_{jl}}$, 则$\sigma_i < \sigma_j$.此时设置$m$个模糊控制子系统就是设置$a$为$m$行矩阵, $a_i$的值越大, 第$i$个子系统隶属度函数的宽度区间越窄, 该模糊控制子系统的局部搜索能力越强, 精度越高; 相反, $a_j$的值越小, 第$j$个子系统隶属度函数的宽度区间越宽, 该模糊子系统的全局搜索能力越强, 搜索范围越大.当离散步长取一个非常小的正数时, 式(5)离散化描述如下:

这里$\eta_4$是离散步长, 一个非常小的正数.式(6)的参数$\omega_0$, $\omega_1$, $\omega_2$, $\omega_3$和$\omega_4$修正如下:

其中, $\sigma=1, 2, \cdots, m$; $l=1, 2, \cdots, N_\sigma$; $\eta_1$, $\eta_2$, 是很小的正数; $m$是子系统个数.

收敛性分析:构造一个Lyapunov函数$g_\sigma(t)=L+f(x_\sigma(t))+0.5\eta_2\left\| {x_\sigma(t)-x_{\rm allbest}} \right\|^2\geq0$, $f(x_\sigma(t))+L\geq0, L=0$.

把式(5)和式(7)代入上式得:

很容易看出$\frac{{\rm d}[g_\sigma(t)]}{{\rm d}t}\leq0$.当$\frac{{\rm d}[g_\sigma(t)]}{{\rm d}t}=0$时, $\frac{\partial f(x_\sigma(t))}{\partial x_\sigma(t)}+\eta_2(x_\sigma(t)-x_{\rm allbest}(t))=0$或$\frac{{\rm d}[x_\sigma(t)]}{{\rm d}t}=0$.若$\frac{\partial f(x_\sigma(t))}{\partial x_\sigma(t)}+\eta_2(x_\sigma(t)-x_{\rm allbest}(t))=0$, 由式(7)得$\omega_0=\omega_1=\omega_2=\omega_3=\omega_4=0$.代入式(5)得$\frac{{\rm d}[x_\sigma(t)]}{{\rm d}t}=0$.综上所述, $\frac{{\rm d}[x_\sigma(t)]}{{\rm d}t}=0$时, 即$x_\sigma(t)$达到一个常值使得$\frac{{\rm d}[x_\sigma(t)]}{{\rm d}t}=0$.由Lyapunov稳定性定理和La Salle's定理知, 变量$x_\sigma(t)$根据式(5)能够收敛到一个常值.

当然变换函数的参数$a$和$b$在迭代中, 可由下式修正.

其中, $\eta_3$是很小的正数.

3.   填充函数法

上述算法改善了粒子群的全局搜索能力, 但也有陷入局部极值点的可能性.为此, 构造一个填充函数使其跳出局部极值点至更小的点, 从而继续迭代搜索直到全局最优值.

填充函数法的主要思想是:如果已经找到了一个局部极小点$x_1$, 但它不是全局最小点, 我们可以在$x_1$处构造一个填充函数使迭代点列离开$x_1$所在的谷域, 找到更好的点$x_2$~(即$x_2$处的目标函数值比$x_1$处的目标函数值更小), 然后以$x_2$为初始点极小化原问题找到更优的极小点.

如图 9所示, 根据目标函数$F(x)$的一个已知极小点$x_1$构造一个函数$P(x)$, 称为填充函数, 它在$x_1$处取得极大值, 而在任何比$x_1$的盆域B1高的盆域中没有极小点或鞍点, 但在某个比$x_1$的盆域低的盆域中有一个极小点或鞍点$x^\ast$, 然后以$x^\ast$作为初始点去极小化$F(x)$, 并找到$F(x)$的一个新的极小点$x_2$, 使得$F(x_2)\leq F(x_1)$, 用$x_2$代替$x_1$, 重复上述过程, 直到找到$F(x)$的全局极小点^[14-15].

图 9  填充函数示意图

Fig. 9  Graph of the filled function

下载: 全尺寸图片 幻灯片

常用的填充函数有以下几种:

这里$r$, $\rho$以及$a$是计算过程中需要调整的参数, 第一个填充函数为双参数, 调节起来比较麻烦、费时, 最后一个为单参数, 可是只有当$a$很大时, 才能成为$x1$处的填充函数, 由于参数$a$在指数上, 当$a$很大时, 计算性能不理想.作者构造了一种新的填充函数:

该填充函数是无参的, 节约冗长的计算步骤及调整参数的时间, 提高算法的效率.

以下说明$P(x, x_1)$满足填充函数的条件:

$1)$如果$x_1$是$f(x)$的一个已知的局部最小点, 则$x_1$是$P(x, x_1)$的一个严格的局部最大点.

证明.设$B_1$是$f(x)$在局部最小点$x_1$的一个$S$型盆域, 即对任意$x\in B_1$有$f(x)\geq f(x_1)$, 则$P(x, x_1)=-\frac{\left\| {x-x_1} \right\|}{1+\left\| {x-x_1} \right\|} < 0=P(x_1, x_1)$, 这表明$x_1$是$P(x, x_1)$的一个严格的局部最大点.

2) 如果$x_1$是$f(x)$的一个已知的局部最小点, 且$x\in S_1=\left\{x|f(x)\geq f(x_1), x\neq x_1\right\}$, 则$P(x, x_1)$在$S_1$域内无固定点或者最小点.

证明.因为$x\in S_1, P(x, x_1)=-\frac{\left\| {x-x_1} \right\|}{1+\left\| {x-x_1} \right\|}$, 则当$f(x)>f(x_1)$和$x\neq x_1$, 有$\nabla p(x, x_1)=-\frac{(x-x_1)}{(1+\left\| {x-x_1} \right\|)^2(\left\| {x-x_1} \right\|)}\neq0$.如果$d(x)=x-x_1$和$x\neq x_1$, 则$d^{\rm T}\nabla p(x)=-\frac{\left\| {x-x_1} \right\|}{(1+\left\| {x-x_1} \right\|)^2} < 0$.因此$d(x)$是对于任意$x\in S_1$, $P(x, x_1)$的一个速降方向.当$f(y)=f(x_1)$, 有$P(y, x_1)=-\frac{\left\| {y-x_1} \right\|}{1+\left\| {y-x_1} \right\|}$, 则对于满足$\left\| {x-x_1} \right\|>\left\| {y-x_1} \right\|$和$\left\| {x_2-x_1} \right\| < \left\| {y-x_1} \right\|$的任意$x, x_2\in S_1$, 有$\nabla p(x, x_1)=-\frac{\left\| {x-x_1} \right\|}{1+\left\| {x-x_1} \right\|} < -\frac{\left\| {y-x_1} \right\|}{1+\left\| {y-x_1} \right\|}=P(y, x_1)$, $\nabla p(x_2, x_1)=-\frac{\left\| {x_2-x_1} \right\|}{1+\left\| {x_2-x_1} \right\|}>-\frac{\left\| {y-x_1} \right\|}{1+\left\| {y-x_1} \right\|}=P(y, x_1)$, 因此, $p(x, x_1)$在$S_1$域内无稳定点或者最小点.

3) 如果$x_1$是$f(x)$的一个已知的局部最小点, 对于任意$x_2\in \Omega$, $x_2\neq x_1$, 如果$f(x_2)\geq f(x_1)$, 则$x_2$是$P(x, x_1)$的一个连续点, 或者$x_2$是$P(x, x_1)$的一个下半连续点.

证明.因为$x_1$是$f(x)$的$B_1$盆域内的一个局部最小点, 所以对于任意$x\in B_1$, 有$f(x_2)\geq f(x_1)$, 则$\lim_{x\rightarrow x_1}P(x, x_1)=\lim_{x\rightarrow x_1}-\frac{\left\| {x-x_1} \right\|}{1+\left\| {x-x_1} \right\|}=0=P(x_1, x_1)$.因此, $P(x, x_1)$在$x_1$点是连续的.如果$f(x_2)\neq f(x_1)$且$x_2\neq x_1$, 即$f(x_2) < f(x_1)$或$f(x_2)>f(x_1)$, 则$p(x_2, x_1)=\frac{\left\| {x_2-x_1} \right\|}{1+\left\| {x_2-x_1} \right\|}$或$p(x_2, x_1)=-\frac{\left\| {x_2-x_1} \right\|}{1+\left\| {x_2-x_1} \right\|}$, 这样存在一个邻域$N(x_2, \delta), \delta>0$使对于任意$x\in N(x_2, \delta)$, $\lim_{x\rightarrow x_1}P(x, x_1)=P(x_2, x_1)$, 即$P(x, x_1)$在点$x_2$是连续的.如果$f(x_2)=f(x_1)$且$x_2\neq x_1$, 则$p(x_2, x_1)=-\frac{\left\| {x_2-x_1} \right\|}{1+\left\| {x_2-x_1} \right\|}$, $\liminf_{x\rightarrow x_2}P(x, x_1)=-\frac{\left\| {x_2-x_1} \right\|}{1+\left\| {x_2-x_1} \right\|}$.因此$x_2$是$P(x, x_1)$的一个下半连续点.

4) 如果$x_1$是$f(x)$的一个已知的局部最小点, $x_2$是满足$f(x_2) < f(x_1)$的另一个与$x_1$邻近的局部最小点, 则存在一点$x^\prime\in S_2=\left\{x, f(x) < f(x_1), x\in \Omega\right\}$使得$P(x, x_1)$最小, 且它在$x_1$和$x^{\prime\prime}$的连接线上($x^{\prime\prime}$属于$x_2$某一个邻域).

证明.已知局部最小点$x_2$满足$f(x_2) < f(x_1)$, 则$p(x_2, x_1)=-{\rm sgn}(f(x_2)-f(x_1))\frac{\left\| {x_2-x_1} \right\|}{1+\left\| {x_2-x_1} \right\|}=\frac{\left\| {x_2-x_1} \right\|}{1+\left\| {x_2-x_1} \right\|}>0$; 对于局部极小点$x_1$, 存在$S$型盆域$B_1$, 存在$x_3\in B_1$满足$x_3 < x_1$, $f(x_3)>f(x_1)$, 则有$p(x_3, x_1)=-{\rm sgn}(f(x_3)-f(x_1))\frac{\left\| {x_3-x_1} \right\|}{1+\left\| {x_3-x_1} \right\|}=-\frac{\left\| {x_3-x_1} \right\|}{1+\left\| {x_3-x_1} \right\|} < 0$; 又$P(x_1, x_1)=0$.由式(3)证明知, 则存在一点$x^{\prime}\in S_2=\left\{x, f(x) < f(x_1), x\in \Omega\right\}$使得$P(x, x_1)$最小, 且它在$x_1$和$x^{\prime\prime}$的连接线上($x^{\prime\prime}$属于$x_2$某一个邻域).

5) 如果$x_1$是$f(x)$的一个已知的局部最小点, $d$是满足$d^{\rm T}(x-x_1)>0$的一个方向, 如果$f(x)>f(x_1)$, 则$d$是$P(x, x_1)$在点$x_1$的一个速降方向.

证明.因为$f(x)>f(x_1)$, $d^{\rm T}(x-x_1)>0$, 则$d^{\rm T}\nabla p(x, x_1)=-\frac{d^{\rm T}(x-x_1)}{(1+\left\| {x-x_1} \right\|)^2(\left\| {x-x_1} \right\|)} < 0$因此$d$是$P(x, x_1)$在点$x_1$的一个速降方向.

由上述式(1) $\sim$ 式(5)的论述证明函数$P(x, x_1)$满足构造函数的条件, 可以作为构造函数使用.令$d=-\eta\nabla P(x)(1+\left\| {x-x_1} \right\|)^2\left\| {x-x_1} \right\|$, 其中$\eta>0$, $x_1$是$f(x)$的一个已知的局部极小值点, $x\neq x_1$.

因此$d=-\eta\nabla P(x)(1+\left\| {x-x_1} \right\|)^2\left\| {x-x_1} \right\|$是$P(x, x_1)$在$x_1$点的一个速降方向, 在$d$方向上具有扩散性, 能够扩散到搜索范围内的所有点, 所以填充函数的迭代公式为

其中, $x_i^{(k+1)}$是第$k+1$次迭代位置, $x_i^{(k)}$是第$k$次的迭代位置, $x_1$是局部极小值点.当寻优过程陷入局部最优时$x_i^{(k)}=x_1$, 不能直接采用式(11)迭代, 而是令$x_i^{(k)}$加一个随机数, 使$x_i^{(k)}\neq x_1$再采用式(11)的形式迭代, 具体迭代公式见算法步骤.

4.   FPSO-TF算法流程

FPSO-TF (Fuzzy partical swarm optimization based on filled function and transformation function)算法流程如图 10所示, 具体步骤如下:

图 10  FPSO-TF算法流程图

Fig. 10  Flowchart of FPSO-TF algorithm

下载: 全尺寸图片 幻灯片

步骤1.给出变量$X_{n\times m}$, 系统输入$R$, 模糊规则数$N_\sigma$, 变换函数的参数$a$, 学习系数$\eta_1, \eta_2, \eta_3, \eta_4$的初值, 算法最大迭代次数$MN$, 进入填充函数的最大次数$TT$, 在填充函数中变量迭代的最大次FMN, 以及填充函数迭代参数$d_1~(0 < d_1 < 1E-5), d_2~(0 < d_2 < 1), d_3~(1 < d_3), d_4~(0 < d_4 < 1), d_5~(1 < d_5)$; 令$L=0, j=0, k=1$;

步骤2.令$h(x)=f(x)$, 由模糊粒子群迭代式(6)更新$X_{n\times m}$, 使$j=j+1$, 并计算$f(x)=[f_1, f_2, \cdots, f_m](f_\sigma(x)=T_\sigma(f(x)), \sigma=1, 2, \cdots, m)$; 如果$f_\sigma < f_{\sigma {\rm best}}(\sigma=1, 2, \cdots, m)$, 则$x_{\sigma {\rm best}}=x_\sigma, f_{\sigma {\rm best}}=f_\sigma$; 如果$f_{\sigma {\rm best}} < f_{\rm allbest}(\sigma=1, 2, \cdots, m)$, 则$x_{\rm allbest}=x_{\sigma {\rm best}}, f_{\rm allbest}=f_{\sigma {\rm best}}$;

步骤3.如果$\left|h_\sigma(x)-f_\sigma(x)\right| < d_1, k\leq FMN$且$f_\sigma(x)\geq f_{\rm allbest}$, 则进入填充函数循环迭代: $k=k+1$, $x_\sigma(i, k)=x_\sigma(i, k)+d_2({\rm rand}(1)-0.5)\exp (d_3(k-1)/FMN), \sigma=1, 2, \cdots m; i=1, 2, \cdots, n$, $x_\sigma(i, k)=x_\sigma(i, k)+d_4({\rm rand}(1)+1)\exp (d_5(k-1)/FMN)(x_\sigma(i, k)-x_{\rm allbest})$, 并计算$f(x)=[f_1, f_2, \cdots, f_m](f_\sigma(x)=T_\sigma(f(x)), \sigma=1, 2, \cdots, m)$; 如果$f_\sigma < f_{\sigma {\rm best}}$ $(\sigma=1, 2, \cdots, m)$, 则$x_{\sigma {\rm best}}=x_\sigma, f_{\sigma {\rm best}}=f_\sigma$; 如果$f_{\sigma {\rm best}} < f_{\rm allbest}(\sigma=1, 2, \cdots, m)$, 则$x_{\rm allbest}=x_{\sigma {\rm best}}, f_{\rm allbest}=f_{\sigma{\rm best}}$; 转至步骤3;否则, 若$k>FMN$, 则$L=L+1$, 转至步骤4;若$\left|h_\sigma(x)-f_\sigma(x)\right|>d_1$或$f_\sigma < f_{\rm allbest}$, 转至步骤4;

步骤4.根据式(7)更新参数$\omega_0, \omega_1, \omega_2, \omega_3, \omega_4$; 如果$j>MN$或$L>TT$, 则输出$f_{\rm allbest}(x_{\rm allbest})$, 算法结束; 否则, 转至步骤2.

5.   仿真

为了测试FPSO-TF算法的寻优能力, 本文主要从两个方面来考核算法的性能: 1)更高维函数测试; 2)均值与方差.本文选取了近年来文献给出的各种改进优化算法作为对比算法, 用来作为评估本论文给出的算法性能的参照.这里表 1给出了用于测试的15个基准函数, 表 2给出了15个测试函数初值, 其中变量$x$的初值是在可行域内随机产生的, 而参数$\omega_0, \omega_1, \omega_2, \omega_3, \omega_4$由式(7)得到. $a$矩阵为$m\times N_\sigma$的矩阵, $N_\sigma$表示第$\sigma$个模糊子系统的模糊规则数, 初值如下:

表 2  参数初值

Table 2  The initial values of the parameters

测试函数	$d_1/d_2/d_3/d_4/d_5$	$\eta_1/\eta_2/\eta_3/\eta_4$	$a_1/R/N_\sigma$	$MN/TT/FMN$
$F1 $	1E-6/1E-18/30/1/12	2E-2/1E-4/1E-4/0.1	1E-3/-1/6	1 000/2 000/100
$F2 $	1E-4/1E-26/40/0.1/30	5E-3/1E-4/1E-5/0.1	1E-4/-1/5	5 000/8 000/100
$F3 $	1E-6/1E-15/25/0.3/18	1E-4/1E-3/1E-3/1E-2	1E-3/-1/3	15 000/9 000/60
$F4$	1E-4/1E-18/35/0.2/15	1.6/0.1/1E-5/0.1	1E-4/-1/5	12 000/20 000/100
$F5 $	1E-4/1E-18/35/0.05/18	1.4/1E-2/1E-5/0.1	1E-3/0/3	8 000/20 000/100
$F6 $	1E-5/1E-11/12/1/13	0.6/1E-4/1E-3/1E-2	1E-3/-1/5	1 000/2 000/60
$F7 $	1E-15/1E-16/19/0.2/10	1/1E-3/1E-4/1E-3	1E-3/-1/4	190 000/40 000/100
$F8 $	1E-5/1E-24/35/1/30	6E-2/1E-3/1E-4/6E-3	1E-4/-1/3	30 000/4 000/100
$F9$	1E-3/1E-18/25/1/25	3E-2/1E-4/1E-4/1E-3	1E-3/-1/5	8 000/5 000/60
$F10$	1E-9/1E-21/35/1/30	0.6/0.3/1E-4/0.8	1E-4/-1/5	50 000/250 000/130
$F11$	1E-9/1E-20/25/1/25	3E-3/1E-4/1E-3/1E-3	1E-4/-1/5	20 000/18 000/100
$F12 $	1E-6/1E-4/10/1/14	6E-4/1E-3/1E-5/0.1	1E-4/-1/5	300/1 000/60
$F13$	1E-15/1E-10/20/1/11	1/1E-3/1E-4/1E-2	1E-4/-1/10	80 000/6 000/60
$F14 $	1E-4/1E-18/35/0.2/15	1.5/0.1/1E-5/0.1	1E-4/-1/3	12 000/20 000/100
$F15 $	1E-4/1E-18/35/0.05/18	1.4/0.01/1E-5/0.1	1E-4/0/3	10 000/20 000/100

下载: 导出CSV 

| 显示表格

模糊子系统个数$m=5$, 对不同的例子$N_\sigma$可以不同, 如$F1, $ $N_\sigma=6$, 而$F2, $ $N_\sigma=5$.当$a=1, b=0, x=6$, 变换函数$T(a, b, x)=6/\sqrt{1+6^2}\approx0.99$, 这表明当$|x|>6$, $T(a_{\sigma l}, b_{\sigma l}, f)$的变化非常小, 为了防止变换函数重叠, 变换函数的参数$b_{\sigma l}$的初值由下式给出.

采用本文提出的基于变换函数与填充函数的模糊粒子群优化算法对上述15个基准函数进行测试, 按算法步骤在Matlab中编程仿真, 每个函数运行30次, 统计结果如下: 表 3给出了每个函数运行30次中的平均值, 最好寻优结果, 最坏寻优结果, 95$\%$置信区间, 成功率和运行时间.这里成功率定义为$N/30$, 其中$N$是运行结果满足$|f_{i}-f_{\min}| < SA$的次数, $f_{i}$为第$i$次运行结果, $f_{\min}$为被测函数的实际最优值.从中可以看出运行每一次都能搜索到最优值, 运行30次的成功率为$100\, \%$; 寻优结果精度高, 误差小, 并且除F11外平均寻优时间短.体现了FPSO-TF算法的高效性.

表 3  仿真结果

Table 3  Results of simulation

Functions	Average value	Best value	Worst value	Confidence interval	N/30	CPU times (s)	SA
$F1 $	0	0	0	0±0	30/30	12.7	1E-15
$F2$	2.37E-16	0	3.55E-15	2.37E-16±2.90724E-31	30/30	58.2	1E-10
$F3 $	0	0	0	0±0	30/30	42.7	1E-15
$F4$	1.57E-32	1.57E-32	1.60E-32	1.57E-32±1.5E-69	30/30	119	1E-15
$F5$	1.35E-32	1.35E-32	1.35E-32	1.35E-32±2.5E-95	30/30	118	1E-15
$F6$	-78.33233141	-78.33233141	-78.33233141	-78.33233141±1.9E-24	30/30	7.3	1E-7
$F7 $	3.23E-22	2.61E-26	4.39E-21	3.23E-22±3.66E-43	30/30	205	1E-15
$F8 $	1.21E-78	2.91E-221	3.60E-77	1.21E-78±1.54E-155	30/30	42.8	1E-15
$F9$	1.64E-52	1.94E-53	4.59E-52	1.64E-52±5.11E-105	30/30	18.2	1E-15
$F10$	2.74E-22	3.35E-32	1.11E-21	2.74E-22±5.12E-44	30/30	219	1E-15
$F11$	2.75E-07	0	7.87E-06	2.75E-07±7.39168E-13	30/30	1 816	1E-5
$F12$	3.82E-04	3.82E-04	3.82E-04	3.82E-04±2.71969E-38	30/30	4.1	1E-8
$F13$	6.36E-04	3.03E-05	2.54E-03	6.36E-04±1.7E-07	30/30	30.9	1E-2
$F14$	1.57E-32	1.57E-32	1.6E-32	1.57E-32±1.55E-69	30/30	126	1E-15
$F15$	1.35E-32	1.35E-32	1.35E-32	1.35E-32±2.5E-95	30/30	119	1E-15

下载: 导出CSV 

| 显示表格

5.1   更高维函数测试

当函数的维数增加时, 也就增加了问题的复杂性.为了进一步考察FPSO-TF算法的更高维寻优能力, 使函数维数从30增至60, 观察该算法的寻优收敛鲁棒性.设置测试函数$F1$\, $\sim$\, $F15$的维数分别为$n=30$, $n=40$, $n=50$, $n=60$, 且参数都使用维数$n=30$的.每个函数的寻优收敛曲线如图 11所示, 对于测试函数$F1$\, $\sim$\, $F8$, $F13$\, $\sim$\, $F15$随着维数从30增至60, 其寻优曲线依然能收敛至足够精度的极值点, 验证了该算法解决某些更高维多峰寻优问题的鲁棒性; 虽然其他测试函数对于更高维没有表现出很好的收敛精度, 但是对于30维均能收敛至高精度的极值点.这表明该算法很好地改善了传统粒子群算法的早熟收敛问题, 同时验证了该算法的高效性.测试函数$F$7, $F$13的寻优曲线在收敛过程中均出现阶梯型平台, 说明寻优过程陷入局部极小值点, 此时体现了填充函数的有效性, 使计算过程跳出局部极小值点继续寻优迭代至更小的极值点.并且大部分测试函数均能在较短时间内收敛至高精度的极值点, 体现了提出的多回路模糊控制系统全局搜索和局部搜索能力.

图 11  (a) $F$1, (b) $F$2, (c) $F$3, (d) $F$4, (e) $F$5, (f) $F$6 (这里曲线在每次迭代时都加78.33233140745), (g) $F$7, (h) $F$8, (i) $F$9, (j) $F$10, (k) $F$11, (l) $F$12 (这里曲线在每次迭代时都加0.000381827), (m) $F$13, (n) $F$14和(o) $F$15的收敛曲线

Fig. 11  Convergence progress of the FPSO-TF on (a) $F$1, (b) $F$2, (c) $F$3, (d) $F$4, (e) $F$5, (f) $F$6 (where curves are obtained by subtracting 78.33233140745 from the true value of $F$6 for each iteration), (g) $F$7, (h) $F$8, (i) $F$9, (j) $F$10, (k) $F$11, (l) $F$12 (where curves are obtained by subtracting 0.000381827 from the true value of $F$12 for each iteration), (m) $F13$, (n) $F14$ and (o) $F15$

下载: 全尺寸图片 幻灯片

5.2   均值与方差

表 4给出了在一定精度内FPSO-TF算法与各种比较算法的均值与方差.这些指标反映了算法的优化能力强弱.其中对比的算法来源于文献[5]的各种PSO算法, 参与对比的测试函数均采用30维, 运行25次得到均值与方差.

表 4  与现有算法的结果比较

Table 4  Comparison with other algorithms

Functions		FPSO-TF	FEP	OGA/Q	CMA-ES	JADE	OLPSO-L	OLPSO-G
$F1$	Mean	0	4.6E-2	0	1.76E+2	0	0	1.07
	SD	0	1.2E-2	0	13.89	0	0	0.99
	Rank	1	2	1	4	1	1	3
	t-test	-	-19.17	0	-63.36	0	0	-8.5
$F2$	Mean	2.37E-16	1.8E-2	4.4E-16	12.124	4.4E-15	4.14E-15	7.98E-015
	SD	2.90724E-31	2.1E-3	3.99E-17	9.28	0	0	2.03E-15
	Rank	1	6	2	7	4	3	5
	t-test	-	-42.8	-5	-6.5	-7.2E+16	-6.5+16	-19
$F3$	Mean	0	1.6E-2	0	9.59E-16	2.E-4	0	4.8E-3
	SD	0	2.2E-2	0	3.5E-16	1.4E-3	0	8.63E-3
	Rank	1	5	1	2	3	1	4
	t-test	-	-3.64	0	-13.7	-0.7	0	-2.7
$F7$	Mean	3.23E-22	5.06	0.75	2.33E-15	0.32	1.26	21.52
	SD	3.66E-43	5.87	0.11	7.7E-16	1.1	1.4	29.92
	Rank	1	6	4	2	3	5	7
	t-test	-	-4.26	-34.1	-15.1	-1.45	-4.5	-1.2
$F8$	Mean	1.21E-78	5.7E-4	0	4.56E-16	1.3E-54	1.11E-38	4.1E-54
	SD	1.54E-155	1.3E-4	0	1.13E-16	9.2E-54	1.3E-38	6.32E-54
	Rank	2	7	1	5	3	6	4
	t-test	-	-21.9	0	-20.4	-0.7	--4.3	-3.2
$F12$	Mean	3.82E-4	14.98	3.03E-2	3.15E+3	7.1	3.82E-4	3.84E+2
	SD	2.71969E-38	52.6	6.45E-4	5.79E+2	28	0	2.17E+2
	Rank	1	4	2	6	3	1	5
	t-test	-	-1.4	-234.8	-27.2	-1.3	0	-8.8
$F13$	Mean	6.36E-4	7.6E-3	6.3E-3	5.92E-2	6.8E-4	1.64E-2	1.16E-2
	SD	1.7E-7	2.6E-3	4.07E-4	1.73E-2	2.5E-4	3.25E-3	4.1E-3
	Rank	1	4	3	7	2	6	5
	t-test	-	-13.5	-71.3	-16.9	-0.88	-26.3	-12.8
$F14$	Mean	1.57E-32	9.2E-6	6.02E-6	1.63E-15	1.6E-32	1.57E-32	1.57E-32
	SD	1.55E-69	3.6E-6	1.16E-6	4.93E-16	5.5E-48	2.79E-48	1.01E-33
	Rank	1	5	4	3	2	1	1
	t-test	-	-12.8	-1.65	-16	-2.7E+18	0	0
$F15$	Mean	1.35E-32	1.6E-4	1.87E-4	1.71E-15	1.4E-32	1.35E-32	4.39E-4
	SD	2.5E-95	7.3E-5	2.62E-5	3.7E-16	1.1E-47	5.6E-48	2.2E-3
	Rank	1	4	5	3	2	1	6
	t-test	-	-10	-35.6	-23	-0.23E+15	0	-0.99
	Ave.rank	1.1	4.7	2.5	4	2.6	4.9	5.6
	Final rank	1	5	2	4	3	6	7

下载: 导出CSV 

| 显示表格

从表 4可以看出, 虽然仿真设定子系统个数为5个远比文献[5]所用的少, 但对于$F$2、$F$7, $F$13用FPSO-TF算法获得的平均值与误差均比其他比较算法更小, 体现了该算法的高效性; 对于$F1$、$F3$、$F12$、$F14$、$F$15用FPSO-TF算法获得的平均值等于其他比较算法中的最优结果, 误差与之相等甚至更小; 只有$F$8的平均值比OGA/Q的差, 仍然比其他算法结果的精度高.

表 4也对这些算法按均值排序, 最终顺序为: FPSO-TF、OGA/Q、JADE CMA-ES、FEP、OLPSO-L、OLPSO-G.同时表 4也给出了显著性水平为0.05, FPSO-TF和比较算法之间的t-test计算结果, 当t-test为负时, 意味着在均值和方差方面, FPSO-TF好于其他比较算法; 反之亦然.

6.   结论

本文结合粒子群优化算法与模糊控制系统, 并采用变换函数与填充函数, 提出了基于变换函数与填充函数的模糊粒子群优化算法(FPSO-TF).该算法通过设置$m$个不同宽度区间的的隶属度函数生成了$m$个模糊控制系统, 从而控制粒子的搜索范围, 使粒子更好地进行全局搜索和局部搜索, 提高搜索效率与精度.同时, 填充函数的引入使算法陷入局部极值点问题得到解决.最后的仿真结果与现有算法的比较结果表明该算法的有效性.